Lagrange-mesh solution of the Schrödinger equation in generalized spherical coordinates

The solution of the multi-dimensional Schrödinger equation in the generalized spherical coordinates is constructed in the Lagrange-mesh method. Laguerre and Jacobi meshes are used to construct matrix elements for the generalized Hamiltonian. The matrix elements are functions of quadrature abscissas and involve few free parameters for each angular dimension. A ring-shaped non-central separable potential, and a system of four linearly coupled anharmonic oscillators are used for illustrations of the efficiency and accuracy of the method. The numerical solutions involving eigenfunctions of the kinetic energy operator display the typical slow convergence with the accuracy improving as the Lagrange-mesh bases size increases. The Lagrange-mesh solutions converge to the exact solution when the parameters of the matrix elements are chosen appropriately.


Introduction
Theoretical studies of quantum mechanical systems rely on solutions to quantum mechanical equations for the systems to extract theoretical predictions of properties of the systems. Several numerical methods have been, and continue to be, developed to generate accurate solutions to the Schrödinger equation efficiently, see [1][2][3][4] for example. However, many of the methods become challenged when the dimensionality of the problem increases. The Lagrange-mesh method [5][6][7] was shown to have several advantages over many of the methods. This method is a variational method that employs the Lagrange basis functions defined on numerical quadratures to construct variational solutions to quantum mechanical equations. Compared to other numerical methods, the method was shown to generate simpler, faster and more accurate solutions for different quantum mechanical systems [8][9][10][11], and is not restricted to one-dimensional or separable problems. In the present work the Lagrange-mesh method is explored for the numerical solution of the Schrödinger equation in more than three dimensions in generalized spherical coordinates. Studies of this nature will benefit ongoing endeavors to find exact solutions to the many-body quantum mechanical problem [4,12], because of the striking similarity in the Hamiltonian of the problems.
The numerical solution of the Schrödinger equation in generalized spherical coordinates in D-dimensions (D>3) has, for a long time, been considered not practical [13]. However, the advent of the Lagrange-mesh method promises to render such a problem tractable. The Lagrange-mesh treatment of the Schrödinger equation in spherical coordinates in one-, two-and three-dimensions has been addressed in the literature [6,14,15]. However, the D-dimensional case is often reduced to a one-dimensional problem [8,16]. Lagrangemesh matrix elements based on the regularized Lagrange-Jacobi functions [7] are one of the recent developments of the method. These matrix elements extend the applicability of the Lagrange-mesh method in solving multi-dimensional quantum mechanical problems in polar coordinates to beyond the Fourier and the Legendre meshes. This paper investigates the accuracy and efficiency of the matrix elements by solving the Schrödinger equation in several dimensions, in spherical coordinates without directly using the usual spherical harmonics. The regularized Lagrange-Jacobi functions employed and, therefore, the associated kinetic energy matrix elements are characterized by free parameters that can be tuned to accelerate the convergence of the solution. In this work, the free parameters are chosen consistent with the eigenfunctions of the angular kinetic energy operator. The discussion of an effective procedure for tuning the free parameters is postponed to beyond this paper. Until now, these matrix elements had not been truly tested for efficiency and accuracy in practical applications.
This paper is organized as follows. Section 2 presents an outline of the D-dimensional Schrödinger equation in spherical coordinates. Principles of the Lagrange-mesh method and the associated matrix elements of the Hamiltonian for the Schrödinger equation are discussed in section 3. Illustrative applications of the matrix elements are discussed in section 4. The central harmonic oscillator and the Coulomb potential in six dimensions, four linearly coupled anharmonic oscillators, and a double ring-shaped non-central separable potential are considered for this purpose. Conclusions are given in section 5.

