The hardwall method of solving the radial Schr\"odinger equation and unmasking hidden symmetries

TL;DR

The paper introduces a direct method for solving the radial Schrödinger equation that efficiently finds all eigenvalues without iterative procedures and reveals hidden symmetries through a 'hardwall' approach.

Contribution

A novel non-iterative technique for solving the radial Schrödinger equation that uncovers hidden symmetries and simplifies eigenvalue determination.

Findings

Automatically determines entire spectrum of eigenvalues

Unmasks accidental degeneracies due to hidden symmetries

Effective for various potentials like Coulomb and Woods-Saxon

Abstract

Solving for the bound state eigenvalues of the Schr\"odinger equation is a tedious iterative process when the conventional shooting or matching method is used. In this work, we bypass the eigenvalue's dependence on the eigenfunction by simply trying out all eigenvalues to a desired accuracy. When the eigenvalue is known, the integration for the eigenfunction is then trivial. At a given energy, by outputting the radial distance at which the wave function crosses zero (the hardwall radius), this method automatically determines the entire spectrum of eigenvalues of the radial Schr\"odinger equation without iterative adjustments. Moreover, such a spherically symmetric "hardwall" can unmask "accidental degeneracy" of eigenvalues due to hidden symmetries. We illustrate the method on the Coulomb, harmonic, Coulomb+harmonic, and the Woods-Saxon potentials.

Tables

Table 1. Table 1: Calculated eigenvalues of the Coulomb potential. The first column is the exact results − 1 / ( 2 ​ n 2 ) 1 2 superscript 𝑛 2 -1/(2n^{2}) with n = l + n r + 1 𝑛 𝑙 subscript 𝑛 𝑟 1 n=l+n_{r}+1 .

$-1/\left(2n^2\right)$	$l=0$	$l=1$	$l=2$	$l=3$	$l=4$
-0.50000	-0.49998
-0.12500	-0.12499	-0.12499
-0.05555	-0.05555	-0.05555	-0.05554
-0.03125	-0.03124	-0.03124	-0.03124	-0.03123
-0.02000	-0.01984	-0.01987	-0.01991	-0.01996	-0.01998

Table 2. Table 2: Calculated eigenvalues of the harmonic potential. The first column is the exact results ( l + 2 ​ n r + 3 / 2 ) 𝑙 2 subscript 𝑛 𝑟 3 2 (l+2n_{r}+3/2) .

$\left(l+2n_r+3/2\right)$	$l=0$	$l=1$	$l=2$	$l=3$	$l=4$	$l=5$	$l=6$
1.50000	1.50002
2.50000		2.50001
3.50000	3.50001		3.50002
4.50000		4.50001		4.50002
5.50000	5.50001		5.50001		5.50001
6.50000		6.50001		6.50001		6.50002
7.50000	7.50001		7.50001		7.50001		7.50002

Table 3. Table 3: Calculated eigenvalues of the Coulomb+harmonic potential. There is now no degeneracy;

$n_r$	$l=0$	$l=1$	$l=2$	$l=3$	$l=4$	$l=5$	$l=6$
0	0.17968
1		1.70903
2	2.50001		2.88224
3		3.80193		3.97553
4	4.63196		4.93068		5.03608
5		5.86036		6.00654		6.07947
6	6.71260		6.96584		7.05815		7.11255

Table 4. Table 4: Calculated eigenvalues as compared to published values GVB WS for the Woods-Saxon potential for angular momentum channels l = 0 , 1 , 2 𝑙 0 1 2 l=0,1,2 . In each channel, only energies of the even radial quantum state n r subscript 𝑛 𝑟 n_{r} are compared. The unit of energy is MeV.

		Exact results			Calculated results
$n_r$	$l=0$	$l=1$	$l=2$	$l=0$	$l=1$	$l=2$
0	-49.457788728	-48.951731623	-48.34981052	-49.4570	-48.9510	-48.3485
2	-46.290753954	-45.237176986	-44.121537377	-46.2905	-45.2370	-44.1215
4	-41.232607772	-39.767208069	-38.253426539	-41.2325	-39.7670	-38.2530
6	-34.672313205	-32.868392986	-31.026820921	-34.6720	-32.8680	-31.0265
8	-26.873448916	-24.794185466	-22.689041510	-26.8730	-24.7940	-22.6890
10	-18.094688282	-15.812724871	-13.522303352	-18.0945	-15.8125	-13.5220
12	-8.676081670	-6.308097192	-3.972491432	-8.6760	-6.3080	-3.9720

