Numerical Solutions of Quantum Mechanical Eigenvalue Problems

1. Introduction

The Schrödinger equation has been a central part of “modern” physics for almost a century. When interpreted broadly, it can be formulated in a multitude of ways [1]. Here we mainly restrict our discussion to the non-relativistic, time independent form,

$\begin{array}{ll}\left[-{\Delta }_{\text{q}}+V(\text{q})\right]\psi (\text{q})=E\psi (\text{q}).& (1)\end{array}$

This constitutes an eigenvalue problem for E (there are many cases where the operator defined by Equation (1) allows for a continuous spectrum of E-values, but this will not directly influence the treatment of finite discretizations of such systems). In Equation (1), q denotes the configuration space coordinate for a system of one or more particles in one or more spatial dimensions, and Δ_q is a Laplace operator on this configuration space. V(q) is the interaction potential, and E the eigenvalue parameter, interpreted as an allowed energy for the quantum system.

Despite its appearance as a single-particle equation, Equation (1) can also be used to model N-particle systems, with q = (r₁, …, r_N) and Δ_q = (c₁Δ₁, …, c_NΔ_N). Here each Δ_k is an ordinary flat space Laplace operator, and c_k is a numerical coefficient inversely proportional to the mass m_k of particle k; this mass may differ from particle to particle. By a suitable scaling of each coordinate r_k, one can mathematically transform all c_k to (for instance) unity. But such transformations may obscure physical interpretations of the coordinates, and make mathematical formulations more error-prone.

How to solve eigenvalue problems like (1)? Fortunately for the rapid initial development of quantum mechanics, for many important physical cases [like the hydrogen atom [2, 3] and harmonic oscillators [4]] it could be reduced to a set of one-dimensional eigenvalue problems, through the separation of variables method. Moreover, the resulting one-dimensional problems could all be solved exactly by analytic methods. The origins for such fortunate states of affairs can invariably be traced to an enhanced set of symmetries. However, not every system of physical interest enjoy a high degree of symmetry. Even most one-dimensional problems of the form (1) have no known analytic solution. A popular and much investigated system is the anharmonic oscillator,

$\begin{array}{ll}\left[-\frac{d^2}{d{x}^2}+\mu x^2+\epsilon x^4\right]\psi (x)=E\psi (x).& (2)\end{array}$

This model has often functioned as a theoretical laboratory [5, 6], for instance to investigate the behavior and properties of perturbative [7, 8] and other [9–12] expansions, and alternative solution methods [13–15].

It this article we describe some attempts to simplify numerical solutions of eigenvalue problems like (1). Our approach relies on standard numerical algorithms, already coded and freely available through Python packages like numpy [16] and scipy [17, 18]. The main aim is to automatize the transformation of (1) to function calls accepted by the numerical eigenvalue solvers. Within the above class of models, the problem is completely defined by the coefficient vector (c₁, c₂, …, c_N) and the real function V(q). In principle, this should be the only user input required for a numerical solution.

In practice some additional decisions must be made, like how a possibly infinite configuration space should be reduced to a region of finite extent, how the boundaries of this region should be treated, and how this region should be further approximated by a finite lattice. Other options involve selection of numerical approaches, like whether dense or iterative sparse matrix solvers should be used. Such decisions have consequences for many “trivial” details of the numerical programs, but they can be provided in the form of parameters and selectors, automatically implemented without further tedious and error-prone human intervention. Even many of the decisions indicated above may ultimately by delegated to artificial intelligence systems, but this is beyond our current scope.

2. Available Python Procedures for Numerical Solution

Numerical approaches to problems like those above are in principle straightforward: The operator

$\begin{array}{l}\text{H}=\text{T}+\text{V}\end{array}$

defined by Equation (1) is approximated by a finite real symmetric matrix

$\begin{array}{l}{M}_{\text{H}}=M_{\text{T}}+M_{\text{V}}\end{array}$

where we have introduced the symbol T = −Δ_q. For densely defined matrices M_H there are several standard numerical eigenvalue solvers available, like eig and eigvals in the scipy.linalg package. A 10⁴ × 10⁴ matrix of double precision numbers requires 800 Mb of storage space; this is indicative of the problem magnitudes that can be handled by dense matrix methods on (for instance) modern laptops. That is, such computers have more than enough memory for numerical treatment of one-dimensional problems, and usually also sufficient memory for two-dimensional ones.

For higher-dimensional problems one may utilize the sparse nature of M_H to find solutions through iterative procedures, like the eigsh eigenvalue solver in the scipy.sparse.linalg package. This solver does not require any explicit matrix construction of M_H, only a LinearOperator function that returns the vector M_Hψ for any input vector ψ. In the representations we consider, M_V is always diagonal, and M_T can be made diagonal by a Fast Fourier Transform (FFT), or some of its discrete trigonometric variants. This opens the possibility it to handle non-sparse matrix problems, where T is replaced by more general expressions of F(T), by the same procedures. For instance functions F that involves fractional and/or inverse powers of its arguments.

3. Required Parameters and Selectors

In this section we describe the additional quantities that a user must input for a full specification of the numerical problem. They assume that configuration space has been modeled by a rectangular point lattice, with a selection of possible boundary conditions.