The Schrödinger equation in many dimensions
The Schrödinger equation in several dimensions is widely discussed in the literature, see for example [17,18]. This equation, for a particle of mass m in a potential field V r  ( ), in a D-dimensional configuration space, has the form (ÿ=m=1) where E is the energy, r Y  ( ) the wave function, and r  the position vector in D dimensions of the particle. This vector is characterized by Cartesian components {x 1 , x 2 , K, x D } and generalized spherical components {r, θ 1 , θ 2 , K, θ D−1 }, which are related by [17] x r x r x r x r x r x r for D2. The generalized spherical coordinates are defined in the domain r 0, (1) is solved for E and r Y  ( ) in a domain determined mainly by the potential.
In the Cartesian coordinates and in the spherical coordinates, the Laplace operator has the form [17] x r D r r L r is the generalized or total angular operator and , , , This angular operator is written in the compact iterative form [17] The iteration procedure, over k, say, is carried out using the relation . The operator (4) has known eigensolutions which, as shown later, are readily constructed. Based on the structure of the operators in (3) with (4), it may be argued that the Laplacian in Cartesian coordinates appears less involved, and seem easier to treat numerically, than in spherical coordinates.
The objective of this work is to solve (1) numerically in spherical coordinates. Bound-state solutions to this equation for a given potential V r  ( ) are determined by the Dirichlet boundary conditions as well as the squareintegrability condition where Ψ * denotes the complex conjugate of Ψ and τ the volume of integration with as the volume element. The presence of the first order derivatives in the Laplacian (3) introduce non-unique non-zero initial conditions on the wave function. To eliminate such first order derivatives in the wave function, the factorization is introduced, where the auxiliary function r F  ( ) satisfies the Dirichlet boundary conditions. This auxiliary function also satisfies the differential equation where the corresponding angular operator is constructed with . Note that the constant terms in (10) at each iteration step for k<D−1 have the form k 2 4 . The same finding holds at the final step k=D−1 only when the total Laplacian is considered. The eigensolutions of (10) are readily constructed in terms of the Jacobi functions which are more advantageous when considering (9) with a general potential V r  ( ). The eigenfunctions of the Laplace operator are readily developed from the iterative structure of the operator. Note that the operator is separable in all the coordinates, which justifies a factorization of the auxiliary wave function in the form where R(r) is the radial function and Θ j (θ j ) the angular functions. As a result, the differential equations for these functions, developed from (9) with V r 0 =  ( ) , are separable and have the form where σ 2 k are the separation constants. These three equations, the solutions of which are known, together describe the motion of a free particle. For a particle in a central potential only the radial equation (14), and its solution, differ from the free-particle case.
Focusing on the angular equations, the differential equation for θ 1 is immediately solved to obtain 2 exp i ; 0, 1, 2, . 15 All the differential equations (13) have a form that is consistent with the Pöschl-Teller Hamiltonian [19]. Such a Hamiltonian admits exact solutions that involve Gegenbauer functions, which are a special case of the Jacobi functions. For this reason, and noting that (13) have singularities at θ k =0 and θ k =π, the solutions to (13) can be shown to have the form are Gegenbauer polynomials of order ℓ k . It, then, follows that the corresponding eigenvalues are given by which defines the separation constants. It can be inferred from this expression that the separation constants obey the recurrence relation and , ] [17]. Note that C P cos cos where the latter are Jacobi polynomials.
Other solutions to (13) can be derived by considering the symmetry of the equations about the centers of the angular domains. Because of this symmetry, the equations may be solved considering only one segment of each domain and then adapt the results to the entire domain. Such a treatment is also necessary when solving (13) with potentials that have singularities that segment the angular domains [19,20]. This approach leads to general solutions of the form P sin cos cos 2 19 where c k are constants. These constants are determined consistent with the general Pöschl-Teller Hamiltonian [19]. The eigenvalues corresponding to (19) are ] . Therefore, solutions to (1) with general potentials may involve combinations of (16) and (19). However, the relation (16) may be more suited for domains in which the potential is even. Note that the kinetic energy operator is an even function of all the angles θ i (i>1).
Illustrative applications of the two solutions discussed above can be found in the literature. The D=3 case, the solution of which involve Θ 2 (θ 2 )Θ 1 (θ 1 ) usually referred to as the ring-shaped functions, is discussed in [21], while several D3 cases can be found in [19,20]. The construction of the numerical solutions Θ k (θ k ) to (13) in the Lagrange-mesh method is discussed in the next section. For this purpose, it should be noted that when