Equations

$$-\frac{{\hslash }^2}{2m}\left(\frac{{\partial }^2}{\partial r^2}-\frac{l(l+1)}{r^2}\right)u(r)+V(r)u(r)=Eu(r),$$

$$-\frac{{\hslash }^2}{2m{a}^2}\left(\frac{{\partial }^2}{\partial r^{\ast 2}}-\frac{l(l+1)}{r^{\ast 2}}\right)u(r^{\ast })+V(r^{\ast }a)u(r^{\ast })=Eu(r^{\ast }).$$

$$-\frac{1}{2}\left(\frac{{\partial }^2}{\partial r^{\ast 2}}-\frac{l(l+1)}{r^{\ast 2}}\right)u(r^{\ast })+v(r^{\ast })u(r^{\ast })=E^{\ast }u(r^{\ast }),$$

$$v(r^{\ast })=\frac{V(r^{\ast }a)}{e_0}=\frac{m{a}^{2}V(r^{\ast }a)}{{\hslash }^2}.$$

$$v(r^{\ast })=-\frac{m{a}^{2}k{e}^2}{r^{\ast }a{\hslash }^2}=-\frac{1}{r^{\ast }},$$

$$v(r^{\ast })=\frac{m{a}^{2}m{\omega }^2{a}^2{r}^{\ast 2}}{2{\hslash }^2}=\frac{1}{2}{r}^{\ast 2},$$

$$\frac{{\partial }^{2}u(r)}{\partial r^2}=f(r)u(r),$$

$$f(r)=\frac{l(l+1)}{r^2}+2V(r)-2E.$$

$$(1-\frac{1}{12}{\Delta r}^2{f}_{n+1})u_{n+1}=(2+\frac{5}{6}{\Delta r}^2{f}_n)u_n-(1-\frac{1}{12}{\Delta r}^2{f}_{n-1})u_{n-1},$$

$$u_{n+1}=\frac{(12-10g_n)u_n-g_{n-1}{u}_{n-1}}{g_{n+1}}.$$

$$V_{\mathrm{WS}}(r)=\frac{u_0}{1+t}+\frac{u_{1}t}{(1+t{)}^2},t=\exp [\frac{r-r_0}{a}],$$

Full text

I Introduction

Solving the radial Schrödinger equation for a central potential $V(r)$,

[TABLE]

where $u(r)=rR_{nl}(r)$, has been considered an important teaching tool in this Journalbol and a standard subject in Koonin’s early computational physics text.koo However, in recent years, this topic is not discussed in most computational physics texts,gou ; gio ; dev ; ves ; lan or only alluded to briefly.gio Perhaps, this is related to the difficulty of implementing the conventional shooting or matching method bol ; gou ; gio ; dev ; ves ; lan of solving the eigenvalue problem. While more sophisticated methodsixa ; GVB WS ; sac ; aug have been suggested in the literature, we introduce here a classroom-tested, start-from-scratch method, that must be the simplest among all known ways of solving the bound state problem of the radial Schrödinger equation.

Our simplification follows from two key ideas: 1) The eigenvalue $E$ is usually determined by requiring the eigenfunction to satisfy certain boundary conditions. The adjustment of the eigenvalue by shooting or matching the eigenfunctiongou ; gio ; dev ; ves ; lan to satisfy the boundary condition is a tedious iterative process. In this work, we completely by-pass this adjustment process by sweeping through all values of $E$ to a desire accuracy. That is, we simply try all values of $E$, up to a certain precision, by brute-force. When $E$ is known, the eigenfunction can be obtained easily. This then breaks the dependence of the eigenvalue to the eigenfunction. 2) To decide which $E$ is the correct eigenvalue, we observe that every $E$ value is an eigenvalue of the given potential plus an infinite potential barrier at some radius $C$. The correct eigenvalue $E$ is then obtained in the limit of $C\rightarrow\infty$. Given $E$, $C$ is the simply the radial distance wherever the wave function vanishes. This is then the “hardwall” condition to be described in Section III. Our method takes advantage of the fact that it is far easier to find $C$ as a function of $E$, than the other way around, as in the matrix method.aug Also in contrast to the matrix method, no basis functions need to be assumed, no matrix elements need to computed and no matrix needs to be diagonalized.