3.1. Lattice Shape

The most basic quantity of the numerical model is the discrete lattice approximating the relevant region of configuration space. For rectangular approximations this is defined by the shape parameter, a Python tuple,

$\begin{array}{ll}\text{shape}=(s_0,s_1,\dots ,s_{d-1}),& (3)\end{array}$

where each s_k is a positive integer specifying the number of lattice points in the k'th direction, and d is the (effective) dimension of configuration space. For models with continuous symmetries (for instance rotational ones) the effective dimension may be chosen smaller than the physical one, by separation of variables. Likewise, discrete symmetries may can used to reduce the size of configuration space that this lattice must approximate.

In Python programs, quantities like the wave function ψ and the potential V are defined as floating point NumPy arrays of shape shape.

3.2. Edge Lengths and Offsets

The geometric extent of the selected region is specified by its edge lengths xe. This is a NumPy array of positive floating point numbers,

$\begin{array}{ll}\text{xe}=\left[e_0,e_1,\dots ,e_{d-1}\right].& (4)\end{array}$

A secondary quantity, derived from xe and shape is the elementary lattice cell,

$\begin{array}{ll}\text{dx}=\text{xe}/\text{shape}=\left[e_0/s_0,e_1/s_1,\dots ,e_{d-1}/s_{d-1}\right].& (5)\end{array}$

The absolute positioning of the region, with respect to some fixed coordinate system, is specified by a NumPy array of floating point numbers,

$\begin{array}{ll}\text{xo}=\left[x_0,x_1,\dots ,x_{d-1}\right].& (6)\end{array}$

This is defined as the position of the “lower left” corner of the selected region. The placement of the lattice points within the region still needs to be specified, as will be discussed below.

3.3. Boundary Conditions

The restriction to finite regions of space requires imposition of boundary conditions. For regions of rectangular shape (generalized to arbitrary dimensions), as considered here, the perhaps simplest choice is periodic boundary conditions in each directions. This may be viewed as a topological property of configuration space itself. Other boundary conditions are really properties of functions defined on this space, as specifications of how the functions should be extended beyond the boundary. Two natural choices are symmetric and anti-symmetric extensions. With a lattice approximation a further distinction can be made, related to how the lattice points are positioned relative to the boundary.

In this connection, it is natural to consider the cases handled by the trigonometric cousins of the fast Fourier transform (FFT). In the one-dimensional case the extension may be symmetric or anti-symmetric with respect to a boundary, which is situated either (i) at a lattice point, or (ii) midway between two lattice points. Thus, at each boundary there is 2 × 2 matrix of possibilities, as indicated by Table 1.

TABLE 1

Table 1. Individual boundary conditions covered by standard discrete trigonometric transforms (DCT and DST).

With two boundaries there are altogether 4 × 4 = 16 possibilities. However, the routines in scipy.fftpack (dct and dst of types I–IV) only implement cases where both options come from the same row of Table 1. With the periodic extension P in addition, one ends up with a set of nine possibilities in each direction:

$\begin{array}{ll}ℬ=\left\{\prime \text{P}{\text{P}}^{\prime }{,}^{\prime }\text{S}{\text{S}}^{\prime }{,}^{\prime }\text{S}{\text{A}}^{\prime }{,}^{\prime }\text{A}{\text{S}}^{\prime }{,}^{\prime }\text{A}{\text{A}}^{\prime }{,}^{\prime }\text{s}{\text{s}}^{\prime }{,}^{\prime }\text{s}{\text{a}}^{\prime }{,}^{\prime }\text{a}{\text{s}}^{\prime }{,}^{\prime }\text{a}{\text{a}}^{\prime }\right\}.& (7)\end{array}$

Hence, the numerical model must be further specified by a Python tuple of two-character strings, defining the selected boundary condition in all directions,

$\begin{array}{ll}\text{bc}=(b_0,b_1,\dots ,b_{d-1})& (8)\end{array}$

with each $b_k\in {B}$ (or in an enlarged set of possibilities).

3.4. Lattice Positions. Dual Lattice Squared Positions

When bc is given, one may automatically calculate the positions of all lattice points

$\begin{array}{ll}\text{xlat}=({\text{X}}_0,{\text{X}}_1,\dots ,{\text{X}}_{d-1}),& (9)\end{array}$

provided shape, xe, and xo are also known. In Equation (9), the property xlat is a tuple of one-dimensional arrays. For illustration, consider the case of a 3-dimensional lattice of shape (s_x, s_y, s_z). Then xlat is a Python tuple (X, Y, Z), where X is a numpy array of shape (s_x, 1, 1), Y is a numpy array of shape (1, s_y, 1), and Z is a numpy array of shape (1, 1, s_z). These are all one-dimensional arrays, but their shape information implies that (for instance) the Python expression X*Y automatically evaluates to a numpy array of shape (s_x, s_y, 1).

A Python function V(x, y, z), defined by an expression that can involve “standard” functions, may then be evaluated on the complete lattice by the short and simple expression V(*xlat). When V depends on all its arguments, the result will be a numpy array of shape (s_x, s_y, s_z).