The Lagrange-mesh matrix elements
Lagrange-mesh bases functions are constructed so as to generate simple forms of matrix elements for common mathematical operators. A set of N Lagrange functions f i (z) are defined on a mesh z j ( j=1, 2, 3, K, N) in a domain [a, b] by the properties [6,7] f z V z f z V z 22 where the Dirac notation denotes integration (here Gauss approximated), V(z) a scalar function and f (k) =d k f/dz k . For a quadrature rule defined by a weight function w(z), the mesh is given by the quadrature abscissas and weights {z i , w i } and λ i =w i /w(z i ). When the domain [a, b] is scaled by some factor h, the properties will also depend on h. These properties hold for all Lagrange functions [2]. . The matrix elements of the operator −d 2 /dr 2 at the Gauss approximation have the exact form as those given in [6] for a regularization by r 3/2 . The Jacobi mesh is defined by N zeros z j of the Jacobi polynomials P α, β N (z) of order N and characterized by the weight function w z z z . This mesh is then used to construct regularized Lagrange-Jacobi functions [7] U z where μ is the regularizing parameter, h α, β N the normalizing constant of the Jacobi polynomials and . These functions are defined in the interval [−1, +1]. Other regularized Lagrange functions defined on meshes based on other classical orthogonal polynomials that are related to the Jacobi polynomials can be deduced from (26). The matrix elements of the Jacobi differential operator, constructed with these functions, at Gauss approximation, are given in [7].
To use the Lagrange-Jacobi function (26) in the construction of the solution to (10), the angular domains are first mapped onto the domains u 1, This can be achieved through a number of transformations, such as u n cos where n k are constants related to the periodicity of the potential. The exact form of the transformations, or values of the n k s, are determined by convenience, related mainly to the angular domain of the potential to be mapped. To simplify the discussion, without loss of generality, the transformations u n cos are considered, and the more simple cases of n k =1 and n k =2 used. For these cases, one readily shows that are independent of k. This paper focuses on the operator (29) because of its relative simplicity. The operator (30) is related to the Chebyshev polynomials of the first kind, which are a special case of the classical Jacobi polynomials [22]. Note that the transformation u cos q = necessarily involve (16) and (18), whereas u cos 2q = requires the use of (19) and (21).
Note that the differential operators in (29) are likely to have the same form at each iteration step. Therefore, the matrix elements of the operator (29) would involve similar, if not the same, u ˆwhich, as a result, may be calculated only once. The matrix elements ij u  for the operator (30), constructed with the regularized Lagrange-Jacobi functions U j (z), are given in [7]. Using any of the classical orthogonal polynomials related to the Jacobi polynomials as the Lagrange functions requires only the appropriate choice of the values for (α, β) in (26). These matrix elements are then used to construct iteratively, from k=1 to k=D−1, the matrix elements of the operator (29).
The solution to (9) is constructed by expanding the reduced wave function r u , for the parameters, where the abbreviation j j j j , The symbol ¢ signifies the exclusion, from the product, of the delta function corresponding to the iteration step index. The elements jj u  ¢ are given in [7] while T j j r 0 0 ¢ are given by a suitable Lagrange-Laguerre mesh [6]. It can be said that the elements H j j j j , 0 0 0 ¢ ¢   are simple and easy to construct. It must be emphasized that the elements H j j j j , )are the Gauss quadrature approximations, respectively, of the kinetic and potential energy matrix elements. These approximations introduce inherent integration errors in the numerical results [23]. Instructive illustrations of the use and accuracy of the matrix elements (32) are presented in the next section.