This brute-force approach is conceptually very simple and can be easily understood by any undergraduate taking a first course in quantum mechanics. However, this approach, by leveraging the power of modern computers, can decipher the entire spectrum of a central potential automatically without piecemeal adjustments. Moreover, this method can unmask hidden symmetries by revealing “accidental degeneracy” of eigenvalues and their patterns of degeneracy.

II Solving the radial Schrödinger equation numerically

To numerically solve the reduced radial Schrödinger equation (1), we will first cast it into a dimensionless form in terms of units most appropriate for the equation. This digression seems appropriate since most current textsgou ; gio ; dev ; ves ; lan do not emphasize this important point.

To arrive at a dimensionless equation, one first sets $r=r^{*}a$, where $a$ is a unit of length to be determined, and $r^{*}$ is just a dimensionless number. The equation then reads

[TABLE]

Since ${\rm e}_{0}\equiv\hbar^{2}/ma^{2}$ is a unit of energy, it is natural to write $E=E^{*}e_{0}$, so that (1) has the dimensionless form

[TABLE]

where the dimensionless potential is given by

[TABLE]

One then chooses $a$ so that $v(r^{*})$ is as simple as possible. For example, in the case of the Coulomb potential with $V(r)=-ke^{2}/r$, the dimensionless potential is

[TABLE]

when one chooses $a=\hbar^{2}/(mke^{2})$, which is the Bohr radius. The unit of energy ${\rm e}_{0}$ is then the Hartree (Ha) and the spectrum of the hydrogen atom is $E=-1/(2n^{2})$ Ha. (Note that 1 Ha=2 Ry (Rydberg)=2(13.6 eV). ) This choice defines the atomic units, most appropriate for dealing with atomic problems. Another example is the the harmonic potential, with $V(r)=(1/2)m\omega^{2}r^{2}$. In this case

[TABLE]

with the natural choice of $a=\sqrt{\hbar/m\omega}$, which is the harmonic length, and ${\rm e}_{0}=\hbar^{2}/ma^{2}=\hbar\omega$. These are then the harmonic units.

The advantage of solving a dimensionless equation is that one solves an entire class of equations at once, not just one particular equation. For example, if one is interested in solving the Coulomb potential with a nuclear charge $Ze$, then there is no need for a new calculation. One simply changes the units from $k\rightarrow kZ$. One then sees immediately that the Bohr radius will shrink by a factor of $Z$ and the energy grow by a factor of $Z^{2}$. If one is interested in the case of muonic atoms, where the electron is replaced by the muon 207 times as massive, then the radius of muonic-hydrogen will be 207 times smaller with binding energy $(207)^{2}$ times greater. Similarly for the harmonic case; the dimensionless equation solves an entire class of problems with any values of $m$ and $\omega$. These important insights are lost when explicit units are used, as done in some computational texts.dev

From this point forward, we will drop all asterisks ∗ in referencing the dimensionless radial equation (2).

The eigenvalue problem associated with the radial equation (2) is that one must determine $E$ subject to the boundary conditions $u(0)=0$ and $u(\infty)=0$. The shooting method starts out at $u(0)=0$, integrates out to some large value of $r=R$, and adjusts $E$ so that $u(R)=0$. (The needed $R$ value is different for different eigenstates.)

In this work, as in the shooting method, we will also start at the origin with $u(0)=0$ and integrate outward, but will not need to impose the condition $u(R)=0$ for an a priori unknown $R$ value. Instead, we will introduce the hardwall condition to be described in the next Section.

The dimensionless radial equation (2) can be further arranged as

[TABLE]

where

[TABLE]

It is well-known that (3) can be efficiently solved by use of the fourth-order Numerovnum algorithm, which has been derived numerous time in this Journal,cho ; bol and elsewhere:Hartree ; koo

[TABLE]

where $u_{n}=u(n{\Delta r})$ and $f_{n}=f(n{\Delta r})$. For a given ${\Delta r}$, one can define $g_{n}=1-\frac{1}{12}{\Delta r}^{2}f_{n}$ so that the algorithm reads simply,

[TABLE]

The above Numerov algorithm is used to obtain all results reported in this work.

III The hardwall condition