In general, when Fourier transforming a periodic function f(x), where x takes values on some discrete lattice, the result becomes another periodic function $\tilde{f}(\text{k})$, where k takes values on another discrete lattice (the dual lattice/reciprocal space). Modulo an overall scaling, a set of k-values (labeling the points of some complete, minimal subdomain to be extended by periodicity) can be defined such that f(x + a) transforms to ${\mathrm{\text{e}}}^{-\mathrm{\text{i}}\text{k}·\text{a}}\tilde{f}(\text{k})$. A natural choice for that minimal domain is, in physicists language, the first Brillouin zone (this choice may still leave a somewhat arbitrary selection of boundary points to be included). On this subdomain of the dual lattice, derivatives can be defined as the multiplication operators −ik. But these operators must still be extended to the full dual lattice by periodicity. The common stensil expressions for lattice derivatives correspond to the lowest Fourier components of the (periodically extended) multiplication operator −ik.

For the other (discrete trigonometric) transformations a complication arises, because a derivation also induces a transposition of the boundary conditions in ${B}$. However, two derivations in the same direction leave the boundary conditions unchanged, and hence can be represented as a multiplication operator q on the transformed functions. Let ∂_k be shorthand notation for ∂/∂x_k. The previous conclusion implies that all operators of the form $F({\partial }_0^2,{\partial }_1^2,\dots ,{\partial }_{d-1}^2)$ can be evaluated through multiplications and fast discrete transforms,

$\begin{array}{ll}F({\partial }_0^2,{\partial }_1^2,\dots ,{\partial }_{d-1}^2)={{T}}^{-1}F(q_0,q_1,\dots ,q_{d-1}){T}.& (10)\end{array}$

We have implemented code that performs ${T}$ and ${{T}}^{-1}$ through a sequence of discrete trigonometric or fast Fourier transforms, dependent on bc and the other parameters. Analogous to the arrays xlat of lattice positions (Equation 9), one may automatically calculate similar arrays of squared positions for reciprocal lattice,

$\begin{array}{ll}\text{qlat}=({\text{Q}}_0,{\text{Q}}_1,\dots ,{\text{Q}}_{d-1}).& (11)\end{array}$

3.5. Lattice Laplacian. Stensil Representations

Instead of relying on FFT type transforms, one may directly construct discrete approximations (stencils) of the Laplace operator, and similar differential operators. The simplest implementation of a lattice Laplacian in one dimension is obtained by use of the formula

$\begin{array}{ll}\frac{d^2\psi }{d{x}^2}(x_n)\approx \frac{\psi (x_n+\delta x)-2\psi (x_n)-\psi (x_n-\delta x)}{\delta x^2},& (12)\end{array}$

where δx is the distance between nearest-neighbor lattice points. The formal discretizations error of this approximation is of order δx². By summing such expression in d orthogonal directions one finds the (2d + 1)-stensil expression for the lattice Laplacian.

A more accurate approximation is the (4d + 1)-stencil,

$\begin{array}{ll}\Delta \phi ({\text{x}}_{\text{n}})\approx {\displaystyle \sum _{k=0}^{d-1}}\frac{-\phi ({\text{x}}_{\text{n}}+2{\delta }_k)+16\phi ({\text{x}}_{\text{n}}+{\delta }_k)-30\phi ({\text{x}}_{\text{n}})+16\phi ({\text{x}}_{\text{n}}-{\delta }_k)-\phi ({\text{x}}_{\text{n}}-2{\delta }_k)}{12|{\delta }_k{|}^2}.& (13)\end{array}$

Here δ_k denotes a vector of length |δ_k| pointing in positive k-direction.

An arbitrary (short-range) position independent operator O can in general be represented by a stensil s_O(b) such that

$\begin{array}{ll}(O\psi )({\text{x}}_{\text{n}})={\displaystyle \sum _{\text{b}}}{s}_O(\text{b})\psi (x_{n-b}).& (14)\end{array}$

When n − b falls outside the lattice, the value of ψ(x_{n − b}) must is interpreted according to the boundary conditions bc. This can again be automatized. We have implemented an algorithm for this, currently only for 5 of the 9 cases in ${B}$ in each direction, but for an arbitrary number of directions.

The various ways to approximate the Laplace operator, or more generally the kinetic energy operator, is made available through the selector ke, whose value is currently limited to the set of options {′2dplus1′,′4dplus1′,′fftk2′}. The last of these options is discussed in section 5.

4. Simple Applications

In this section we will demonstrate some applications of our automatic code. The main requirement is that in each case only a set of parameters and selectors should be provided, with no coding required by the application itself. This should be sufficient to generate eigenvalues E_n as requested, and optionally also the associated eigenfunctions (an issue which we have not yet tested).

4.1. Example: One-Dimensional Harmonic Oscillator

Consider the eigenvalue problem of the one-dimensional harmonic oscillator,

$\begin{array}{ll}-{\psi }_n^{″}(x)+x^2{\psi }_n(x)=E_n{\psi }_n(x).& (15)\end{array}$

The eigenvalues are E_n = 2n + 1 for n = 0, 1, …, and the extent of the wavefunction ψ_n(x) can be estimated from the requirement that a classical particle of energy E_n is restricted to $x^2\le E_n$. A quantum particle requires a little more space than the classically restricted one.

