Applications of the Numerov method to simple quantum systems using Python

Abstract

Numerov’s numerical method is developed in a didactic way by using Python in its Jupyter Notebook version 6.0.3 for three different quantum physical systems: the hydrogen atom, a molecule governed by the Morse potential and a quantum dot. After a brief introduction to the Numerov method, the complete code to calculate the eigenfunctions and eigenvalues of the hydrogen atom is presented. The necessary code changes to compute the other two examples are also provided in the sequel.

Keywords:
Hydrogen Atom; Morse Potential; Quantum Dot; Numerov method; Python

1. Introduction

In the last decades, the teaching of computing techniques has become more present and also increasingly essential in the development of students from all areas, and, therefore, it would not be different for physics teaching. Having said that, it is of paramount importance that we always produce new teaching materials for new technologies such as the Python programming language, which, despite of being relatively new, has already dominated the market becoming one of the most important languages today.

We will disclose here a powerful numerical calculation method originally developed by Boris Vasil’evich Numerov [¹,²,³] – see also [⁴,⁵,⁶,⁷] – applying it to the time-independent Schrödinger equation describing physical systems like the hydrogen atom, a diatomic molecule governed by the Morse potential and a particular model for the quantum dot atom [⁸,⁹,¹⁰]. These three examples will be solved and the parameters needed for each solution using Numerov’s method will be shared in their respective sections.

In summary, this work aims to provide the complete code developed in Python with the Jupyter Notebook for Numerov’s numerical method. However, it is important to emphasize that we do not aim to teach Python to the reader, who must have a basic knowledge of programming to be able to keep up with the examples.

The code used in this article can be found in full and downloaded at the following link: https://1drv.ms/u/s!Ai_Lqkgh1kiskp5_cfOtwCfqX-LpTw?e=srdmzS, as well as in our repository on Github https://github.com/Felipecrdss/Numerov.

2. Numerov’s Method

Numerov’s initial motivation was to be able to calculate corrections to the trajectory of Halley’s comet. Therefore, Numerov’s method was initially developed to determine solutions to eigenvalue problems associated with second order ordinary differential equations of celestial mechanics, which did not contain terms involving the first derivative of an unknown function $y(x)$, that is, equations of the form

(1) $$\frac{{\text{d}}^{2}y}{\text{d}{x}^2}=f(y,x).$$

Every differential equation equal to equation (1) can be replaced by the following system of first order equations