Illustrations
The following illustrations of the application of the Lagrange-mesh matrix elements involve the instructive Coulomb and harmonic oscillator potentials. Only the Gauss approximations of the matrix elements for the Hamiltonian of the problems are used in all the calculations presented. All the Gaussian quadrature abscissas and weights used are determined with accuracy no better than 3×10 −14 . It is noted that the matrix elements involve free parameters, three at each iteration step in the construction of the total Hamiltonian matrix. The values of these parameters are determined as follows. A closer look at (26), and comparing with (19), suggests that all the regularization parameters μ k for the angular functions may be set to  (15). However, the parameters (α k , β k ) may be optimized to accelerate the convergence of the bases expansions.

Central harmonic oscillator and Coulomb potentials
Simple, yet instructive, applications of the matrix elements presented are provided by the central harmonic oscillator and Coulomb potentials. The Schrödinger equations with these potentials can be solved exactly [20]. The central harmonic oscillator (ho) and Coulomb (C) potentials have the forms where n i =0, 1, 2, Kfor all i and σ D−1 is the parameter for the associated Laguerre polynomials which are part of the analytical solution to the central Coulomb problem [20]. This parameter coincides with that of the Laguerre mesh for this problem. The Schrödinger equations for these potentials are solved in six dimensions, as a one-dimensional radial problem considering only (14), and directly as a six-dimensional problem using (32). =´where N 0 =3 for the oscillator and N 0 =30 for the Coulomb potential. The radial domain for the Coulomb potential was scaled by a parameter h=1/(σ D−1 +1) to focus on the ground state [6], while the domain for the oscillator was unscaled (h=1). The same basis size N 0 for the radial functions was used in both the one-dimensional and the six-dimensional treatments.
It can be seen, in the table, that the energies for the two potentials are reproduced to machine accuracy in the one-dimensional treatment. In the six-dimensional treatment the results are reproduced to close to machine accuracy with the degeneracy of the states correctly determined. In general, the accuracy is the same for all the calculated states and appears to be independent of the bases size in the case of the oscillator. The maximum error in the results of the direct six-dimensional treatment is ò n ∼10 −15 in the case of the oscillator potential, which is an order of magnitude greater that the error in the one-dimensional treatment. The error difference in the two treatments is two orders of magnitude in the case of the Coulomb potential. However, note that in the Coulomb case a larger radial basis was required to generate results of reasonable accuracy. This difference in the errors can be attributed to the use of the matrix elements ij u  in the solution of the problems. The error differences are, however, very small considering the dimensional size of the problem. The results of the Coulomb potential display a similar trend in accuracy as the oscillator potential except that, in this case, the accuracy deteriorates as the level of the states increases. Errors in the results of problems involving only central potentials are expected to display a similar trend.