The hardwall method replaces the $u(\infty)=0$ boundary condition with the following hardwall condition: for a set of given $l$ and $E$ values, $u(r)$ is iterated outward using (4) at a chosen ${\Delta r}$, with $u_{0}=0$ and $u_{1}$ an arbitrary but small real number, out to a large value of $r=R$, depending on the number of eigenvalues desired. Whenever $u_{n}$ crosses zero, i.e, $u_{n-1}u_{n}<0$, we output $(C,E)$ where $E$ is the given energy and $C=(n-1/2){\Delta r}$ is the hardwall radius. That is, whenever $u_{n}$ crosses zero, $E$ is the exact energy corresponding to the original potential $v(r)$ plus an infinite potential (the “hardwall”) at $r=C$. In other words, every value of $E$ is an eigenvalue for $v(r)$ plus a hardwall potential at some values of $C$. For each value of $l$, one does the above outward integration of $u_{n}$ over a set of prescribed $E$ values spanning the bound state range of $v(r)$. When $E$ versus $C$ is plotted after one has swept through an interval of $E$, one will see that $E$ converges to the eigenvalue of $v(r)$ in the limit of large $C$.

Let’s see how the hardwall method works for the null-potential case of $v(r)=0$. In this case, the solutions to the radial Schrödinger equation are spherical Bessel functions $j_{l}(k_{ln}r)$, with energy eigenvalues $E_{ln}=(1/2)k_{ln}^{2}.$ For $l=0$, $j_{0}=\sin(k_{0n}r)/k_{0n}r$. For an infinite wall at the hardwall radius $C$, the wave function must vanish, forcing $k_{0n}C=n\pi$ and $E=(1/2)(n\pi/C)^{2}.$ In Fig.1, we compare the outputted $(C,E)$ of the hardwall method with this analytical result. We shall refer to this $(C,E)$ plot of the hardwall method as a “$C$-scan” of the potential. For this calculation, we set $l=0$ and iterate over an outer loop of 50 values of $E_{i}=i\Delta E$, with $\Delta E=0.1$ and an inner loop of 1100 values of $r=j{\Delta r}$ with ${\Delta r}=0.01$ out to $C=11$ (so that only 10 eigenvalues are visible). This calculation demonstrates the essential characteristic of the method: 1) At any fixed value of $C$, the $C$-scan gives the correct eigenvalues of the hardwall potential at $r=C$ when $\Delta E$ is sufficiently fine. The power of the method resides in the fact at a fixed $E$, multiple values of $C$ are found simultaneously. If one does the opposite, fixing $C$ then finding $E$, then $E$ must be determined one by one. 2) As $C\rightarrow\infty$, all eigenvalues approach zero, as they should. However, higher eigenvalues approach zero more slowly at larger $C$. Thus when $v(r)\neq 0$, the $C$-scan will approach the eigenvalues of $v(r)$, also with higher eigenvalues converging at larger values of $C$. 3) A single sweep in energy produces all the eigenvalues at the same time. An accuracy of 5-6 significant digits in determining $E$ is easily achievable on a laptop computer.

We will now see how effective is this method when applied to physical potentials.

IV Hidden symmetries in the Coulomb and harmonic potentials

For the Coulomb potential $v_{\rm C}(r)=-1/r,$ the $C$-scan showing the convergence of the eigenvalues is shown in Fig.2. One immediately notices something unusual. In the large $C$ limit, distinct higher energy levels corresponding to different radial quantum numbers $n_{r}$ and $l$ coalesce into a single level. Since the hardwall is spherically symmetric, its addition will not alter the spherical symmetry of the Coulomb potential. For spherical symmetry, each energy level with quantum number $l$, regardless of $n_{r}$, is $2l+1$ degenerate, meaning that all $2l+1$ states with azimuthal quantum number $m=-l,\cdots,-1,0,1,\cdots,l$ have the same energy. However, as the hardwall is gradually removed in letting $C\rightarrow\infty$, further degeneracy is revealed, due to the hidden $SO(4)$ symmetrybay of the Coulomb potential.

The pattern of this degeneracy is also visible. The ground state $(l,n_{r})=(0,0)$ is non-degenerate. The next two states $(0,1)$, $(1,0)$ are degenerate. The next three states $(0,2)$, $(1,1)$, $(2,0)$ are degenerate, etc., with the energy only depending on the sum of $l$ and $n_{r}$ as $E=-1/(2n^{2})$, with $n=l+n_{r}+1$. Thus the hardwall method not only computes the eigenvalues, its $C$-scan also reveal the hidden symmetry and the pattern of degeneracy of the Coulomb potential. The resulting eigenvalues are shown in Table 1.

