Solution of the 1d Schrödinger Equation for a Symmetric Well

Abstract

We suggest a mathematical potential well with spherical symmetry and apply to the 1d Schrödinger equation. We use some well known techniques as Stationary Perturbation Theory and WKB to gain insight about the solutions and compare them each other. Finally, we solve the 1d Schrödinger equation using a numerical approach with the so-called Numerov technique for comparison. This can be a good exercise for undergrad students to grasp the above cited techniques in a quantum mechanics course.

Keywords:
Potential Well; WKB; Numerov

1. Introduction

In Quantum Mechanics books we usually find trivial examples when Stationary Perturbation Theory (SPT), Wentzel-Kramers-Brillouin (WKB) and even other techniques are discussed. The interesting applications are left to some complicated exercises at the end of the chapter [¹- ⁶]. Most of the books apply those techniques to the simple harmonic oscillator with $V(x)~x^2$ or at most to the $x^4$ potential in 1d. Increasingly, the computer is becoming part of the physics courses and it would be very interesting to have certain classes of problems to be solved in a Quantum Mechanical course.

Nowadays, several numerical techniques has been successfully employed to solve the Schrödinger equation to obtain both the energy levels and the respective wave functions. In particular, for this kind of differential equation a powerful method among others is that proposed by Boris Vasil’evich Numerov ^[7]. This method takes the advantage of the fact that the Schrödinger equation is an eigenvalue equation to handle the wave function subjected to certain boundary conditions in order to minimize the energy for each level. Its algorithm is very efficient and converges very fast with at least $O(h^6)$ of precision. Under this perspective, it is a good exercise to solve the 1d-Schrödinger equation using this technique to gain good insights about the Schrödinger equation during the classes ^[8].

In view of that, we suggest a mathematical potential well, which can be expanded in even exponent power series inside the well, becoming a good exercise to treat it perturbatively or using the WKB method or other well-know Quantum Mechanics techniques to compare with the powerful numerical results. In what follows we compare some states order by order to a good numerical calculation in determining the solutions of the 1d-Schrödinger equation.

This article is sketched as follows: in section 2 we present the mathematical potential well. In section 3 we apply the potential well to the 1d-Schrödinger equation and analyses some particular aspects for both analytical and WKB approximation. In section 4 we present the Stationary Perturbation Theory to calculate the energy levels and some wave functions to compare with. In section 5 we discuss about the Numerov's numerical approach to solve the 1d-Schrödinger equation. Finally, in section 7 we draw conclusions and perspectives.

2. The potential well

Now we will introduce a mathematical potential whose series expansion is interesting because its symmetry and divergence at specific points. In 3D this mathematical potential is given by

(1) $$V(r)=V_0\left(\frac{1-kr\cot (kr)}{{(kr)}^2}\right)$$

where k has m${}^{-1}$ units, $r=\sqrt{x^2+y^2+z^2}$ in m, and $V_0$ in J. k must be such that when r equals some value, say a, the product $kr$ must be equal to $\pi$. As a consequence $k=\pi /a$ and $V(r)$ diverges at $r=a$. Outside this interval the cotangent function makes the potential oscillates between regions with negative and positive divergences with $V_0$ positively defined, which difficult the analysis and we will not consider this situation in the present work.

In figure 1 we display the 3D plot for z-coordinate equals zero. If one varies z, the diameter of the circle at the bottom of the figure becomes narrower as $\left|z\right|$ grows up to a maximum value $z=a$, making $x=0$ and $y=0$ such that $r=a$, which is a point of divergence at the center of the circle. For pedagogical applications the z-coordinate will be set zero from now on. Note that the surface figure generated is formed basically by symmetrical curves with minimum values.

Figure 1
Surface plot of the mathematical potential in Eq. (1).

In particular, we are interested in those curves with minimum values given by

(2) $$\underset{r\to 0}{\lim }V(r)=\frac{V_0}{3}$$

i.e., those curves that passes through the point $(0,0,0)$. Let us take the curve belonging to the plane $(x,0,0)$, which has the limit (2) above and, consequently, $r=\left|x\right|$. However, as we already mentioned, Eq. (1) is spherically symmetric, which means that if we substitute $r=\left|x\right|$ by $-r=-\left|x\right|$, the potential remains the same. As a consequence we can drop the absolute value of $\left|x\right|$ and write only x to have a 1d version of this mathematical potential

(3) $$V(x)=V_0\left(\frac{1-kx\cot (kx)}{{(kx)}^2}\right)$$

This 1d potential is plotted in figure 2. As one can see, in the symmetrical interval $kx\in [-\pi ,\pi ]$, it takes a minimum value at $V_0/3$ and diverges at $\pm \pi$. This is the potential we will use to solve the 1d Schrödinger equation.

Figure 2
1d plot of the potential.

Figure 3
Assume that the interval is discretized $x\in [-\pi ,\pi ]$ in N equal parts $h=x_{k+1}-x_k$ or equivalently $h=2\pi /N$.

To assess the singularities of this potential we firstly consider the series expansion of the cotangent function around the origin given by

(4) $$\begin{array}{ll}\cot kx\hfill & =\frac{B_0}{kx}+{\displaystyle \sum _{n=1}^{\infty }\frac{{(-1)}^n{2}^{2n}{B}_{2n}{(kx)}^{2n-1}}{(2n)!}}{B}_{2n},\hfill \\ \hfill & 0<\left|kx\right|<\pi \hfill \end{array}$$