For a numerical analysis we provide the parameters

$\begin{array}{l}\text{shape}=(128,),\text{bc}={(}^{\prime }{\text{a}}^{\prime },{\text{a}}^{\prime }),\text{xe}=(25,),\text{xo}=(-12.5),\\ \text{V}=\text{lambda}x:x**2,\end{array}$

selects the 3-stensil approximation for T (default choice), and the dense matrix solver eigvalsh (default choice). This instantly returns 128 eigenvalues as plotted in Figure 1. We may easily change shape to (1024), for a much better result. The potential for additional explorations, without any coding whatsoever, should be obvious.

FIGURE 1

Figure 1. The 128 lowest eigenvalues of Equation (15), computed with the standard 3-stensil approximation for the Laplace operator (here the kinetic energy T). The parameters are chosen to illustrate two typical effects: With the bc = (a, a) boundary conditions the harmonic oscillator potential is effectively changed to V = ∞ for x ≥ 12.5, thereby modifying the behavior of extended (highly exited) states. The effect of this is to increase the eigenenergies of such states, to a behavior more similar to a particle-in-box. This is visible for n ≳ 80. The effect of using the 3-stensil approximation for T is to change the spectrum of this operator from k² to (the slower rising) (2/δx)²sin²(kδx/2). This is visible in the sub-linear rise of the spectrum for N = 2⁷.

For a better quantitative assessment of the accuracy obtained we plot some energy differences, $E_n^{(\mathrm{\text{exact}})}-E_n$, in Figure 2.

FIGURE 2

Figure 2. The discretizations error of energy eigenvalues when using the standard 3-stensil approximation for the one-dimensional Laplace operator (here the kinetic energy T). There is no improvement in E₉₀ beyond a certain lattice size N, because the corresponding oscillator state is too large for the geometric region. Hence, for improved accuracy of higher eigenvalues one should instead increase the xe, while maintaining xo = −xe/2. For the other states the improvement is consistent with the expectation of an error proportional to δx². This predicts an accuracy improvement of magnitude 2¹² = 4, 096 when the number of lattice sites increases from N = 2⁷ to N = 2¹³ for a fixed geometry. The eigenvalues are computed by the dense matrix routine eigvalsh from scipy.linalg.

This brute force method leads to a dramatic increase in memory requirement with increasing lattice size. For a lattice with N = 2^m sites, the matrix requires storage of 4^m double precision (8 byte) numbers. For m = 13 this corresponds to about $\frac{1}{2}\mathrm{\text{ Gb}}$ of memory, for m = 14 about 2Gb. The situation becomes even worse in higher dimensions.

Assuming that we are only interested some of the lowest eigenvalues, an alternative approach is to calculate these by the iterative routine eigsh from scipy.sparse.linalg. This allows extension to larger lattices, as shown in Figure 3.

FIGURE 3

Figure 3. The discretizations error computed by the routine eigsh from scipy.sparse.linalg. For a fixed lattice size the discretizations error is essentially the same as with dense matrix routines. However, with a memory requirement proportional to the lattice size (instead of its square) it becomes possible to go to much larger lattices. This figure also demonstrates (E₇₀) that the error can be limited by boundary effects instead of the finite discretization length δx.