For the harmonic potential $v_{\rm H}(r)=(1/2)r^{2},$ the corresponding $C$-scan is shown in Fig.3. Again, one observes energy degeneracy at large value of $C$. One can therefore conclude that the 3D harmonic oscillator also has a hidden symmetry (that of $SU3$).fra Let’s see whether one can also deduce its pattern of degeneracy. There are two states $(l,n_{r})=(0,0),(1,0)$ that are non-degenerate. There are two states that are twice degenerate: [(0,1), (2,0)] and [(1,1), (3,0)], thrice degenerate: [(0,2), (2,1), (4,0)] and [(1,2), (3,1), (5,0)], and four-fold degenerate, etc. The degeneracy pattern is therefore for $l$ even, [(0,0)], [(0,1), (2,0)], [(0,2), (2,1), (4,0)], [(0,3), (2,2), (4,1), (6,0)] etc.. And for $l$ odd, [(1,0)], [(1,1), (3,0)], [(1,2), (3,1), (5,0)], [(1,3), (3,2), (5,1), (7,0)] etc.. (Degenerate levels are grouped together by square brackets.) Both patterns of degeneracy can be accounted for if the energy only depends on $l+2n_{r}$. Therefore, one can deduce that the energy spectrum must be given by $l+2n_{r}+3/2$ where $3/2$ is the ground state energy of (0,0). This analytical result is compared to the computed spectrum given in Table 2.

The hidden symmetry of the Coulomb and the harmonic potential is unique in each case. If one simply adds the Coulomb and the harmonic potential together, then the hidden symmetry is destroyed for the combined potential. This is illustrated in the $C$-scan of Fig.4. One sees that energy levels which were $k$-fold degenerate in the harmonic oscillator case, now split into $k$ energy levels in the large $C$ limit. The computed eigenvalues are shown in Table 3. Since the Coulomb potential is only strong near the origin, it acts as a perturbation, and only modifies the low-lying spectrum of the harmonic oscillator.

V Woods-Saxon potential

The Woods-Saxon potential WS orig ; GVB WS describes the effective average force confining each nucleon to the interior of the nucleus in the form of

[TABLE]

with $u_{0}$ and $u_{1}$ fixing the potential well depth, $a$ the surface thickness of the nucleus, and $r_{0}$ the nuclear radius proportional to the mass number. We use the parameterization by Vanden Berghe, Fack, and De MeyerGVB WS with $u_{0}=-50$ MeV, $r_{0}=7$ fm, $a=0.6$ fm, and $u_{1}=-u_{0}/a$. In this calculation, we use the dimensionless form of the equation but restate results in units as given above for comparison. The convergence of the eigenvalues, with ${\Delta r}=0.001$, $\Delta E=0.0005$, is shown in Fig.5. Since the Woods-Saxon potential has no special symmetry, its $C$-scan shows close-by energy levels, but no actual degeneracy. The converged eigenvalues, a selected set of 21, are compared with exact results in Table 4. Our eigenvalues agree with published results up to five significant digits.

VI Conclusion

In this work, we have demonstrated an extremely simple method of solving the radial Schrödinger equation by trying out all eigenvalues systematically, up to a certain precision. The hardwall method relies on the insight that any trial eigenvalue $E$ is an exact eigenvalue for the given potential plus a hardwall at some radius $C$. By plotting $E$ verses $C$, the eigenvalues can be determined to five or more significant digits in the large $C$ limit. Moreover, this $C$-scan of the spectrum can detect hidden symmetries by revealing additional degenerate energy levels beyond spherical symmetry and their patterns of degeneracy. The method is therefore a valuable tool, not just for solving the radial Schrödinger equation, but also for teaching the emergence of “accidental degeneracies” from hidden symmetries.

Finally, the method can also be applied to a strictly 1D Schrödinger equation with a symmetric potential $v(x)=v(-x)$. By setting $l=0$, the same starting condition $u_{0}=0$ now gives the odd-state solutions. The even-state solutions can be obtained by starting with $u_{0}=1$ and $u_{1}=1$.