where $B_{2n}$ are the Bernoulli numbers: $B_0=1$, $B_1=-\frac{1}{2}$, $B_2=\frac{1}{6}$, $B_4=-\frac{1}{30}$, $B_6=\frac{1}{42}$, $B_8=\frac{-1}{30}$, etc. As expected, this function exhibits a singularity at the origin. If one multiplies both sides by $kx$, this series can be rewritten as

(5) $$kx\cot (kx)=B_0+{\displaystyle \sum _{n=1}^{\infty }\frac{{(-1)}^n{(2kx)}^{2n}}{(2n)!}}{B}_{2n}$$

where the singularity is now removed. If we use the value of $B_0=1$ and divide both sides by ${(kx)}^2$ with $x\ne 0$, we can rewrite this series as

(6) $$\begin{array}{ll}\left(\frac{1-kx\cot (kx)}{{(kx)}^2}\right)=-4{\displaystyle \sum _{n=1}^{\infty }}\hfill & \frac{{(-1)}^n{(2kx)}^{2n-2}}{(2n)!}{B}_{2n},\hfill \\ \hfill & \left|kx\right|<\pi \hfill \end{array}$$

This series converges for all the values in the interval $kx\in ]-\pi ,\pi [$, as one can check taking the limit $x\to 0$ in Eq. (6). Thus, the 1d potential can be expressed as

(7) $$V(x)=\{\begin{array}{ll}-4V_0{\displaystyle \sum _{n=1}^{\infty }\frac{{(-1)}^n{(2kx)}^{2n-2}}{(2n)!}}{B}_{2n}\hfill & \text{if}\left|kx\right|<\pi \hfill \\ \infty \hfill & \text{if}\left|kx\right|=\pi \hfill \end{array}$$

The first terms of the infinite series are given by

(8) $$\begin{array}{l}{\displaystyle \sum _{n=1}^{\infty }\frac{{(-1)}^n{(2kx)}^{2n-2}}{(2n)!}}{B}_{2n}=-\frac{B_2}{2}+\frac{B_4}{6}{(kx)}^2-\hfill \\ \frac{B_6}{45}{(kx)}^4+\frac{B_8}{630}{(kx)}^6+\cdots \hfill \end{array}$$

Notice that the potential will be given in even exponents of the series, which favors the use of the harmonic oscillator Hamiltonian as the unperturbed Hamiltonian when perturbation theory is employed. The potential plot in the $O(x^2)$ approximation is also displayed in figure 2. We will see the employment of approximated methods in the next sections.

3. 1d Schödinger equation

The 1d Schrödinger equation with is given by

(9) $$-\frac{{\hslash }^2}{2m}\frac{d^2\psi (x)}{d{x}^2}+V(x)\psi (x)=E\psi (x)$$

For this particular potential the 1d Schrödinger equation can be solved analytically in a very special situation. To see that, let us consider the potential given by Eq. (3) and define a new variable $\xi =kx$ such that Eq. (9) can be rewritten as

(10) $$\frac{d^2\varphi (\xi )}{d{\xi }^2}-2{\upsilon }_0\left(\frac{1-\xi \cot (\xi )}{{\xi }^2}\right)\varphi (\xi )=-2\epsilon \varphi (\xi )$$

where $\varphi (\xi )\equiv \psi (\xi /k)$, ${\upsilon }_0=\frac{m{V}_0}{{\hslash }^2{k}^2}$, and $\epsilon =\frac{mE}{{\hslash }^2{k}^2}$. We want to solve this equation for $\xi$ with the boundary conditions $\varphi (\pm \pi )=0$, where the potential diverges. Note that ${\upsilon }_0$ is dimensionless. $V_0$ can assume any positive value in dimensions of ${\hslash }^2{k}^2/m$ provided that E is greater than $V_0/3$. For simplicity we consider this positive value equals 1 from now on.

It is straightforward and important to notice that,

(11) $$\frac{1-\xi \cot (\xi )}{{\xi }^2}=-\frac{1}{\sin (\xi )}\frac{d}{d\xi }\left(\frac{\sin (\xi )}{\xi }\right),$$

then we plug the above identity in Eq. (10) to obtain

(12) $$\frac{d^2\varphi (\xi )}{d{\xi }^2}+2{\upsilon }_0\frac{\varphi (\xi )}{\sin (\xi )}\frac{d}{d\xi }\left(\frac{\sin (\xi )}{\xi }\right)+2\epsilon \varphi (\xi )=0.$$

As we can note, if we choose a Spherical Bessel Function $j_0(\xi )=\sin (\xi )/\xi$ we obtain,

(13) $$\frac{d^2}{d{\xi }^2}{j}_0(\xi )+\frac{2{\upsilon }_0}{\xi }\frac{d}{d\xi }{j}_0(\xi )+2\epsilon j_0(\xi )=0.$$

In this sense we must impose ${\upsilon }_0=1$ and $2\epsilon =1$ to obtain a common solution for Spherical Bessel Function of the first kind of the zeroth-order. Consequently we have

(14) $$V_0=\frac{{\hslash }^2{k}^2}{m}\text{and}E=\frac{{\hslash }^2{k}^2}{2m},$$

valid for a wave function $\varphi (\xi )=A\sin (\xi )/\xi$. We shall demonstrate a posteriori that this is actually the ground state wave function with eigenvalue $E=0.5{\hslash }^2{k}^2/m$ using a numerical approach to solve the 1d Schrödinger equation (9).

Now it is a good exercise to employ some approximate approach to have an idea about the energy levels for the bound states without too much effort. The Wentzel-Kramers-Brillouin (WKB) approximation is perfect for such a situation, because the bounded particle in this case has a sufficiently high momentum around the center of the potential as well as we can expect that its wave function varies rapidly with position, much more rapidly than the potential ^[1], as is the case when one considers the particle in a box.

To this end we can approximate the potential to $O(x^2)$ with the aid of Eq.'s (7) and (8), as plotted in figure 2. The boundary condition for the potential well in Eq. (9) with WKB method is such that $V(\pm \pi )\to \infty$ implies $\varphi (\pm \pi )=0$. In the pedagogical exercise of one-dimensional square well potential with perfectly rigid walls in a symmetrical interval $-a\le x\le a$ where $V(x)=0$ and $V(\pm a)=\infty$, the WKB method gives us an exact solution of energy levels [³- ⁵]. Clearly the potential (7) has rigid walls in $-\pi \le kx\le \pi$ in $O(x^2)$(as can be seen in figure 2), so that we can use WKB method to obtain approximated energy levels. Thus, in the classical region of the potential $E>V(x)$, the WKB method gives the following solution

(15) $$\psi (x)\propto \frac{C}{\sqrt{p(x)}}{e}^{\pm i\sigma (x)}$$

where $\sigma (x)=\frac{1}{\hslash }{\displaystyle \int p}(x)dx$, and $p(x)$ is the momentum given by

(16) $$p(x)=\sqrt{2m[E-V(x)]}.$$

The linear combination of solutions (15) is also a solution and

$$\psi (x)\approx \frac{1}{\sqrt{p(x)}}\left[C_1{e}^{i\sigma (x)}+C_2{e}^{-i\sigma (x)}\right].$$

Due to the symmetry of the well $-\pi \le kx\le \pi$, it admits even wave functions

(17) $${\psi }_{even}(x)\approx \frac{A}{\sqrt{p(x)}}\cos [\sigma (kx)],$$

and odd wave functions

(18) $${\varphi }_{odd}(x)\approx \frac{B}{\sqrt{p(x)}}\sin [\sigma (kx)].$$

As already mentioned, the potential (7) has rigid walls in $kx=\pm \pi$ and the wave function must be zero at these points. For even wave function, ${\varphi }_{even}(\pi )=0$, we have

(19) $$\sigma (\pi )=\frac{\pi }{2},\frac{3\pi }{2},\frac{5\pi }{2},\cdots$$

while for the odd wave function, ${\varphi }_{odd}(\pi )=0$, we have

(20) $$\sigma (\pi )=\frac{2\pi }{2},\frac{4\pi }{2},\frac{6\pi }{2},\cdots$$

so we assign

$$\sigma (\pi )=\frac{n\pi }{2}\text{and}\sigma (-\pi )=-\frac{n\pi }{2}\text{where}n=1,2,3,\cdots$$

It follows that

$$\sigma (\pi )-\sigma (-\pi )=\frac{1}{\hslash }{\displaystyle {\int }_{-\pi /k}^{\pi /k}p}(x)dx,$$

then it is straightforward to obtain

(21) $$\frac{1}{\hslash }{\displaystyle {\int }_{-\pi /k}^{\pi /k}p}(x)dx=n\pi \text{with}n=1,2,3,\cdots$$

Now one can write explicitly the potential (7) to $O(x^2)$

(22) $$V(x)\approx V_0\left(\frac{1}{3}+\frac{1}{45}{(kx)}^2\right),$$

where in the symmetrical interval $[-\pi ,\pi ]$ we have $O(1)\gg O(x^2)\gg O(x^4)$, such that the momentum (16) can be written as

$$p(x)\approx \sqrt{2m\left(E-\frac{1}{3}\frac{{\hslash }^2{k}^2}{m}\right)}\left(1-\frac{\left(\frac{{\hslash }^2{k}^2}{m}\right){(kx)}^2}{90\left(E-\frac{1}{3}\frac{{\hslash }^2{k}^2}{m}\right)}\right).$$

and the integral (21) can now be calculated

$$\frac{1}{\hslash }\sqrt{2m\left(E-\frac{1}{3}\frac{{\hslash }^2{k}^2}{m}\right)}\left[\frac{2\pi }{k}-\frac{(\frac{{\hslash }^2{k}^2}{m})\frac{{\pi }^3}{k}}{135\left(E-\frac{1}{3}\frac{{\hslash }^2{k}^2}{m}\right)}\right]=n\pi ,$$

it follows that

(23) $$E_n=\left[\frac{1}{8}\left(\frac{4{\pi }^2}{135}+\frac{n^2}{2}+\sqrt{\frac{4{\pi }^2}{135}{n}^2+\frac{n^4}{4}}\right)+\frac{1}{3}\right]\frac{{\hslash }^2{k}^2}{m}.$$

where we put the label n in $E_n$ to indicate that the energy depends on n. Note that for $n\to \infty$, $E_n~n^2$ as we can expect for a particle in a box. This approximation gives a good idea about the energy levels although the corrections do not differ significantly from the energy levels of the particle in a box. For example, consider $n=1$ in Eq. (23). One obtains $E_1\approx 0.5244{\hslash }^2{k}^2/m$, while for a particle in a box with an energy shift of $V_0/3$ and a length a gives $0.4583{\hslash }^2{k}^2/m$. For $n=8$ one obtains $E_8=8.406{\hslash }^2{k}^2/m$ and $8.333{\hslash }^2{k}^2/m$ for a particle in a box.

It is straightforward and important to notice that the ground state energy, for $n=1$, is approximately $E=0.5\frac{{\hslash }^2{k}^2}{m}$ as seen in Eq. (14). The WKB approach indicates the approximate values of energy of a quantum system, but as we have seen it is necessary to solve integrals like equation (21), with $p(x)=\sqrt{2m[E-V(x)]}$, where there are some difficulties to calculate it when the potential is a function expanded in a power series.

In fact, one must employs the WKB approach calculating correctly the turning points for each energy, which is a tough task once theses points depend on the energy, and the condition (21) does not apply anymore. But this will be left to the section results only for comparison reasons. For now, the resulting equation (23) gives us a good idea about the energy levels for our problem. In the next section we shall see the Stationary Perturbation Theory where we obtain approximated values of energy and the wave functions for some specific cases.

4. Stationary Perturbative Theory

We have already encountered from equations (7) and (8) that the expansion of potential contains terms like $x^2$, $x^4$, and so on. Thus, it is useful to make an approach with the Stationary Perturbation Theory to study the effect of the perturbed potential with $x^4$ and $x^6$ on the energy levels into the one-dimensional harmonic oscillator Hamiltonian [², ⁴- ⁶]. As we shall see, the Stationary Perturbation Theory is also a good technique to compute the first levels of energy and the respective wave functions ^[3]. However, the calculations can be very tough as soon as we go to higher levels of energy. Needless to say that the respective wave functions are also cumbersome to compute. This is due to the fact that the potential (7) is an infinite series and we can only probe some terms of the series. To circumvent this inconvenience one can resort to numerical techniques as we will see in the next section. Nevertheless, it remains important to understand the basic physics of the approximate solutions even before we explore a more accurate numerical method. Let us write the potential well (7) expanded to sixth order,

$$V(x)\approx V_0\left(\frac{1}{3}+\frac{1}{45}{(kx)}^2+\frac{2}{945}{(kx)}^4+\frac{1}{4725}{(kx)}^6\right),$$

as we already set $V_0=\frac{{\hslash }^2{k}^2}{m}$, we write the unperturbed Hamiltonian as

(24) $$-\frac{{\hslash }^2}{2m}\frac{d^2\psi }{d{x}^2}+\frac{1}{45}\frac{{\hslash }^2{k}^2}{m}{(kx)}^2\psi =H_0\psi$$

where $H_0$ is the unperturbed Hamiltonian, which is the Hamiltonian for an 1d Harmonic Oscillator with $\frac{1}{45}\frac{{\hslash }^2{k}^4}{m}=\frac{1}{2}m{\Omega }^2$.

For the perturbing potential, W, we write

(25) $$W(x)\psi =\left(\frac{1}{3}+\frac{2}{945}{(kx)}^4+\frac{1}{4725}{(kx)}^6\right)\frac{{\hslash }^2{k}^2}{m}\psi .$$

Using the transformation $\xi =kx$ we now write the 1d Schödinger equation as

(26) $$-\frac{d^2\varphi (\xi )}{d{\xi }^2}+\frac{1}{2}{\omega }^2{\xi }^2\varphi (\xi )+\left(\frac{1}{3}+\frac{2{\xi }^4}{945}+\frac{{\xi }^6}{4725}\right)\varphi (\xi )=\epsilon \varphi (\xi )$$

where we have $\omega =\frac{m\Omega }{\hslash k^2}\equiv \sqrt{\frac{2}{45}}$ and $\epsilon =\frac{mE}{{\hslash }^2{k}^2}$ dimensionless constants.

Now, it is convenient to define the creation and annihilation operators respectively [⁶,⁹],

(27) $$a^{†}=\sqrt{\frac{\omega }{2}}\left(\xi -\frac{ip}{\omega }\right),a=\sqrt{\frac{\omega }{2}}\left(\xi +\frac{ip}{\omega }\right),$$

where we define the operator $p=-i\frac{d}{d\xi }$. These operators obey the rules below

(28) $$a^{†}a=N,[a,a^{†}]=1,[N,a^{†}]=a^{†},[N,a]=-a,$$

also we have

(29) $$a|n\rangle =\sqrt{n}|n-1\rangle \text{and}{a}^{†}|n\rangle =\sqrt{n+1}|n+1\rangle .$$

The unperturbed Hamiltonian of the one-dimensional harmonic oscillator yields energy values

(30) $${\epsilon }_n^{(0)}=\left(n+\frac{1}{2}\right)\omega ,$$

with ${\epsilon }_n^{(0)}\equiv \frac{m{E}^0}{{\hslash }^2{k}^2}$, where the Stationary Perturbation Theory will add corrections to (30).

It is straightforward to obtain the coordinate operators using equations (27)

(31) $$\xi =\sqrt{\frac{1}{2\omega }}\left(a+a^{†}\right),$$

and with the aid of commutation relations (28) we can get

(32) $$\begin{array}{ll}{\xi }^4={\left(\frac{1}{2\omega }\right)}^2\hfill & [a^4+a^{†4}+(4N+6)a^2+(4N-2)a^{†2}\hfill \\ \hfill & +(6N^2+6N+3)]\hfill \end{array}$$

and

(33) $$\begin{array}{l}{\xi }^6={\left(\frac{1}{2\omega }\right)}^3[a^6+a^{†6}+(6N+15)a^4+(6N-9)a^{†4}\hfill \\ +(15N^2+45N+45)a^2+(15N^2-15N+15)a^{†2}\hfill \\ +(20N^3+30N^2+40N+15)].\hfill \end{array}$$

According to Stationary Perturbation Theory, the first-order correction to the energy is simply equals the mean value of perturbation term W in the unperturbed state $|n\rangle$,

(34) $${\epsilon }_n^{(1)}=\langle n\left|W\right|n\rangle .$$

where we set the perturbing term as

(35) $$W=\frac{1}{3}+\frac{2{\xi }^4}{945}+\frac{{\xi }^6}{4725}$$

We can see that the only non-vanishing terms in the first-order correction to the energy in Eq. (34)(when we use the equations (32) and (33)) are $6N^2+6N+3$ and $20N^3+30N^2+40N+15$, such as

$$\begin{array}{lll}{\epsilon }_n^{(1)}\hfill & =\hfill & \frac{1}{3}+\frac{2}{945}{\left(\frac{1}{2\omega }\right)}^2(6n^2+6n+3)\hfill \\ \hfill & \hfill & +\frac{1}{4725}{\left(\frac{1}{2\omega }\right)}^3(20n^3+30n^2+40n+15).\hfill \end{array}$$

Replacing ${\omega }^2=\frac{2}{45}$ into the above equation we can write

(36) $$\begin{array}{lll}{\epsilon }_n^{(1)}\hfill & =\hfill & \frac{1}{3}+\frac{1}{84}(6n^2+6n+3)\hfill \\ \hfill & \hfill & +\sqrt{\frac{45}{2}}\left(\frac{1}{1680}\right)(20n^3+30n^2+40n+15).\hfill \end{array}$$

As we mentioned elsewhere, the calculations for higher levels in Stationary Perturbation Theory is a tough task with a great quantity of integral calculations. There is no reason to go beyond once several accurate numerical techniques are available to solve the Schödinger equation. As a consequence, we are limited to the first levels of energy where the first terms in the potential expansion are important. For this reason we will disregard corrections due to the sixth-order term ${\xi }^6$, we will simplify the calculations to second-order of correction to the energy using just only the term ${\xi }^4$ of (35) for both energy calculations and the respective wave functions, keeping this contribution only in the first-order correction as shown in Eq. (36). We therefore have

(37) $${\epsilon }_n^{(2)}={\displaystyle \sum _{n\ne n^{\prime }}\frac{|\langle n|W_{{\xi }^4}|n^{\prime }\rangle {|}^2}{{\epsilon }_{n^{\prime }}^{(0)}-{\epsilon }_n^{(0)}}}.$$

The terms of $W_{{\xi }^4}=\frac{2{\xi }^4}{945}$ that contribute to second-order are ${\left(\frac{1}{2\omega }\right)}^2[a^4+a^{†4}+(4N+6)a^2+(4N-2)a^{†2}]$, so that it follows

(38) $$\begin{array}{l}\langle n-2\left|W_{{\xi }^4}\right|n\rangle =\frac{1}{84}\left[(4n-2)\sqrt{n(n-1)}\right]\hfill \\ \langle n+2\left|W_{{\xi }^4}\right|n\rangle =\frac{1}{84}\left[(4n+6)\sqrt{(n+1)(n+2)}\right]\hfill \\ \langle n-4\left|W_{{\xi }^4}\right|n\rangle =\frac{1}{84}\left[\sqrt{n(n-1)(n-2)(n-3)}\right]\hfill \\ \langle n+4\left|W_{{\xi }^4}\right|n\rangle =\frac{1}{84}\left[\sqrt{(n+1)(n+2)(n+3)(n+4)},\right]\hfill \end{array}$$

The values of the matrix elements $\langle n\left|W_{{\xi }^4}\right|n^{\prime }\rangle$ necessary to compute the energy corrections as well as the respective wave functions are listed in Appendix in the supplementary material. In what follows, when we calculate the second-order correction for the ground state energy we find

(39) $$\begin{array}{lll}{\epsilon }_0^{(2)}\hfill & =\hfill & \frac{|\langle 2|W_{{\xi }^4}|0\rangle {|}^2}{{\epsilon }_0^{(0)}-{\epsilon }_2^{(0)}}+\frac{|\langle 4|W_{{\xi }^4}|0\rangle {|}^2}{{\epsilon }_0^{(0)}-{\epsilon }_4^{(0)}}\hfill \\ \hfill & =\hfill & \frac{\frac{2}{196}}{(-2)\sqrt{\frac{2}{45}}}+\frac{\frac{6}{1764}}{(-4)\sqrt{\frac{2}{45}}}\hfill \\ \hfill & =\hfill & -\mathrm{0.028234621.}\hfill \end{array}$$

We can continue with this method and obtain some values for second-order correction values such as

(40) $$\begin{array}{llllll}{\epsilon }_1^{(2)}\hfill & =\hfill & -0.2218,10,\hfill & {\epsilon }_2^{(2)}\hfill & =\hfill & -0.8269\hfill \\ {\epsilon }_3^{(2)}\hfill & =\hfill & -2.1176,10,\hfill & {\epsilon }_4^{(2)}\hfill & =\hfill & -4.3683\hfill \\ {\epsilon }_5^{(2)}\hfill & =\hfill & -7.8533,10,\hfill & {\epsilon }_6^{(2)}\hfill & =\hfill & -12.8468\hfill \\ {\epsilon }_7^{(2)}\hfill & =\hfill & -19.6231\hfill & \hfill & \hfill & \hfill \end{array}$$

Therefore we can compute the energy values using (30), (36), (39), and (40) where we have

(41) $${\epsilon }_n={\epsilon }_n^{(0)}+{\epsilon }_n^{(1)}+{\epsilon }_n^{(2)}$$

with ${\epsilon }^{(i)}=\frac{m{E}^{(i)}}{{\hslash }^2{k}^2},i=0,1,2$.

Now we move on to obtain the approximate wave function. The first-order correction of the state vector is a linear superposition of all the unperturbed states

(42) $${\psi }_n^{(1)}(x)={\psi }_n^{(0)}(x)+{\displaystyle \sum _{n\ne k}\frac{\langle {\psi }_k^{(0)}\left|W\right|{\psi }_n^{(0)}\rangle }{E_n^{(0)}-E_k^{(0)}}}{\psi }_k^{(0)}.$$

The eigenfunctions of the unperturbed one-dimensional harmonic oscillator are

(43) $${\psi }_n^{(0)}(u)={\left(\frac{m\omega }{\pi \hslash }\right)}^{1/4}\frac{1}{\sqrt{2^{n}n!}}{H}_n(u)e^{-u^2/2},$$

where $u=\sqrt{\frac{m\omega }{\hslash }}x\equiv {\left(\frac{2}{45}\right)}^{1/4}\xi$ and $H_n(u)$ are the Hermite polynomials,

(44) $$\{\begin{array}{l}{H}_0=1\hfill \\ H_1=2u\hfill \\ H_2=4u^2-2\hfill \\ H_3=8u^3-12u\hfill \\ H_4=16u^4-48u^2+12\hfill \\ H_5=32u^5-160u^3+120u.\hfill \end{array}$$

Thus, we can calculate the ground state wave function using equation (42), it follows that

$$\begin{array}{l}{\psi }_0^{(1)}={\psi }_0^{(0)}+\frac{\langle {\psi }_2^{(0)}\left|W_{{\xi }^4}\right|{\psi }_0^{(0)}\rangle }{E_0^{(0)}-E_2^{(0)}}{\psi }_2^{(0)}\hfill \\ +\frac{\langle {\psi }_4^{(0)}\left|W_{{\xi }^4}\right|{\psi }_0^{(0)}\rangle }{E_0^{(0)}-E_4^{(0)}}{\psi }_4^{(0)}\hfill \\ ={\psi }_0^{(0)}+\frac{\langle 2\left|W_{{\xi }^4}\right|0\rangle }{E_0^{(0)}-E_2^{(0)}}{\psi }_2^{(0)}+\frac{\langle 4\left|W_{{\xi }^4}\right|0\rangle }{E_0^{(0)}-E_4^{(0)}}{\psi }_4^{(0)}\hfill \\ ={\psi }_0^{(0)}+\frac{\frac{\sqrt{2}}{14}}{(-2)\sqrt{\frac{2}{45}}}{\psi }_2^{(0)}+\frac{\frac{\sqrt{6}}{42}}{(-4)\sqrt{\frac{2}{45}}}{\psi }_4^{(0)}\hfill \\ ={\psi }_0^{(0)}-0.239578711{\psi }_2^{(0)}-0.069160416{\psi }_4^{(0)},\hfill \end{array}$$

we can now put the values of unperturbed functions of (43) into above result and obtain

(45) $$\begin{array}{lll}{\psi }_0^{(1)}\hfill & =\hfill & {\left(\frac{1}{\pi }\sqrt{\frac{2}{45}}\right)}^{1/4}{e}^{-u^2/2}[1-0.084703865(4u^2-2)\hfill \\ \hfill & =\hfill & -0.003529328(16u^4-48u^2+12)].\hfill \end{array}$$

In the same way we can obtain the first excited wave function as fallows

(46) $$\begin{array}{l}{\psi }_1^{(1)}={\left(\frac{1}{\pi }\sqrt{\frac{2}{45}}\right)}^{1/4}{e}^{-u^2/2}[1.414213562u\hfill \\ -0.099824463(8u^3-12u)\hfill \\ -0.002495612(32u^5-160u^3+120u)].\hfill \end{array}$$

where these two wave functions are not normalized yet. As we already mentioned, the Perturbative Method is unproductive to obtain accurate results because it is necessary many calculation that converge slowly. For this problem with potential well (7), we will employ an efficient numerical method which can be implemented in any computer language.

5. Numerical Method

For simplicity, we will write the 1D Schrödinger Equation exactly like we did in Eq. (10), where ${\upsilon }_0\equiv \frac{m{V}_0}{{\hslash }^2{k}^2}$ and $\epsilon \equiv \frac{mE}{{\hslash }^2{k}^2}$. Unlike what we did in the last sections, we will use from now on x as dimensionless variable instead of $\xi$ for convenience.

(47) $$-\frac{1}{2}\frac{d^2}{d{x}^2}\psi +\upsilon (x)\psi =\epsilon \psi$$

The Numerov Method ^[7] can always be applied if the differential equation is in the format

(48) $$\frac{d^2\psi }{d{x}^2}=-g(x)\psi +s(x)$$

Then we want to solve the equation (47) when the potential $V(x)$ is given by

(49) $$\upsilon (x)={\upsilon }_0\left(\frac{1-x\cot (x)}{x^2}\right)$$

we will rewrite it as

(50) $$\frac{d^2\psi }{d{x}^2}=-2\left[\epsilon -\upsilon (x)\right]\psi$$

and if we make $g(x)=2\left[\epsilon -\upsilon (x)\right]$ then

(51) $$\frac{d^2\psi }{d{x}^2}={\psi }^{″}=-g(x)\psi$$

in this case the function $s(x)$ appearing in Eq. (48) is zero. With this in mind, $x_k=-\pi +hk$, and $k=0,1,2,\cdots ,N$. Also assume that at $x_k$ we have ${\psi }_k=\psi (x_k)$. Since Eq. (51) is an homogeneous linear differential equation. Therefore, if ${\phi }_1(x)$ is solution to the equation, then is as well ${\phi }_2(x)=A{\phi }_1(x)$, where A is some real constant.

It is worth to mention the special case when $A=-1$ then ${\phi }_2(x)=-{\phi }_1(x)$. This fact will be important because in the Numerov Method we need two consecutive values for the solution, ${\psi }_k$ and ${\psi }_{k+1}$. So let us assume we start with ${\Psi }_0=0$ then we can guess ${\Psi }_1=\pm ϵ$.

From the Boundary Conditions, at $x_0=-\pi$ we must have ${\Psi }_0=0$, in other words, this point $(-\pi ,0)$ has to be a value for the solution. Guessing a positive or negative value for ${\Psi }_1$ will only drive us to $+\Psi (x)$ or $-\Psi (x)$ then it doesn't matter if we guess a positive or negative variation value in the neighborhood of ${\psi }_0$, since we will get either $\Psi$ or $-\Psi$ but both are solutions. Assume that the amplitudes of some solution $\psi$ when x is discretized in the interval $[-\pi ,\pi ]$ be given by $\psi =\{{\psi }_0,{\psi }_1,{\psi }_2,\cdots ,{\psi }_N\}$. Accordingly, $\Psi =A\psi =\{A{\psi }_0,A{\psi }_1,A{\psi }_2,\cdots ,A{\psi }_N\}$ constitute also a set of amplitudes which are also amplitudes for another solution of the Schrödinger Equation (47).

We are free to pick any value for constant A, we can always adjust it such that $A{\psi }_1=ϵ$, for instance, $ϵ=0.0001$, or in other words: $A={10}^{-4}/{\psi }_1$. That is reasonable-since ${\Psi }_0=0$ and because the solution is assumed to be a smooth function-then we won't expect huge variations from one value to the ones in some neighborhood. If we can calculate the next value ${\Psi }_2$ based on ${\Psi }_0$ and ${\Psi }_1$ we can repeat the process and find all other values of ${\Psi }_k$. We can use the Numerov's Method to find ${\Psi }_{k+1}=f({\Psi }_k,{\Psi }_{k-1})$, which is derived in the Appendix ??. This method has order 4 of convergence ^[7].

The selection goes as following, for those values of E, starting with ${\Psi }_0=0$ and ${\Psi }_1=ϵ$ which finish with ${\Psi }_N\approx 0$ will constitute a set of valid values for the solution of the Schrödinger equation, in this way complying with the two Boundary Conditions in $\pm \pi$ or $\Psi (\pi )=\Psi (-\pi )=0$.

6. Results

Unlike what was done in section 3 we must implement the full WKB method to calculate the energies and the respective wave functions ^[3]. There we approximated the potential to $O(x^2)$ and derived an expression for $E_n~n^2$. As a consequence the energy levels started from the value of $n=1$. This occurred due to the considerations which led to the expression (21). To obtain the energies and the wave functions in WKB for the potential (3) we must find the turning points for each energy. This is an additional difficulty because the turning points define the energy. However, one can circumvent this difficulty by applying the following WKB constraint

(52) $$\frac{1}{\hslash }{\displaystyle {\int }_{x_1}^{x_2}\sqrt{2m[E-V(x)]}}dx=\left(n+\frac{1}{2}\right)\pi \text{with}n=0,1,2,\cdots$$

where $x_1\equiv x_1(E)$ and $x_2\equiv x_2(E)$ are the turning points. To obtain the energy and the respective turning points we passed the right hand side in Eq. (52) to the left hand side, defining a new function. Now the task is to determine the respective values of the energy which makes this new function equals zero, i.e., for a given level n one must only find roots, remembering that the turning points are also dependent of the energy. We did that in Scilab by means of Newton-Raphson algorithm to find roots, varying the energy with 0.0001${\hslash }^2{k}^2/m$ for each step, calculating the respective turning points making $E=V(x)$, which in turn are used to calculate the integral in Eq. (52). Thus, using the approximated expression for the energy levels in section 3 we could assign an initial value bellow these values to find the energy which minimizes this new function for a given n. The resulting energies can be visualized in first column in Table 1. These are the energy values for the full WKB implementation considering the turning points in each level of energy. We will not describe the computation of the respective wave functions here because it is well-established in the text books [³,⁴].

Table 1
Energies in units of $\frac{{\hslash }^2{k}^2}{m}$ calculated using WKB, Stationary Perturbation Theory (SPT) and Numerov's numerical approach for the 1d Schrödinger equation with the suggested potential V(x).

Level n	WKB	SPT	Numerov
0	0.4579	0.4886	0.5000
1	0.8408	0.9028	0.9424
2	1.4307	1.5566	1.6108
3	2.2465	2.5143	2.5175
4	3.2976	3.8414	3.6672
5	4.5890	5.6014	5.0628
6	6.1237	7.8593	6.7055
7	7.9035	10.6795	8.5963

Now we call attention to the first level of Numerov's algorithm calculation. The value of 0.500 ${\hslash }^2{k}^2/m$ is exactly the ground state value obtained analytically in section 3. This method produces good results in solving 1d-Schrödinger equation like pointed out by F. Caruso and V. Oguri ^[8]. The energy results are summarized in the last column and will serve to compare with other approximated methods. As one can see, the WKB method gives the best results if compared with Stationary Perturbation Theory (SPT) for the higher levels. One of the reasons that the SPT gives some discrepancies is the fact that we only go up to $O(x^4)$ in the potential series to the second-order corrections as we mentioned elsewhere. However, these discrepancies are only perceptive after $n=4$, giving reasonable results for the energy levels bellow that. As expected, the WKB method approximates to the exact values as soon as n becomes big.

In figure 4 we plot $|\psi {|}^2$ for some wave functions generated by each method. In figure 4 (a) we plot the analytic expression for the ground state wave function from section 3 with that obtained using SPT and the Numerov methods. As one can see, Numerov is quite accurate while the SPT, although not so accurate, seems to be a reasonable approximation. However, we cannot say the same in figure 4 (b) for the first excited state. The SPT result in comparison with Numerov's solution is quite different due to the increase in the difference of the energies between the methods as can be seen in Table 1. Finally, in figures 4 (c) and (d) we plot the WKB results for $n=6$ and $n=7$ with the respective Numerov's results. We can see that there is a difference in amplitude and in the phase which can be explained by the difference in energy as can be seen in Table 1.

Figure 4
(a) is the plot for $|\psi {|}^2$ of the ground state wave function using the analytic expression, and SPT and Numerov methods. (b) is the plot for $|\psi {|}^2$ of the first excited state for both SPT and Numerov methods. (c) is the plot for $|\psi {|}^2$ of the sixth excited state for both WKB and Numerov methods. (d) is the plot for $|\psi {|}^2$ of the seventh excited state for both WKB and Numerov methods.

7. Conclusion

We have suggested a mathematical potential as a good exercise in employing different Quantum Mechanical techniques to solve the 1d-Schödinger equation. We applied both WKB and Stationary Perturbation Theory to get the energy states of the 1d-Schödinger equation and compared them with an accurate numerical solution described elsewhere. One can explore other applications of such a potential to other physical situations like to solve the 3d-Schödinger equation for certain conditions and even to solids. Other possibilities will be explored in future works.

Supplementary material

The following online material is available for this article:

Appendix - Perturbative Matrix Elments

References

• ^[1]
  S. Weinberg, Lectures on Quantum Mechanics (Cambridge University Press, Cambridge, 2013).
• ^[2]
  J.J. Sakurai, Modern Quantum Mechanics (Addison-Wesley Publishing Company, Boston, 1995).
• ^[3]
  J.R.P. Mahon, Mecânica Quântica: Desenvolvimento Contemporâneo com Aplicações (LTC, São Paulo, 2011).
• ^[4]
  D. Griffths, Introduction to Quantum Mechanics (Cambridge University Press, Cambridge, 2016), 2ª ed.
• ^[5]
  L.I. Schiff, Quantum Mechanics (McGraw-Hill Book Company, New York, 1968).
• ^[6]
  C. Cohen-Tannoudji, B. Diu and F. Laloë, Quantum Mechanics (John Wiley & Sons, New Jersey, 1977), v. 2.
• ^[7]
  E. Hairer, S.P. Nørsett and G. Wanner, Solving ordinary differential equations I: Nonstiff problems (§III.10) (Springer-Verlag, Berlin, 1993).
• ^[8]
  F. Caruso and V. Oguri, Revista Brasileira de Ensino de Fisica 36, 2310 (2014).
• ^[9]
  R.P. Feynman, Statistical Mechanics-A Set of Lectures (Westview Press, Colorado, 1998).

Publication Dates

• Publication in this collection
  25 Apr 2019
• Date of issue
  2019

History

• Received
  28 Dec 2018
• Accepted
  05 Mar 2019