Ring-shaped potential
An interesting application is provided by a more practical problem of a non-central ring-shaped potential coupled to a central Coulomb potential in three-dimensions (D = 3). The potential has the form [ where a i are constants. This potential leads to a separable Schrödinger equation that admits exact solutions with eigenvalues The regularized Laguerre and the Jacobi meshes were used to solve the Schrödinger equation for this potential. The potential parameters were set to a 0 =a 1 =2 and a a a b = ( ) ( ) are tested for convergence. The former set of values corresponds to the case of central potentials while the latter corresponds to a symmetric Pöschl-Teller in u 2 . The radial domain was treated as in the case of the central Coulomb potential, where the radial basis was set to N 0 =30, however, here the scaling parameter is set to h=1. Since the basis size corresponding to u 1 could also be fixed, only the u 2 basis size N 2 was varied. The convergence of the energies, with increasing N 2 , for the first two states of the potential are shown in table 2. As can be seen, in the table, the accuracy of the calculated energies change from relative errors of ∼10 −4 involving N 2 =10 basis functions to ∼10 −6 with N 2 =40 basis functions. This signifies slow convergence of the numerical solutions, wich is characteristic of the eigenfunction expansion approach [25], determined with both sets of parameters. It also appears that the results for one set of the trial values of (α 2 , β 2 ) converge from above while those of the other converge from below, though at comparable rates, to the exact energies. This trend can be used to tune the trial values of (α 2 , β 2 ) to improve convergence. Recalling that the parameter σ D−1 depends on (α 2 , β 2 ), the variation of the radial basis size N 0 may affect the accuracy of the results.

Linearly coupled anharmonic oscillators
Another interesting application is offered by the non-central and non-separable system of four linearly coupled anharmonic oscillators. The potential for this system in the generalized spherical coordinates has the form . The Schrödinger equation for this potential is non-separable but has exact solutions with eigenvalues [26] where u u u u 1 a n d 1 .
The modified Laguerre and the Jacobi meshes were used to solve the Schrödinger equation for this potential, and the radial domain was unscaled (h = 1). The potential does not have symmetries that are as readily noticeable as in Cartesian coordinates. As a result, optimal values of all the parameters (α k , β k ) and σ D−1 cannot be easily identified for trial. For this reason, the parameters are determined through trial-and-error. It is observed, in the table, that the energies of the first two levels, calculated using parameter set (a), are obtained with relative errors of ∼10 −3 with N i =6 by the current method, which are similar to the accuracy for the Cartesian treatment. A slight change in the parameter values, set (b), improved the accuracy of the calculated energy of the ground state energy by two orders of magnitude. Increasing the bases sizes to N i =10 does not improve the accuracy significantly. However, the results for the Cartesian treatment improved by four orders of magnitude. It is seen that the calculated energy E 2 =3.346 of second excited state is not successfully located and, instead, the nearby third excited state E 3 =3.565 is located, in the present treatment. Nonetheless, the rate of convergence with both bases sets is a lot slower than in the Cartesian-coordinate calculations. This is not surprising since, unlike in the present work, the results of [26] were generated with optimized bases. The rate of convergence of the results of this work is, therefore, expected to improve when the parameters (α k , β k ) and σ D−1 are optimized. A systematic optimization procedure for this purpose still needs to be explored.

Conclusions
The Lagrange-mesh method was used to construct numerical solutions to the multi-dimensional Schrödinger equation in spherical coordinates. The Laguerre mesh was used to treat the radial domain while Jacobi meshes were used to treat all the parametrized angular domains. Standard Laguerre-mesh matrix elements were used for the radial component, while the newly developed Jacobi-mesh matrix elements were used for the parametrized angular components of the kinetic energy operator. Simple, yet instructive, applications involving central, noncentral separable, and non-central non-separable Coulomb as well as harmonic oscillator-type potentials in three, four and six dimensions were treated. Results for the Lagrange-mesh solutions of the two central potentials problems in six dimensions were reproduced within machine accuracy. These results showed practically no error contributions from the direct treatments involving the use of the parametrized angular kinetic energy matrix elements. Results for the solutions of the non-central potentials displayed the characteristic slow convergence of the eigenfunction expansion approach when the matrix elements corresponding to the eigenfunctions of the kinetic energy operator were used. However, the Lagrange-mesh results with smaller bases sizes for the anharmonic oscillator potential in four dimensions had the same accuracy as those of a Cartesian treatment of the same potential. The results also showed that Lagrange-mesh matrix elements with optimized free parameters may improve the convergence rate. The numerical solution of the Schrödinger equation in generalized spherical coordinates, in more than three dimensions, is realized through the Lagrange-mesh method. The matrix elements of the Hamiltonian for the problem are simple and require no computer time to evaluate. The results of this work has demonstrated that the multi-dimensional Schrödinger equation in spherical coordinates can be numerically solved with the same simplicity as in Cartesian coordinates. However, a more rigorous optimization procedure for determining the free parameters, of the Lagrange-mesh matrix elements associated with the kinetic energy operators, needs to be investigated. Currently such parameters are optimized empirically through trial-and-error. The independence of the potential matrix elements from the free parameters renders many standard parameter optimization techniques inapplicable.