$$\{\begin{array}{cc}\frac{\text{d}z}{\text{d}x}=f(x,y),\hfill & \\ z=\frac{\text{d}y}{\text{d}x}.\hfill & \end{array}$$

Traditional methods for numerically solving this system of equations, such as those of Euler or Runge-Kutta, consider that the values of $y(x)$ and of $\text{d}y/\text{d}x$ are known at a given point in the domain $[a,b]$ of system validity, i.e., are suitable for the so-called seed problems.

In non-relativistic quantum mechanics, more specifically in bound state problems involving a particle of mass m confined in a well of potential $V(x)$, in a given interval $a<x<b$, the allowed energies $(E)$ and the corresponding wave functions $\psi (x)$ that describe these steady states satisfy the Schrödinger’s eigenvalue equation

(2) $$\frac{{\text{d}}^2\psi }{\text{d}{x}^2}+k^2(x)\psi =0$$

where $k=\sqrt{2m[E-V(x)]}/\mathrm{\hslash }$ and $\mathrm{\hslash }\simeq 1.055\times {10}^{-34}J.s$ is the reduced Planck constant.

In these cases, as the value of the first derivative of the wave function is not known, the Euler and Runge-Kutta cannot be employed. Nonetheless, it is possible to establish continuity conditions for the values of $\psi$ and $\text{d}\psi /\text{d}x$ at two or more points of the domain of the wave function, which characterizes the so-called boundary value problems.

In addition, making the transformation of a second order differential equation in a first order system, the Numerov method allows the simultaneous determination of the particle energy spectrum and the eigenfunctions associated with each energy value.

Like any iterative numerical method, the solution of equation (2) is constructed by successive integrations. In Numerov’s method, initially, the solution is considered to be known at two subsequent points of the interval $[a,b]$; for example, at $\psi (x-\delta )$ and $\psi (x)$, where $\delta$ is an arbitrarily small quantity, called the integration step. Next, we try to establish an algorithm to determine the solution at the next point, $\psi (x+\delta )$.

The starting point for establishing this algorithm is the expansion of $\psi (x\pm \delta )$ in Taylor series, up to fourth-order derivatives, that is,

(3) $$\psi (x\pm \delta )=\psi (x)\pm \delta {\psi }^{\prime }(x)$$ $$+{\displaystyle \frac{{\delta }^2}{2}}{\psi }^{\prime \prime }(x)\pm {\displaystyle \frac{{\delta }^3}{6}}{\psi }^{\prime \prime \prime }(x)$$ $$+{\displaystyle \frac{{\delta }^4}{24}}{\psi }^{iv}(x)$$

Adding the terms $\psi (x+\delta )$ and $\psi (x-\delta )$, only the derivatives of even order survive and, therefore, a relationship between the values of a function in three points and its second derivative is reached, given by

(4) $${\displaystyle \frac{\psi (x+\delta )+\psi (x-\delta )-2\psi (x)}{{\delta }^2}}$$ $$={\psi }^{\prime \prime }(x)+{\displaystyle \frac{{\delta }^2}{12}}{\psi }^{iv}(x)$$ $$\equiv \left(1+{\displaystyle \frac{{\delta }^2}{12}}{\displaystyle \frac{{\text{d}}^2}{\text{d}{x}^2}}\right){\psi }^{\prime \prime }x$$

Writing the unidimensional Schrödinger equation, equation (2), in a more convenient form

(5) $$\left(1+{\displaystyle \frac{{\delta }^2}{12}}{\displaystyle \frac{{\text{d}}^2}{\text{d}{x}^2}}\right){\psi }^{\prime \prime }(x)=-k^2(x)\psi (x)$$ $$-{\displaystyle \frac{{\delta }^2}{12}}{\displaystyle \frac{{\text{d}}^2}{\text{d}{x}^2}}\left[k^2(x)\psi (x)\right]$$

and using equation (4) to replace the terms that contain second order derivatives, one obtains

(6) $$\frac{\psi (x+\delta )+\psi (x-\delta )-2\psi (x)}{{\delta }^2}$$ $$=-k^2(x)\psi (x)-\frac{{\delta }^2}{12}$$ $$\times \left[\begin{array}{c}\hfill \frac{\begin{array}{c}{k}^2(x+\delta )\psi (x+\delta )\hfill \\ +k^2(x-\delta )\psi (x-\delta )\hfill \\ -2k^2(x)\psi (x)\hfill \end{array}}{{\delta }^2}\hfill \end{array}\right]+𝒪({\delta }^4)$$

Regrouping the terms, we obtain the Numerov difference formula for the problem of a particle under action of a one-dimensional potential:

(7) $$\left[1+{\displaystyle \frac{h^2}{12}}{k}^2(x+\delta )\right]\psi (x+\delta )$$ $$=2\left[1-{\displaystyle \frac{5{\delta }^2}{12}}{k}^2(x)\right]\psi (x)$$ $$-\left[1+{\displaystyle \frac{{\delta }^2}{12}}{k}^2(x-\delta )\right]\psi (x-\delta )$$

Since the problem of interest is an eigenvalue problem, the numerical integration technique of the one-dimensional Schrödinger equation for a particle in a well depends on attaching arbitrary values conveniently to eigenvalues and the respective (possible) eigenfunctions in 2 points of the domain of the problem. But how to do it? Regarding the choice of the initial value for the energy (first eigenvalue), just remember that, according to Heisenberg’s uncertainty relation, the energy E of a particle in a well of potential $V(x)$ must be greater than the minimum value of the well. Thus, it is considered, initially, that $E_{initial}=V_{min}+\mathrm{\Delta }E$, with $\mathrm{\Delta }E>0$.

The choice of an energy value, determines two turning points, x and x, where the energy value is equal to the potential energy value, whose motion obeys the classical Newtonian mechanics. That is, from the point of view of classical mechanics, the movement of the particle is restricted only to the region $[x_{\mathrm{\ell }},x_r]$, in which the energy is greater than or equal to the potential energy. The regions $x<x_{\mathrm{\ell }}$ and $x>x_r$ are called classically prohibited regions, and are indicated in Figure 1.

Figure 1
Meeting points between the potential curve and the initial energy, also called turning points.

As the Schrödinger equation admits solutions for these classically prohibited regions, for each energy value, initially, values are assigned to a possible eigenfunction at two points of the classically prohibited regions, in which the function practically cancels itself. In general, these are the boundary points a and $b>a$ of the function’s integration domain.

However, the implementation of Numerov’s method to solve the problem still requires an iteration scheme that uses the Numerov formula in two steps: from a, or to the left of a from the classic turning points, hereinafter called match point (x_match), and from b, or to the right of the match point.

Thus, arbitrarily taking an initial value for the energy, and two successive arbitrary values for the solution, starting from the lower extremes and upper part of the integration interval $[a,b]$, one can implement the method’s iteration scheme in the two senses, such as:

• •

  Solution to the left of match point $(x<x_{match})$

Being $E_{initial}=V_{min}+\mathrm{\Delta }E\left(\frac{\mathrm{\Delta }E}{|V_{min}|}\ll 1\right)$ an arbitrary value for the energy of the particle. Also arbitrating values for the wave function, in two successive points, from a,

$$\{\begin{array}{cc}{\psi }^{\mathrm{\ell }}(a)=0,\hfill & \\ {\psi }^{\mathrm{\ell }}(a+\delta )={\delta }^{\mathrm{\ell }},({\delta }^{\mathrm{\ell }}\ll 1)\hfill & \end{array}$$

and using the formula of differences, equation (7), the solution on the left is built sequentially until match point $(x_{match})$, in what ${\psi }^{\mathrm{\ell }}(x_{match})={\psi }_{match}^{\mathrm{\ell }}$.

• •

  Solution to the right of match point $(x>x_{match})$

In a similar way, for the same values E_initial, arbitrating

$$\{\begin{array}{cc}{\psi }^r(b)=0,\hfill & \\ {\psi }^r(b+\delta )={\delta }^r,({\delta }^r\ll 1)\hfill & \end{array}$$

the solution to the right, from b, is constituted sequentially until the points x and $x=x_{match}-\delta$, as

$$\{\begin{array}{cc}{\psi }^r(x_{match})={\psi }_{match}^r,\hfill & \\ {\psi }^r(x_{match}-\delta )={\psi }_{match-1}^r,({\delta }^r\ll 1)\hfill & \end{array}$$

To guarantee the boundary condition of the solution, we redefine the solution to the left according to equation (8) given below, and the boundary condition of the first derivatives, according to equation (9).

The procedure is repeated step by step, in the two ways, $a\rightleftharpoons b$. Starting from a, using the Numerov recurrence formula associated with an equation, if we build the solution ${\psi }^{\mathrm{\ell }}$ until the classic rewind point, nearest of b, where $E=V(x_{match})$, called the match point. Then, from b, the analog is made, building a solution ${\psi }^r$ to the match point. In principle, the possible solutions ${\psi }^{\mathrm{\ell }}$ and ${\psi }^r$ will not necessarily be equal in this x stitch. To ensure the continuity of the solution redefines itself ${\psi }^{\mathrm{\ell }}$ like

(8) $${\psi }^{\mathrm{\ell }}(x)\to {\psi }^{\mathrm{\ell }}(x)\frac{{\psi }^r(x_{match})}{{\psi }^{\mathrm{\ell }}(x_{match})}\mathit{}(a\le x\le x_{match})$$

Finally, it is verified how close are the values of the respective first derivatives of ${\psi }^r$ and the new function ${\psi }^{\mathrm{\ell }}$ It is staggered, at match point. To test the boundary condition of the derivatives first, taking into account Taylor’s series for $\psi (x+\delta$ and $\psi (x-\delta )$, up to the first order, you can write

(9) $$\{\begin{array}{cc}{\frac{\text{d}{\psi }^{\mathrm{\ell }}}{\text{d}x}|}_{x_r}=\frac{{\psi }_{match+1}^{\mathrm{\ell }}-{\psi }_{match-1}^{\mathrm{\ell }}}{2\delta },\hfill & \\ {\frac{\text{d}{\psi }^r}{\text{d}x}|}_{x_r}=\frac{{\psi }_{match+1}^r-{\psi }_{match-1}^r}{2\delta },\hfill & \end{array}$$

in what, ${\psi }_{match\pm 1}=\psi (x_{match\pm \delta })$.

If the difference between these values is less than the values of a predefined error, the process is interrupted, confirming the searched eigenvalue and the respective eigenfunction as being

$$\{\begin{array}{cc}{\psi }^{\mathrm{\ell }}(x),(a\le x<x_{match})\hfill & \\ {\psi }^r(x),(x_{match}\le x\le b)\hfill & \end{array}$$

If the continuity condition of the derivatives is not satisfied, the value of the energy is increased and we restarted a search for a new value which is really an eigenvalue of the problem, and its respective eigenfunction. The process can be repeated until the desired number of eigenvalues and eigenfunctions of the problem.

Because it is based on Taylor’s serial expansion to the fourth order, the error in Numerov’s method is much smaller than the errors that come out from the expansion-based methods in lower order, like that of Runge-Kutta.

3. Hydrogen Atom

Numerical solutions of hydrogen atom were previously obtained in [¹] using Numerov’s method. The program was written in C++ for the ROOT cint compiler.

Although it was originally developed for second order linear and homogeneous ordinary differential equations that do not contain terms of the first derivative, the Numerov’s method can be generalized to cover the presence of terms that contain the first derivative in the differential equation, so that eigenvalue problems can also be considered.

In fact, in the case of linear equations, every second order differential equation of the type

$$\frac{{\text{d}}^{2}y}{\text{d}{x}^2}+P(x)\frac{\text{d}y}{\text{d}x}+Q(x)y=0$$

can be written in its normal form

$$\frac{{\text{d}}^{2}y}{\text{d}{x}^2}+q(x)y=0$$

where

$$q(x)=Q(x)-\frac{1}{4}{P}^2(x)-\frac{1}{2}\frac{\text{d}P}{\text{d}x}$$

Schrödinger’s radial equation for a particle of mass m under the action of a Coulombic electric field, like the electron in the hydrogen atom, can be written as

(10) $${\displaystyle \frac{{\text{d}}^{2}R(r)}{\text{d}{r}^2}}+{\displaystyle \frac{2}{r}}{\displaystyle \frac{\text{d}R(r)}{\text{d}r}}$$ $$+{\displaystyle \frac{2m}{{\mathrm{\hslash }}^2}}\left[E+{\displaystyle \frac{e^2}{r}}-{\displaystyle \frac{{\mathrm{\hslash }}^2}{2m}}{\displaystyle \frac{\mathrm{\ell }(\mathrm{\ell }+1)}{r^2}}\right]R(r)=0$$

Making the substitution r = xa_B with $a_B={\mathrm{\hslash }}^2/(m{e}^2)$ being the Bohr radius, equation (10) can be rewritten, for a new function $y(x)=R(r)$, as

(11) $$\frac{{\text{d}}^{2}y}{\text{d}{x}^2}=-\frac{2}{x}\frac{\text{d}y}{\text{d}x}-\left[\epsilon -V(x)\right]y(x)$$

where, $\epsilon =\frac{E}{e^2/(2a_B)}$ and $V(x)=\frac{\mathrm{\ell }(\mathrm{\ell }+1)}{x^2}-\frac{2}{x}$ are, respectively, the energy and the so-called effective potential in atomic units. So, in possession of equation (11), we can start building our program code.

First of all, we must import the functions matplotlib.pyplot (for plotting) and NumPy (for Mathematics and working with arrays); also, to avoid any mistakes, we need to import the function random. We then declare the effective potential $V(x)$. During the construction of all the examples we will use the value $\mathrm{\ell }=0$.

The equation, in this case, that is intended to be solved by Numerov’s method presents a term involving the first derivative, and can be expressed by

(12) $${\psi }^{\prime \prime }(x)=-p(x){\psi }^{\prime }(x)-s(x)\psi (x),$$

where

$$\{\begin{array}{cc}p(x)=\frac{2}{x}\Rightarrow p^{\prime }(x)=-\frac{2}{x^2},\hfill & \\ s(x)=\epsilon -V(x).\hfill & \end{array}$$

From the Taylor expansion, equation (3), we can rewrite equation (12) as

(13) $$\left(1+{\displaystyle \frac{{\delta }^2}{12}}{\displaystyle \frac{{\text{d}}^2}{\text{d}{x}^2}}\right){\psi }^{\prime \prime }(x)=-p(x){\psi }^{\prime }(x)-s(x)\psi (x)$$ $$-{\displaystyle \frac{{\delta }^2}{12}}{\displaystyle \frac{{\text{d}}^2}{\text{d}{x}^2}}\left[p(x){\psi }^{\prime }(x)+s(x)\psi (x)\right].$$

In a similar way to the previous case, according to equation (4), you can write the term on the right side of the equation (13) which contains derivatives of order 2, such as

$${\displaystyle \frac{{\text{d}}^2}{\text{d}{x}^2}}\left[p(x){\psi }^{\prime }(x)+s(x)\psi (x)\right]$$ $$={\displaystyle \frac{1}{{\delta }^2}}[p(x+\delta ){\psi }^{\prime }(x+\delta )+s(x+\delta )\psi (x+\delta )$$ $$+p(x-\delta ){\psi }^{\prime }(x-\delta )+s(x-\delta )\psi (x-\delta )$$ $$+-2p(x){\psi }^{\prime }(x)-2s(x)\psi (x)]$$

Replacing first order derivatives with approximations

$$\{\begin{array}{c}{\psi }^{\prime }(x)=\left[\psi (x+\delta )-\psi (x-\delta )\right]/(2\delta ),\hfill \\ {\psi }^{\prime }(x+\delta )=\left[\psi (x+\delta )-\psi (x)\right]/(\delta ),\hfill \\ {\psi }^{\prime }(x-\delta )=\left[\psi (x)-\psi (x-\delta )\right]/(\delta ),\hfill \end{array}$$

we obtain

$${\displaystyle \frac{{\text{d}}^2}{\text{d}{x}^2}}\left[p(x){\psi }^{\prime }(x)+s(x)\psi (x)\right]$$ $$={\displaystyle \frac{1}{{\delta }^2}}\{[{\displaystyle \frac{p(x+\delta )-p(x)}{\delta }}+s(x+\delta )]\times \psi (x+\delta )$$ $$+\left[{\displaystyle \frac{p(x)-p(x-\delta )}{\delta }}+s(x-\delta )\right]\psi (x-\delta )$$ $$+2[{\displaystyle \frac{p(x-\delta )-p(x+\delta )}{2\delta }}+s(x)]\psi (x)\}$$

or

(14) $$\frac{{\text{d}}^2}{\text{d}{x}^2}\left[p(x){\psi }^{\prime }(x)+s(x)\psi (x)\right]$$ $$=\frac{1}{{\delta }^2}\{[p^{\prime }(x)+s(x+\delta )]\psi (x+\delta )$$ $$+[p^{\prime }(x)+s(x-\delta )]\psi (x-\delta )$$ $$-2[p^{\prime }(x)+s(x)]\psi (x)\}$$

Taking into account that the left side of the equation (13) is equal to

$$[\psi (x+\delta )+\psi (x-\delta )-2\psi (x)]/{\delta }^2$$

we can write

$${\displaystyle \frac{\psi (x+\delta )+\psi (x-\delta )-2\psi (x)}{{\delta }^2}}$$ $$=-p(x)\left[{\displaystyle \frac{\psi (x+\delta )-\psi (x-\delta )}{2\delta }}\right]$$ $$-s(x)\psi (x)+{\displaystyle \frac{1}{12}}\left[p^{\prime }(x)+s(x+\delta )\right]\psi (x+\delta )$$ $$-{\displaystyle \frac{1}{12}}\left[p^{\prime }(x)+s(x-\delta )\right]\psi (x-\delta )$$ $$+{\displaystyle \frac{1}{6}}\left[p^{\prime }(x)+s(x)\right]\psi (x)$$

Regrouping the terms, and making

$$\{\begin{array}{c}\psi (x-\delta )={\psi }_0,\hfill \\ r\psi (x)={\psi }_1,\hfill \\ \psi (x+\delta )={\psi }_2,\hfill \end{array}$$

one obtains the Numerov difference equation for the problem, suitable for the propagation of the solution from the integration limits interval:

(15) $${\psi }_2={\displaystyle \frac{\begin{array}{c}2\left\{1-\left[s(X)-\frac{p^{\prime }(x)}{5}\right]\frac{5{\delta }^2}{12}\right\}{\psi }_1\hfill \\ -\left\{1-p(x)\frac{\delta }{2}+\left[s(x-\delta )+p^{\prime }(x)\right]\frac{{\delta }^2}{12}\right\}{\psi }_0\hfill \end{array}}{\left\{1+p(x)\frac{\delta }{2}+\left[s(x+\delta )+p^{\prime }(x)\right]\frac{\delta }{12}\right\}}}$$

From this formula, a procedure analogous to the previous case can be implemented for the construction of solutions of the radial Schrödinger equation in the interval $(0,\mathrm{\infty })$.

Now, in order to introduce the Numerov difference formula (15), we first need to insert equation (7) in our code, for that, let’s break it down into different pieces p0, p1 and p2, with $h=\delta$, $q_0=s(x)$. Thus, equation (15) is now called y2, where ${\psi }_i=y_i$.

And, for the case $x-\delta$, we repeat the previous step, changing the appropriate sign.

We now create a list for each variable up to the value of dim or nmax which will also be defined in a future step.

In the next steps, the program will determine the interaction for the eigenvalue candidates and determine their solutions and will print the parameters found on the screen.

So far, the only thing we need to change in our program is the equation that defines our effective potential $V(x)$. The rest of the code will be the same for any equation that is in the same form as equation (9). From now on, we must change the code whenever we are looking for solutions with a different potential.

The parameter M indicates the eigenvalue that we are determining, in the next step, in order to avoid that the program needs to sweep the entire potential well in search of solutions, we can give increments between one eigenvalue and another in the form of multiples of the $\text{d}E$ parameter, which will be duly defined in a next step.

Usually, this process involves trial and error, where we adjust the $\text{d}E$ multiplier for each case, until we find the desired eigenvalue. However, as we are developing this first example for the hydrogen atom, which already has its energy eigenvalues well defined, this task becomes much easier. Before we find the ground state (M=0) we must use the step $24\times \text{d}E$. Thus, we find for ground state energy the value of $\epsilon =-0.99453$, which must be compared with the known value of the ground state energy for the hydrogen atom $\epsilon =-1$. And the next multipliers so that we can find at least the first three solutions which are:

If you are interested in finding other solutions, you should add the steps for large values of M.

We finally reach the part of the program where we must introduce the initial values to proceed with the solution. Whenever we start to develop a new code, these must be the first values to be changed, right after choosing the effective potential $V(x)$.

The parameters a and b correspond to the boundary points mentioned in the second section. At first, our wave function should have its initial value (a) equal to 0; thus, to avoid divisions by 0 throughout the program, we should start from a relatively small value ($a=0.001$). In the piece of code below h represents the $\delta$ increment, kmax is the maximum number of interactions for each eigenvalue, nmax is the number of eigenvalues that we are trying to find, Ein is the minimum energy of de effective potential that we are trying to solve, and Vmax is the maximum energy of the same effective potential.

With the inputs given in the previous step, we can now print the initial values of our program on the screen.

Finally, at this point the program was able to determine the eigenvalues and eigenfunctions associated with the first three states of the hydrogen atom. Initially, the program prints the eigenvalues (Eingen) on the screen according to the Figure 2.

Figure 2
Output of eigenvalues displayed by the code.

Table 1 shows the comparison between our results, extract from the Figure 2, and the well known analytical values for the hydrogen atom, in Rydberg units, given by the formula $E_n=-1/n^2$.

Table 1
Comparison between the energy values, in Rydberg units, found by the Numerov’s method and the analytical ones for $\mathrm{\ell }=0$ and n=1, 2 and 3.

n	Numerov’s Energy	Analytical value
1	$-$0.995	$-$1
2	$-$0.245	$-$0.25
3	$-$0.110	$-$0.111

From now on we will introduce the codes necessary to generate the graphs, as well as calculate the normalization of the wave functions. First, let’s plot the effective potential graph. During the process of setting the code for different potentials, it is important to know the effective potential in order to adjust the parameters accordingly.

The code above plots the effective potential, as we can see in Figure 3

Figure 3
Effective potential of equation (11) for $\mathrm{\ell }=0$.

And the graph for the eigenfunctions with an arbitrary normalization, is built by

Which generates the plot shown in the Figure 4.

Figure 4
First three eigenfunctions of equation (11), with $\mathrm{\ell }=0$, calculated with the Numerov method.

As a last step, let’s include in our code, the calculation of the normalization of the wave functions, given by:

Thus, we were able to obtain our final graph, composed of the first three wave functions for the hydrogen atom (Figure 5).

Figure 5
First three normalized eigenfunctions of equation (11), with $\mathrm{\ell }=0$, calculated with the Numerov method.

4. Morse Potential

The Morse potential is a common model for the interatomic interaction of a diatomic molecule [¹¹,¹²,¹³]. In this section, in order to learn how we can use the Phyton code in other problems, we will see what we must change in the code that was made available in the previous section so that the program will be able to solve equation (2) for the quantum number $\mathrm{\ell }=0$ and the Morse potential with arbitrary parameters given by:

(16) $$V(x)=16{(1-e^{-2x})}^2$$

In this example, we will calculate the first two bounded states, for that, we must adjust the $\text{d}E$ multiplier for each case as:

And, the most important part which is the adjustment of the initial data of our problem. Analyzing the effective potential, we can verify that the value os parameter a must be negative, and its minimum value is zero, so the code must be set as follows:

From now on, the program is able to find the energies of the first two states, which are respectively $7.1380$ and $15.0380$. As well as already determined the wave functions also for the first two states.

Then, all that remains is to get all the aesthetic part of the code right, adjusting the limits of graphics, subtitles and other factors. Thus we first get the effective potential, as shown in Figure 6 and finally arrive at the normalized wave functions shown in Figure 7.

Figure 6
Effective Potential given by equation (16).

Figure 7
Normalized Eigenfunctions of equation (11), with the potential given by equation (16) with $\mathrm{\ell }=0$.

In Table 2, we can compare the eigenvalues found by the Numerov’s method with the analytical values given by $E_n=16\left[\left(n+\frac{1}{2}\right)-\frac{1}{4}\left(n+\frac{1}{2}\right)\right]$.

Table 2
Comparison between the energy values for the Morse potential, in Rydberg units, found by the Numerov’s method and the analytical ones for $\mathrm{\ell }=0$ and n=0 and 1.

n	Numerov’s Energy	Analytical value
0	7.1380	7.0
1	15.0380	15.0

5. Quantum Dot

The development of technology based on quantum dots is quite recent, but it is already showing signs that it is the next great technology, when we talk about a kind of optical devices. In a simple model for a quantum dot composed of two electrons, they can be described with an external harmonic oscillator potential of frequency $\mathrm{\Omega }=2\omega$. Following the steps of reference [¹⁰], we have the effective potential to be introduced into equation (2) for the quantum number $\mathrm{\ell }=0$ is given by:

(17) $$V(x)=\frac{1}{x}+{\omega }^{2}x-\frac{0.25}{x^2}$$

Introducing this potential into our code, let’s now calculate the first five solutions that the program is capable of finding. We can observe that in our code used so far, we only have up to three wave functions, in this example we will show how we included two more solutions in the code. The procedure is very simple, just add the yy4 and yy5 functions along the code and all the other parts that are related to them, for example,

Now, as in the previous example, let’s include the dE multipliers, remembering to include the new wave functions.

And finally, we must include the initial conditions for our problem, as shown below.

Thus, adjusting the rest of the code, we obtain the first five eigenvalues. The effective potential and the eigenfunctions can be seen in Figures 8 and 9.

Figure 8
Effective Potential given by equation (17).

Figure 9
Normalized Eigenfunctions of equation (17), with the potential given by equation (16) with $\mathrm{\ell }=0$.

In Table 3, we can compare the eigenvalues found by the Numerov’s method with the analytical values for the quantum dot given by ${\eta }_{n\mathrm{\ell }}=2(n+\mathrm{\ell }+1)\omega$.

Table 3
Comparison between the energy values for the Quantum dot, in Rydberg units, found by the Numerov’s method and the analytical ones for $\mathrm{\ell }=0$. Results are shown for different values of n and fixed $\omega$ ($\omega =0.01$ Ha).

n	Numerov’s Energy	Analytical value
4	0.1046	0.10
6	0.1403	0.14
8	0.1760	0.18
10	0.2134	0.22
12	0.2507	0.26

6. Conclusion

Analyzing the results arranged in Tables 1, 2 and 3 we can conclude that the method used here is able to reproduce the analytical results within tiny errors.

Therefore, it is evident that the Numerov numerical method is a powerful tool, easy to use, which can help in the development of not only new knowledge in the area of programming, but also in solving the Schrödinger equations outside the usual results found in the examples of modern physics textbooks. We hope that this short highlight to the code, which can be found in full at the links provided in the Introduction, will open doors for students to create their own versions, increasingly improving the versatility of this tool.

Acknowledgments

One of us (FS) was financed in part by the Coordenação de Aperfeiçoamento de Pessoal de Nível Superior – Brazil (CAPES), Finance Code 001. We would like to thank two anonymous referees who contributed to a better presentation of this paper.

References

• ^[1]
  F. Caruso and V. Oguri, Rev. Bras. Ens. Fis. 36, 2310 (2014).
• ^[2]
  B.V. Numerov, Not. Roy. Astro. Soc. 84, 592 (1924).
• ^[3]
  B.V. Numerov, Astronomische Nach. 230, 359 (1927).
• ^[4]
  J.M. Blatt, Jour. Comput. Phys. 1, 382 (1967).
• ^[5]
  A.C. Allison, Jour. Comput. Phys. 6, 378 (1970).
• ^[6]
  J.P. Leroy and R. Wallace, Jour. Phys. Chem. 89, 1928 (1985).
• ^[7]
  A. Bağcı and Z. Güneş, arXiv:2111.11379 (2021).
• ^[8]
  F. Caruso, J. Martins and V. Oguri, Ann. Phys. 347, 130 (2014).
• ^[9]
  F. Caruso, J. Martins, V. Oguri and F. Silveira, Ann. Phys. 377, 518 (2017).
• ^[10]
  F. Caruso, V. Oguri and F. Silveira, Braz. J. Phys. 49, 432 (2019).
• ^[11]
  P.M. Morse and E.C.G. Stueckelberg, Phys. Rev. 33, 932 (1929).
• ^[12]
  P.M. Morse, Phys. Rev. 34, 57 (1929).
• ^[13]
  M.N. Srnec, U. Shiv and D.M. Jeffry, J. Chem. Edu. 94, 669 (2017).

Publication Dates

• Publication in this collection
  10 June 2022
• Date of issue
  2022

History

• Received
  24 Mar 2022
• Reviewed
  02 May 2022
• Accepted
  04 May 2022