With a sparse eigenvalue solver the calculation becomes limited by available computation time, which in many cases is a much weaker constraint: With proper planning and organization of calculations, the relevant timescale is the time to analyze and publish results (i.e., weeks or months). The computation time is nevertheless of interest (it shouldn't be years). We have measured the wall clock time used to perform the computations for Figures 2, 3, performed on a 2012 Mac Mini with 16 Gb of memory, and equipped with a parallelized scipy library. Hence, the eigvalsh and eigsh routines are running with four threads. The results are plotted in Figure 4.

FIGURE 4

Figure 4. The wall clock time used to find the lowest 128 eigenvalues, for various systems and methods. We have also used the dense matrix routine eigvalsh to compute the eigenvalues of a 2⁷ × 2⁷ (N = 2¹⁴) two-dimensional lattice; not unexpected it takes the same time as for a 2¹⁴ one-dimensional lattice. Somewhat surprisingly, with eigsh it is much faster to find the eigenvalues for two-dimensional lattice than for a one-dimensional with the same number of sites, and somewhat faster to find the eigenvalues for a three-dimensional lattice than for the two-dimensional with the same number of sites.

Here we have used the eigsh routine in the most straightforward manner, using default settings for most parameters. This means, in particular, that the initial vector for the iteration (and the subsequent set of trial vectors) may not be chosen in a optimal manner for our category of problems. It is interesting to observe that eigsh works better for higher-dimensional problems. The (brief) scipy documentation [17] says that the underlying routines works best when computing eigenvalues of largest magnitude, which are of no physical interest for our type of problems. It is our experience that the suggested strategy, of using the shift-invert mode instead, does not work right out-of-the-box for problems of interesting size (i.e., where dense solvers cannot be used). We were somewhat surprised to observe that the computation time may decrease if the number of computed eigenvalues increases (cf. Figure 5).

FIGURE 5

Figure 5. One may think that it takes longer to compute more eigenvalues. This is not always the case when the number of eigenvalues is small, as demonstrated by this figure. The default choice of eigsh is to compute k = 6 eigenvalues. For our two- and three-dimensional problems this looks close to the optimal value, but it is too low for the one-dimensional problem.

4.2. Example: 2- and 3-Dimensional Harmonic Oscillators

The d-dimensional harmonic oscillator

$\begin{array}{ll}\left[-\Delta +r^2\right]{\psi }_n(r)=E_n{\psi }_n(r),& (16)\end{array}$

has eigenvalues E_n = (d + 2n), for n = 0, 1, …. The degeneracy of the energy level E_n is g_n = (n + 1) in two dimensions, and $g_n=\frac{1}{2}(n+1)(n+2)$ in three dimensions¹. This degeneracy may be significantly broken by the numerical approximation. For a numerical solution we only have to change the previous parameters slightly:

$\begin{array}{l}\text{shape}=(128,)*\text{dim},\text{bc}=({(}^{\prime }{\text{a}}^{\prime }{,}^{\prime }{\text{a}}^{\prime }),)*\text{dim},\\ \text{xe}=(25,)*\text{dim},\text{xo}=(-12.5,)*\text{dim},& (\mathrm{17a})\end{array}$

$\begin{array}{ll}\text{V}=\text{lambda}x,y:x**2+y**2(\text{dim}=2),& (\mathrm{17b})\end{array}$

$\begin{array}{ll}\text{V}=\text{lambda}x,y,z:x**2+y**2+z**2(\text{dim}=3),& (\mathrm{17c})\end{array}$

for dim = 2, 3.

As already discussed, the routine eigsh works somewhat faster in higher dimensions than in one dimension (for the same total number N of lattice points). The corresponding discretizations errors are shown in Figures 6, 7.

FIGURE 6

Figure 6. The discretization error of energy eigenvalues when using the standard 5-stensil approximation for the two-dimensional Laplace operator. Exactly, the states E₇₈ and E₉₀ are the two edges of a 13-member multiplet with energy 26, and the state E₁₂ is the middle member of a 5-member multiplet with energy 10. With the chosen parameters all states considered a well confined inside the geometric region; hence we do not observe any boundary correction effects.

FIGURE 7

Figure 7. The discretization error of energy eigenvalues when using the standard 7-stensil approximation for the three-dimensional Laplace operator. Exactly, the states E₅₆ and E₈₃ are the two edges of a 28-member multiplet with energy 15, and the state E₁₅ is the middle member of a 10-member multiplet with energy 9.

The discretizations error continues to scale like δx². This means that a reduction of this error by a factor 2² = 4 requires an increase in the number of lattice points by a factor 2^d in d dimensions. This means that is becomes more urgent to use a better representation of the Laplace operator in higher dimensions. Fortunately, as we shall see in the next sections, better representations are available for our type of problems.

5. FFT Calculation of the Laplace Operator

One improvement is to use the reflection symmetry of each axis (x → −x, y → −y, etc.) to reduce the size of the spatial domain. This reduces δx by a half, without changing the number of lattice points.

A much more dramatic improvement is to use some variant of a Fast Fourier Transform (FFT): After a Fourier transformation, $\psi (\text{r})\to \tilde{\psi }(\text{k})$, the Laplace operator turns into multiplication operator,

$\begin{array}{l}(-\Delta \psi )(r)\to {\text{k}}^2\tilde{\psi }(\text{k}).\end{array}$

This means that application of the Laplace operator can be represented by (i) a Fourier transform, followed by (ii) multiplication by k², and finally (iii) an inverse Fourier transform. Essentially the same procedure works for the related trigonometric transforms.

For rectangular lattices, these options can also be implemented as practical procedures, due to the existence of efficient and accurate² algorithms for discrete Fourier and trigonometric transforms. The time to perform the above procedure is not very much longer than the corresponding stensil operations. The benefit is that the Laplace operator becomes exact on the space of functions which can be represented by the modes included in the discrete transform.

We have coded this FFT-type representation of the Laplace operator for the various types of boundary conditions listed in Table 1. This possibility can be chosen as an option for the kinetic energy selector, ke. The obtainable accuracy through this option increases dramatically, as illustrated in Figures 8–11. As shown in Figure 8, a decrease of the lattice length δx does not necessarily lead to a more accurate result. We attribute this to an enhanced amplification of roundoff errors.

FIGURE 8

Figure 8. With a FFT representation of the Laplace operator the discretization error drops exceptionally fast with δx∝N⁻¹. When it becomes “small enough” the effect of numerical roundoff becomes visible; the latter leads to an increase in error with δx. The results in this figure is for a one-dimensional lattice, but the behavior is the same in all dimensions. The lesson is that we should make δx “small enough” (which in general may be difficult to determine a priori), but not smaller. It may also be possible to rewrite the eigenvalue problem to a form with less amplification of roundoff errors.

FIGURE 9

Figure 9. Accuracy of computed eigenvalues for a 1D oscillator, using the FFT approximation for kinetic energy T. This figure may suggest that an increase in the number of lattice size N will lead to a accurate treatment of states with higher n. Our findings are that this is not the case: The results for N = 2⁷ and N = 2⁸ have essentially the same behavior as for N = 2⁶.

FIGURE 10

Figure 10. Accuracy of computed eigenvalues for a 2D oscillator, using the FFT approximation for kinetic energy T. As can be seen, a large number of the lowest eigenvalues can be computed to an absolute accuracy in the range 10⁻¹⁴–10⁻¹² with a lattice of size 2⁶ × 2⁶. We observe not improvement by going to 2⁷ × 2⁷ lattice, but no harm either (except for an increase in the wall clock execution time from about 3 to 30 s for each combination of boundary conditions).

FIGURE 11

Figure 11. Accuracy of computed eigenvalues for a 3D oscillator, using the FFT approximation for kinetic energy T. As can be seen, a large number of the lowest eigenvalues can be computed to an absolute accuracy in the range 10⁻¹⁴ to 10⁻¹² with lattice of size 2⁶ × 2⁶ × 2⁶. We observe no improvement by going to 2⁷ × 2⁷ × 2⁷ lattice, but no harm either (except for an increase in the wall clock execution time from about 6 to 95 min for each combination of boundary conditions).

It might be that the harmonic oscillator systems are particularly favorable for application of the FFT representation. One important feature is that the Fourier components of the harmonic oscillator wave functions vanishes exponentially fast, like e^−k2/2, with increasing wave-numbers k². This feature is shared with all eigenfunctions of polynomial potential Schrödinger equations, but usually with different powers of k in the exponent, which quantitatively leads to a somewhat different behavior. Furthermore, the onset of exponential decay will occur for larger values of k² for the more excited states (i.e., with larger eigenvalue numbers).

For systems with singular wavefunctions the corresponding Fourier components may vanish only algebraically with k². An equally dramatic increase in accuracy cannot be expected for such cases.

6. Anharmonic Oscillators

Our general setup allows for any computable potential, by simply changing the definition of the function assigned to V (This does not mean that every potential will lead to a successful calculation of eigenvalues). For demonstration and comparison purposes, like here, one encounters the problem that the exact answers are no longer known. This makes it more difficult to assess the accuracy and other qualities of the code. As an example where some instructive comparisons are possible, we consider the two-dimensional anharmonic oscillator obeying the Schrödinger equation,

$\begin{array}{ll}\frac{1}{2}\left(-\frac{d^2}{d{x}^2}-\frac{d^2}{d{y}^2}+x^4+6x^2{y}^2+y^4\right){\Psi }_E(x,y)=E\Psi (x,y).& (18)\end{array}$

By construction, this problem has separable solutions of the form

$\begin{array}{ll}{\Psi }_E(x,y)={\psi }_m(\xi ){\psi }_n(\eta ),\text{with }\xi =(x+y)/\sqrt{2},\eta =(x-y)/\sqrt{2},& (19)\end{array}$

where the factors ψ obey a one-dimensional equation,

$\begin{array}{ll}\left(-\frac{d^2}{d{\xi }^2}+{\xi }^4\right){\psi }_m(\xi )={\epsilon }_m{\psi }_m(\xi ),& (20)\end{array}$

and E = ε_m + ε_n. As mentioned in the introduction, equations like the latter have been quite intensely studied in the literature. Accurate values for the even parity eigenvalues of Equation (20) can for instance be found in Table 1 of [9]. In Table 2, we list the 10 lowest eigenvalues to 30 decimals precision, calculated by the very-high-precision method described in [14]. Hence, for practical purposes all ε_m of interest can be considered exactly known. This means that the eigenvalues E of Equation (18) can also be considered exactly known. We list the 22 lowest ones of them in Table 3. These are the values we want to compare against the standard solution methods. The latter make no use of the separability property of the problem, which anyway is destroyed by the lattice approximation.

TABLE 2

Table 2. The 10 lowest eigenvalues of the quantum anharmonic oscillator, calculated to high precision by the method described in [14], from the Schrödinger equation $\left(-\frac{d^2}{d{\xi }^2}+{\xi }^4\right){\psi }_n(\xi )={\epsilon }_n{\psi }_n(\xi )$.

TABLE 3

Table 3. The 22 lowest eigenvalues E of the two-dimensional quantum anharmonic oscillator, as defined by the Schrödinger equation $\frac{1}{2}\left(-\frac{d^2}{d{x}^2}-\frac{d^2}{d{y}^2}+x^4+6x^2{y}^2+y^4\right){\Psi }_E(x,y)=E{\Psi }_E(x,y),$ displayed to 30 decimals accuracy.

The first column of Table 3 associates a symmetry classification (P_x, P_y) to each eigenvalue E, or rather to its corresponding eigenfunction Ψ_E(x, y). Since Equation (18) are invariant under reflections,

$\begin{array}{l}{{P}}_x:x\to -x\text{or }{{P}}_y:y\to -y,\end{array}$

all eigenfunctions can be constructed to transform symmetrically (S) or anti-symmetrically (A) under such reflections. For m < n, such a construction is

$\begin{array}{ll}{\Psi }_{mn}^{(\pm )}(x,y)=\frac{1}{\sqrt{2}}\left[{\psi }_m(\xi ){\psi }_n(\eta )\pm {\psi }_n(\xi ){\psi }_m(\eta )\right].& (\mathrm{21a})\end{array}$

For m = n there is only one possibility,

$\begin{array}{ll}{\Psi }_{mm}(x,y)={\psi }_m(\xi ){\psi }_m(\eta ).& (\mathrm{21b})\end{array}$

By use of the properties that

$\begin{array}{r}{\psi }_m(-\xi )={(-1)}^{m-1}\psi (\xi ),{{P}}_x:(\xi ,\eta )\to -(\eta ,\xi ),\text{and }\\ {{P}}_y:(\xi ,\eta )\to (\eta ,\xi ),\end{array}$

we find that

$\begin{array}{l}{\Psi }_{mn}^{(\pm )}(-x,y)=\pm {(-1)}^{m+n}{\Psi }_{mn}^{(\pm )}(x,y),\text{and}\\ {\Psi }_{mn}^{(\pm )}(x,-y)=\pm {\Psi }_{mn}^{(\pm )}(x,y).& (22)\end{array}$

and further that Ψ_mm(−x, y) = Ψ_mm(x, − y) = Ψ_mm(x, y). The conclusion is that in an exact calculation the states Ψ_mn will be double degenerate when m ≠ n, with parities (P_x, P_y) equal to (S, S) and (A, A) when m, n are both even or both odd, otherwise with parities (S, A) and (A, S). The states Ψ_mm are singlets with parities (S, S). The first column of Table 3 is constructed according to these rules.

Table 4 displays the results of some standard numerical solutions to Equation (18), “automagically” generated in the same way as the previous treatments of the harmonic (linear) oscillators. In the second column we show the results of using the minimal 5-point stensil approximation of the Laplace operator on a 1, 024 × 1, 024 lattice (approximating the whole space). The resulting numerical problem is solved with the eigsh sparse solver. The numerical accuracy is indicated by an underscore of the first inaccurate position, when taking proper roundoffs into account: The exact and numerical results are rounded off to the same number of digits, and compared; the underscore indicates the first position where the results differ.

TABLE 4

Table 4. Numerical calculations of the lowest eigenvalues of the two-dimensional quantum anharmonic oscillator, by various approximations and lattice sizes.

As can be seen, the results are less than impressive, taking into account the amount computational work invested. One straightforward improvement is to utilize the reflection symmetries of the problem to reduce the magnitude of the problem (with the same lattice cell size δx²) by a factor 4, or to reduce the lattice cell size δx² (with the same problem magnitude) by a factor 4. Another option is to use a higher order stensil approximation like (13). However, as already discussed in section 5, an even better option (for this class of problems) is to use a FFT type of approximation of the Laplace operator. The resulting eigenvalues are listed in columns 3–5, for various lattice sizes approximating the upper right quadrant (x ≥ 0, y ≥ 0) of space. For each lattice size the problem must be solved 4 times, with symmetric (S) and anti-symmetric (A) boundary conditions at the axes x = 0 and y = 0.

By symmetry under interchange, x ↔ y, we expect the (S, A) and (A, S) to give identical results (as long as the lattice approximation respects this symmetry). As can be seen, the numerical results satisfy the symmetry within a numerical accuracy of few × 10⁻¹², regardless how close the results are to the exact values. The degeneracy of states with (S, S), respectively, (A, A) symmetry cannot be deduced in the same way from the lattice approximated problem. In the infinite space formulation the problem is separable, which in turn implies this degeneracy. However, the lattice approximation introduces boundaries that are non-factorizable in the (ξ, η)-coordinates. This means that the problem is no longer separable in the lattice approximation. As a result the degeneracy of the (S, S) and (A, A) energies are split by a much larger amount, of the same order as the difference between exact and numerical results. (In this case, the lattice problem could be made separable by a rotation of the lattice orientation by 45 degrees.)

We observe that even a 2⁴ × 2⁴ lattice with in the “FFT approximated” Laplace operator provide almost equally accurate results as a 2¹⁰ × 2¹⁰ lattice with the 5-stensil approximation. The results from a 2⁵ × 2⁵ lattice seems more than sufficient for practical purposes (say compared to experimental obtainable accuracy), with little to be gained by further decrease of the lattice length δx.

The computation times for the “FFT approximation” are about 0.06, 0.8, and 75 s for respectively 16 × 16, 32 × 32, and 128 × 128 lattice sizes. For the same number of lattice points, the 5-stensil formulation may lead to somewhat shorter computation times. But this is completely offset by the need to work with a much larger number of lattice points: The computation time for the 1, 024 × 1, 024 stensil approximation was about 30 min.

The Python package described in this paper is available at [19].

Data Availability Statement

The raw data [19] supporting the conclusions of this article will be made available by the authors, without undue reservation.

Author Contributions

All authors listed have made a substantial, direct and intellectual contribution to the work, and approved it for publication.

Conflict of Interest

The authors declare that the research was conducted in the absence of any commercial or financial relationships that could be construed as a potential conflict of interest.

Acknowledgments

An early version of this work was presented at the IAENG International MultiConference of Engineers and Computer Scientists, Hong Kong March 18–20, 2015, as documented in [20].

AM would like to thank Mathematics Teaching and Learning Research Group within The department of mathematics, Bodø, Nord University for partial support.

Footnotes

1. ^The general formula is $g_n=\left(\frac{d-1+n}{d-1}\right)$.

2. ^The error of a back-and-forth FFT is a few times the numerical accuracy, i.e., in the range 10⁻¹⁴ to 10⁻¹⁵. with double precision numbers. However, when an error of this order is multiplied by k² it can be amplified by several orders of magnitude. Hence, the range of k²-values should not be chosen significantly larger than required to represent ψ(r) to sufficient accuracy.

References

1. Dirac PAM. The fundamental equations of quantum mechanics. Proc R Soc A. (1925) 109:642–53. doi: 10.1098/rspa.1925.0150

CrossRef Full Text | Google Scholar

2. Schrödinger E. An undulatory theory of the mechanics of atoms and molecules. Phys Rev. (1926) 28:1049–70. doi: 10.1103/PhysRev.28.1049

CrossRef Full Text | Google Scholar

3. Pauli W. Über die Wasserstoffspektrum vom Standpunkt der neuen Quantenmechanik. Zeitsch Phys. (1926) 36:336–63.

Google Scholar

4. Heisenberg W. Über quantentheoretische Umdeutung kinematischer und mechanischer Beziehungen [Quantum-theoretical reinterpretation of kinematic and mechanical relations]. Zeitsch Phys. (1925) 33:879–93.

Google Scholar

5. Bender CM, Wu TT. Analytic structure of energy levels in a field-theory model. Phys Rev Lett. (1968) 21:406–9. doi: 10.1103/PhysRevLett.21.406

CrossRef Full Text | Google Scholar

6. Bender CM, Wu TT. Anharmonic oscillator. Phys Rev. (1969) 184:1231–60. doi: 10.1103/PhysRev.184.1231

CrossRef Full Text | Google Scholar

7. Bender CM, Wu TT. Large-order behavior of perturbation theory. Phys Rev Lett. (1971) 27:461–5. doi: 10.1103/PhysRevLett.27.461

CrossRef Full Text | Google Scholar

8. Bender CM, Wu TT. Anharmonic oscillator. II. A study of perturbation theory in large order. Phys Rev D. (1973) 7:1620–36. doi: 10.1103/PhysRevD.7.1620

CrossRef Full Text | Google Scholar

9. Bender CM, Olaussen K, Wang PS. Numerological analysis of the WKB approximation in large order. Phys Rev D. (1977) 16:1740–8. doi: 10.1103/PhysRevD.16.1740

CrossRef Full Text | Google Scholar

10. Zinn-Justin J, Jentschura UD. Multi-instantons and exact results I: conjectures, WKB expansions, and instanton interactions. Ann Phys. (2004) 313:197–267. doi: 10.1016/j.aop.2004.04.004

CrossRef Full Text | Google Scholar

11. Zinn-Justin J, Jentschura UD. Multi-instantons and exact results II: specific cases, higher-order effects, and numerical calculations. Ann Phys. (2004) 313:269–325. doi: 10.1016/j.aop.2004.04.003

CrossRef Full Text | Google Scholar

12. Noreen A, Olaussen K. Quantum loop expansion to high orders, extended borel summation, and comparison with exact results. Phys Rev Lett. (2013) 111:040402. doi: 10.1103/PhysRevLett.111.040402

PubMed Abstract | CrossRef Full Text | Google Scholar

13. Janke W, Kleinert H. Convergent strong-coupling expansions from divergent weak-coupling perturbation theory. Phys Rev Lett. (1995) 75:2787–91. doi: 10.1103/PhysRevLett.75.2787

PubMed Abstract | CrossRef Full Text | Google Scholar

14. Mushtaq A, Noreen A, Olaussen K, Øverbø I. Very-high-precision solutions of a class of Schrödinger type equations. Comput Phys Commun. (2011) 182:1810–3. doi: 10.1016/j.cpc.2010.12.046

CrossRef Full Text | Google Scholar

15. Noreen A, Olaussen K. High precision series solutions of differential equations: ordinary and regular singular points of second order ODE's. Comput Phys Commun. (2012) 183:2291–7. doi: 10.1016/j.cpc.2012.05.015

CrossRef Full Text | Google Scholar

16. Van Der Walt S, Colbert SC, Varoquaux G. The NumPy array: a structure for efficient numerical computation. Comput Sci Eng. (2011) 13:22–30. doi: 10.1109/MCSE.2011.37

CrossRef Full Text | Google Scholar

17. Jones E, Oliphant T, Peterson P. SciPy: Open Source Scientific Tools for Python. (2014). Available online at: http://www.scipy.org/

18. Oliphant TE. Python for scientific computing. Comput Sci Eng. (2007) 9:10–20. doi: 10.1109/MCSE.2007.58

CrossRef Full Text | Google Scholar

19. Mushtaq A, Noreen A, Olaussen A. Python Package for Numerical Solutions of Quantum Mechanical Eigenvalue Problems (2020). doi: 10.6084/m9.figshare.127655

CrossRef Full Text

20. Noreen A, Olaussen A. A python class for higher-dimensional schrödinger equations. In: Proceedings of the International MultiConference of Engineers and Computer Scientists IMECS 2015. Hong Kong (2015). p. 206–11.

Google Scholar