Accessibility / Report Error

General recursive solution for one-dimensional quantum potentials: a simple tool for applied physics

Solução recursiva geral para potenciais quânticos unidimensionais: uma ferramenta simples para física aplicada

Sérgio L. Morelhão André V. Perrotta About the authors


A revision of the recursive method proposed by S.A. Shakir [Am. J. Phys. 52, 845 (1984)] to solve bound eigenvalues of the Schrödinger equation is presented. Equations are further simplified and generalized for computing wave functions of any given one-dimensional potential, providing accurate solutions not only for bound states but also for scattering and resonant states, as demonstrated here for a few examples.

one-dimensional quantum potentials; numerical solutions; computational methods; Schrödinger equation

Uma revisão do método recursivo proposto por S.A. Shakir [Am. J. Phys. 52, 845 (1984)] para solucionar autovalores de estados ligados da equação de Schrödinger é apresentado. As equações são ainda mais simplificadas e generalizadas para computar funções de onda em qualquer potencial unidimensional, fornecendo soluções acuradas não somente para estados ligados mas também para problemas de espalhamento e ressonância, como demonstrado aqui em alguns exemplos.

potenciais quânticos unidimensionais; soluções numéricas; métodos computacionais; equação de Schrödinger


General recursive solution for one-dimensional quantum potentials: a simple tool for applied physics

Solução recursiva geral para potenciais quânticos unidimensionais: uma ferramenta simples para física aplicada

Sérgio L. Morelhão11 E-mail:; André V. Perrotta

Instituto de Física, Universidade de São Paulo, São Paulo, SP, Brazil


A revision of the recursive method proposed by S.A. Shakir [Am. J. Phys. 52, 845 (1984)] to solve bound eigenvalues of the Schrödinger equation is presented. Equations are further simplified and generalized for computing wave functions of any given one-dimensional potential, providing accurate solutions not only for bound states but also for scattering and resonant states, as demonstrated here for a few examples.

Keywords: one-dimensional quantum potentials, numerical solutions, computational methods, Schrödinger equation.


Uma revisão do método recursivo proposto por S.A. Shakir [Am. J. Phys. 52, 845 (1984)] para solucionar autovalores de estados ligados da equação de Schrödinger é apresentado. As equações são ainda mais simplificadas e generalizadas para computar funções de onda em qualquer potencial unidimensional, fornecendo soluções acuradas não somente para estados ligados mas também para problemas de espalhamento e ressonância, como demonstrado aqui em alguns exemplos.

Palavras-chave: potenciais quânticos unidimensionais, soluções numéricas, métodos computacionais, equação de Schrödinger.

1. Introduction

Predicting physical properties of quantum systems has basically been treated as an eigenvalue problem. Since from single atoms and molecules up to nanostructured devices, chemical and electronic properties are determined by finding the eigenstates of the system's Hamiltonian. Methods of solving eigenvalue/eigenfunction problems are therefore an inevitable part of the quantum mechanic theory. In real systems, even one-dimensional problems can be quite difficult to solve analytically and often approximative, such as in perturbation theory, and/or numerical methods are required [1-6]. Semiconductor heterostructures are good examples of one-dimensional systems of technological interest where numerical solutions of the Schrödinger equation (SE) are often employed for self-consistent studies of electronic devices [7-11]. Studies of vibrational bound states of molecules in physical chemistry, as well as of other oscillatory systems in atomic and nuclear physics are treated as one-dimensional problems when summarized in the radial SE; they consist another class of important examples in which numerical methods are also often employed [12-15].

Although many procedures and algorithms are available in nowadays for solving the SE numerically - see for instance Refs. [16, 17], and references therein -, they are suitable to address specific problems on distinct fields of applied and theoretical physics. There is still a need for a general treatment that can be summarized into a single routine capable to solve equally well problems of confined modes by arbitrary potential as well as scattering problems in open systems. Moreover, all available numerical methods are essentially mathematical solutions of a differential equation, and hence the level of expertise required for dealing with these methods has kept them out of most physics and chemistry classrooms, even on graduate courses, where the demonstrative examples are still limited to a few cases of quantum potentials for which the analytical solutions are known.

More than 20 years ago, S.A. Shakir [18] presented a recursive method to find bound eigenvalues of the SE. It was a relatively simple method, closely related to those used in the field of optics for calculating the Fresnel amplitudes. In despite of its simplicity, and even of the fact of have been rediscovered a few years latter [19-21], this method has passed unnoted from the applied research fields until very recently [22].

In this work, instead of concerning with mathematical methods to solve the SE for bound and scattering states, we are concerned with alternative methods capable of providing better physical insight about the quantum world. In other words, methods that provide solutions of the SE without having to solve it mathematically. One class of alternative methods, which we call physical methods, is possible in one dimension. They are based on the elemental action of a given stationary force field on the wave function of the particle. Action that occurs at every step of a discretized potential energy function and it is qualitatively always the same, reflecting and transmitting the incident wave as in the Shakir's method [18]. However, here, the original formalism is simplified and generalized into a procedure for computing wave functions of any given potential whatever related to scattering, resonance, tunneling, or bound states, as demonstrated for a few examples. Physical method at their actual stage of development might require higher computational cost than mathematical methods, but offers versatility and exactness in return, as well as the fact that they can be explained and used even by undergraduate students and no-experts in numerical methods.

2. Theory

For a given potential energy function U(x), it is always possible to define the interval X = [x0, xN] such that outside its boundaries, i.e. for x Ï X, either the potential is constant or the amplitude of the wave function vanishes before appreciable variation of the potential. In this interval, the potential has been discretized into N steps, and we seek for stationary states with wave functions in the form

in which

and dj(x) = 1 for x Î [xj,xj+1] and 0 otherwise. The under arrow indicates that this format will be used here for left-hand solutions, as better explained below.

Incremental distances dxj = xj+1-xj, are small enough to guarantee a good numerical solution when taking U(x) = U(xj)dj(x), as shown in Fig. 1, so that

where j = . E and m stand for energy and mass of the particle, respectively.

To calculate the Aj and Bj amplitudes recursively, initial values at the boundaries of the interval X have to be provided. For instance, normalized solutions for incidence from the left-hand side, i.e. incidence from x < x0, must have BN = 0 and A0 = 1. Conservation of the probability current imposes that

at every step of the potential, i.e. at every xj. The prime symbol indicates the first derivative in x. Left-hand solutions are obtained by applying these conditions, Eqs. (4), first at xN where BN = 0 and

and then at xN-1 where


and then at xN-2 where

and AN-2 is analogous to Eq. (7) when replacing N by N-1.

This is a recursive procedure, repeatable since j = N down to j = 1 as summarized by



All amplitudes are proportional to A0, and within the condition assumed when defining the interval X, RN+1 = 0.

Solutions for incident particles from the right-hand side, i.e. incidence from x > xN, are obtained in analogous procedure, but for the sake of simplification in the final format of the recursive equations it is convenient to redefine the amplitudes as follow

and hence the wave function Y(x) of the stationary states will be calculated by using

instead of j(x).

By applying the continuity of

j and 'j, as in Eqs. (4), analogous deduction of that in Eqs. (5)-(8) can be carried out but now starting from x1 where C1 = 0 up to xN+1 where DN+1 = 1 (for normalized solutions). It leads to the recursive equations

of the right-hand solutions where


for j = 1,2,..., N and 0 = 0.

At certain energies E = n, quantum confinements can occur at local minima of the potential in which U(x) < n; given that each minimum is bounded on both sides by non-tunnelable barriers. Around one of such local minima, stationary wave solutions only exist when the waves reflected at both bounding ''walls'' interfere constructively to each other at any instant of time. It means that left-hand and right-hand scattering solutions must coexist in the confinement region, i.e. j(x) = j(x), which takes us back to Eqs. (11). Since Bj = AjRj+1 and Cj+1 = Dj+1

j, as given by Eqs. (9) and (13), all amplitudes cancel each other out in Eq. (11) so that jRj+1 = exp(2ikjdxj). This latter relationship provides a criterion for numerical determination of the allowed modes in any confinement region of the potential since a minimization function such as

can be defined with j running over all values where U(xj) < E. Then, the condition f(E) = 0 is fulfilled only for a set of discrete n values of energy, corresponding to the eigenvalues of the one-dimensional SE at a given local minimum of U(x).

Eigenfunctions Yn, of the energy eigenvalues n, are computed as either j(x) or j(x) solutions in the classical allowed region around the minimum. But both solutions are required to compute the evanescent parts of the eigenfunctions inside the bounding walls and then it is necessary to match the amplitudes of these solutions at some point. If the hth step, for which U(xh) < n, is taken as the matching point, h = h and hence Dh = Ah(1+Rh+1)/(1+h-1). By providing an arbitrary value to Ah, such as Ah = 1, the non-normalized eigenfunction

can be calculated over the entire interval X. Since S = |Yn (xj)|2dxj has to be equal to one for a normalized eigenfunction, is the normalization factor.

3. Examples and discussions

A simple application of this recursive method is to finding transmission coefficients,

across arbitrary potential barriers. It is known as being able to provide exact results even for truncated potentials where the most common approximative methods are used to fail, as compared elsewhere [23]. Besides exactness, another advantage of the present formalism is that it is not just computing the values of T, but the entire wave function instead. Hence, for more complex potential barriers, as those usually find in heterostructures and nanodevices where resonance can take place inside the barriers amplitudes [17, 22, 24, 25], and evanescence times of resonant modes are simultaneously computed.

As an example, consider the potential given in Fig. 2, which is a double barrier with V-shaped well the middle. By scanning the energy values E, Eq. (3), in a chosen range from Emin to Emax in NE steps of width dE so that En = Emin + (n-1)dE and = Emax, a T(E) curve is obtained as the one shown in Fig. 2(a). By storing the coefficients Aj(En) and Bj(En) when calculating the T values, high-resolution plots of the wave functions in the interval X also as a function of energy are already available, i.e.

given that dxj << lj = 2p/|kj|. Hence, amplitudes and attenuations of the waves through the entire barrier can be visualized. For this particular potential, it shows that the mode undergoing resonant-tunneling, at E = 0.406eV, has two maxima in the well as can be seen in Fig. 2(b).

Optionally one can compute only the Rj's and Tj's, Eq. (10), to obtain T = |Tj|2 and either stop the calculation there if only the T(E) curve is desired or compute the Aj's and Bj's, Eqs. (9), later for a chosen set of energies. The options on how the present formalism can be programmed for optimizing computational tasks are detailed in Appendix A.

Although the stored set of solutions, Eq. (18), could be used to study the scattering of wave packets in the barriers, obtaining normalized wave packets would not be as straightforward as possible because the solutions are evenly spaced in energy, not in momentum. Therefore, if one wants to compute wave solutions only once, and also study wave packets composed of these solutions as well as their transmission coefficients, it is better to first define the form of the initial wave packet.

Gaussian wave packets has become a benchmark in computing time evolution of this sort [17], hence let use

as reference for our initial free-particle wave packet in t = 0, centered at x0, and described by a single mode of energy E0 so that k0 = j. It is normalized since ò|Y(x,0)|2dx = 1, and its full width at half maximum (FWHM), W, can be chosen by providing sx = W/@ 0.6 W.

Numerov-type solutions of the time-dependent SE can handled the time evolution of such wave packet, Eq. (19), through a given potential barrier [17], but in the present recursive method it is not possible to start with such expression of Y(x,0) since it is an artificial expression that happens to have the same shape, at t = 0, of the actual wave packet


and sk = sx-1. For sake of normalization, the wave vectors kn must be equally spaced, i.e. kn+1-kn = dk, so that

Since the cn's out of the range kn-k0 = ±3.5sk have negligible amplitude, the NE modes with wave vectors k1, k2,..., are taken in this range to compose the wave packet, whose spectrum of energy is given by En = (kn/j)2. Energy difference between adjacent modes,

increases as n runs toward NE. Hence, the maximum interval of time by which the evolution of the wave packet can be studied is

when the energy range E0 ± DE comprises all En's.

To be more illustrative, let construct an initial guassian wave packet within a previously chosen range of energy, for instance, with En in the narrow range E0 ± 0.058 eV around the energy E0 = 0.406 eV of the resonant mode in Fig. 2. It leads to sk @ jDE/7 = 0.0666 nm-1 and to W = 25.0 nm as the FWHM of |Y(x,0)|2. By composing the wave packet with 101 modes, our limit of time for studying its evolution is Dt < 1654fs, Eq. (23). Figure 3 shows the scattering process of such wave packet by the double barrier (Fig. 2) as computed via Eq. (20). The reflection of the entire wave packet occurs in a time interval no longer than 170fs, Fig. 3(a) through Fig. 3(c), but the resonant mode trapped in the well survives for much longer than that, Fig. 3(d). Its evanescence in time is dictated by an exponential decay of time constant of 184.5fs regarding the probability of finding the particle inside the well after t = 232.3fs, Fig. 3(c).

In case of potential with localized minima not reachable by tunneling from either sides, computing both left-hand and right-hand solutions is necessary. The procedure is very similar to the previous one used to obtain the Y(x,En) set of wave functions evenly spaced in energy where En+1-En = dE. But, most of these unilateral solutions do not exist under confinement. Hence, prior to calculate Aj's, Bj's, and Tj's from Eqs. (9) and (10a), we first compute the Rj's and j's, Eqs. (10b) and (14b), to generate the f(E) plot according to Eq. (15). In Fig. 4(a) there are examples of such plot for a Lennard-Jones potential, a type of potential used to describe diatomic molecules and that has a single minimum. Only after providing an energy eigenvalue, n, all other required coefficients for plotting the eigenfunction Yn(x), Eq. (16), are calculated. For n = 1, 2, and 3, the eigenfunctions are shown in Fig. 4(b) where the one-dimensional variable x stands for the nuclei distance in the molecule. Although reduced mass, equilibrium distance, and ionization energy values are close of those for the H2 molecule, the width of the depression in U(x) is much narrow than in the actual molecule [26]. It means a stronger restoring force and hence larger energy gaps in between adjacent vibrational states, which is used here for illustrative purposes only. In this hypothetical molecule, the non-constancy of the gap values is more evident, as can be seen in Fig. 4(a). It indicates how the exact solution would differ from the harmonic oscillator approximation.

Computational task involved in finding the energy eigenvalues E = n in which f(E) = 0 can be enormous depending on number and accuracy of the desired eigenvalues. There will be no difficult in determining a few eigenvalues with reasonable accuracy, for instance 1 = -3.4416 eV, 2 = -1.8793 eV, 3 = -0.8660 eV, and 4 = -0.2991 eV in Fig. 4(a), are obtained at once in just a few seconds and with an accuracy of dE/2 = 0.0022 eV (NE = 1000); x0 = 0.002 nm, xN = 0.2 nm, and dxj = 0.0005 nm (N = 396). The problem relies in finding many eigenvalues, with high accuracy, occurring at different densities along a broad energy range. But practical needs of finding all possible eigenvalues are unusual, although it can be realized by developing search routines to reduce computational time since the method does not have any intrinsic limitation; as for instance the frequency distortion, or phase-lag, errors of the Numerov-type methods [12, 13].

Potential functions with singularities do not compromise the applicability of the presented method since only solutions with null wave functions at the singularities are those with physical meaning. In other words, when representing the potential by a discretized function, a finite value has to be assigned to the singularity. It compromises mostly the modes with non-null wave functions at the singularity, but these are artificial modes that do not exist in the actual physical system. Occurrence of artificial modes has been observed when solving the potential U(x) = -e2/[|x|+|e|] with e ® 0. In this case, the odd solutions with respect to the singularity provide the well-known eigenvalues of the hydrogen atom.

Another situation that can also be addressed by the recursive method occurs when the potential function has closely spaced minima, as for instance in the potential shown in Fig. 5. In the present formalism, the two minima (or wells) of this potential can be treated separately after establishing distinct intervals in X for calculating the function f(E). For each En Î [Emin,Emax], we first compute the Rj's and j's in the entire interval X that contains the wells, i.e. from x0 to xN. Then, by selecting a position xb somewhere in between the wells as indicated in Fig. 5(b), two f(En) values are obtained: one for xj Î [x0,xb] and another for xj Î [xb+1,xN], and always obeying that U(xj) < En. Both f(E) curves thus obtained are shown in Fig. 5(a). As tunneling across the barrier separating the wells become relevant, the minima of the f(E) curves are coincidental. It means that the energy levels, or eigenvalues n, in each well are independent from each other only if tunneling is negligible. It is more evident when visualizing the wave functions of the confined modes shown in Figs. 5(b), 5(c), and (d). In practice, potential with two minima are found for Cooper pairs (2 spin-coupled electrons) in Josephson-junction of superconducting materials [27].

In semiconductor heterostructures, the validity of this recursive method in its present formalism holds true when the charge carrier's effective mass can be approximated by a constant value through out the structure. However, accounting for effective mass variation across heterojunctions has only been possible in cases of abrupt interfaces and by using appropriated boundary conditions [28, 29].

Computing wave functions recursively is limited to problems treatable in one dimension, such as spherically symmetric potentials, or cases where the potentials in orthogonal directions are independent from each other, as for instance U(x,y) = UX(x)+UY(y), in which the wave functions are obtained by the product of independent one-dimensional solutions. In all other cases of two and three dimensions, it is not possible to extend this recursive method since the number of unknown amplitude coefficients are more than the number of equations written by matching the wave function and its derivative at every steps of the discretized potentials. Hence only mathematical methods for solving 2D and 3D differential equations are available [3, 5, 29, 30].

4. Conclusions

This work has demonstrated that any one-dimensional quantum potential, wherever representing an open or closed system, can be solved by a single method and within good numerical accuracy. It is a simple tool ready to be used in applied physics or just in developing illustrative examples of quantum mechanics. It should also be emphasized that, in the presented formulation, the quantum potentials are solved by a physical treatment where the particle's wave function is molded by the same elemental action of the given force field. This is a completely different approach of the one that is usually taken where, once the potential is known, the problem became the mathematical solution of a differential equation.


During the time dedicated to this work, financial supports were received from CNPq (301617/1995-3 and 304916/2006-4), FAPESP (2002/10387-5), and CAPES-GRICES (140/05).

Recebido em 28/11/2006; Aceito em 8/5/2007

  • 1
  • There are several ways to program the recursive formalism presented in this work, depending on what information is desired. The flowchart in Fig. 6 describes the basic structure of the routines used to elaborate the examples in Figs. 2 to 5. Once the U(x) function and the discretization of X is defined, the user can choose either route or to scan a given energy range in fixed increments of energy, dE, or wavevector, dk, respectively. The latter route is suitable to study time evolution of wave packets by ending through route , i.e. . Eigenvalues {1, 2,...}, of confined modes are obtained from the f(E) curve in route , and used as input to provide the eigenfunctions Yn(x) through route . The straight route are exclusive for users interested only in transmission coefficients across arbitrary potential barriers.

    • [1] D.R. Hartree, The Calculation of Atomic Structures (Pergamon, London, 1957).
    • [2] D. Bohm, Quantum Theory (Dover, New York, 1989).
    • [3] C.S. Lent and D.J. Kirkner, J. Appl. Phys. 67, 6353 (1990).
    • [4] A.R. Lee and T.M. Kalotas, Phys. Scr. 44, 313 (1991).
    • [5] H.W. Jang and J.C. Light, J. Chem. Phys. 102, 3262 (1995).
    • [6] T.E. Simos, Computers Math. Applic. 33, 67 (1997).
    • [7] N.C. Kluksdahl, A.M. Kriman, D.K. Ferry and C. Ringhofer, Phys. Rev. B 39, 7720 (1989).
    • [8] W.R. Frensley, Rev. Mod. Phys. 62, 745 (1990).
    • [9] T.B. Boykin and G. Klimeck, Eur. J. Phys. 25, 503 (2004).
    • [10] N.B. Abdallah and O. Pinaud, J. Comput. Phys. 213, 288 (2006).
    • [11] M. Dehghan, Math. Comput. Simul. 71, 16 (2006).
    • [12] T.E. Simos and G. Tougelides, Computers Chem. 20, 397 (1996).
    • [13] T.E. Simos and P.S. Williams, Computers Chem. 25, 77 (2001).
    • [14] H. Van der Vyver, Appl. Math. Comput. 171, 1025 (2005).
    • [15] S.A. Chin and P. Anisimov, J. Chem. Phys. 124, 054106 (2006).
    • [16] V. Ledoux, M. van Daele and G. vanden Berghe, ACM Trans. Math. Softw. 31, 532 (2005).
    • [17] C.A. Moyer, Am. J. Phys. 72, 351 (2004).
    • [18] S.A. Shakir, Am. J. Phys. 52, 845 (1984).
    • [19] T.M. Kalotas and A.R. Lee, Am. J. Phys. 59, 48 (1991).
    • [20] T.M. Kalotas and A.R. Lee, Am. J. Phys. 59, 1036 (1991).
    • [21] P. Senn, Am. J. Phys. 60, 776 (1992).
    • [22] S.K. Biswas, L.J. Schowalter, Y.J. Jung, A. Vijayaraghavan, P.M. Ajayan and R. Vajtai, Appl. Phys. Lett. 86, 183101 (2005).
    • [23] A. Zhang, Z. Cao, Q. Shen, X. Dou and Y. Chen, J. Phys. A: Math. Gen. 33, 5449 (2000).
    • [24] C.V.-B. Tribuzy, S.M. Landi, M.P. Pires, R. Butendeich and P.L. Souza, Brazilian J. Phys. 32, 269 (2002).
    • [25] Q.S. Zhu, Z.T. Zhong, L.W. Lu and C.F. Li, Appl. Phys. Lett. 65, 2425 (1994).
    • [26] Y.K. Syrkin and M.E. Dyatkina, Structure of Molecules and the Chemical Bond (Interscience Publishers, Nova Iorque, 1950).
    • [27] D.B. Josephson, Rev. Mod. Phys. 46, 251 (1974).
    • [28] G. Bastard, Wave Mechanics Applied to Seminconductor Heterostructures (John Wiley and Sons, Nova Iorque, 1991).
    • [29] J. Wang, E. Polizzi and M. Lundstrom, J. Appl. Phys. 96, 2192 (2004).
    • [30] T.-M. Hwang, W.-W. Lin, W.-C. Wang and W. Wang, J. Comput. Phys. 196, 208 (2004).

    Appendix A

    1 E-mail:

    Publication Dates

    • Publication in this collection
      07 Dec 2007
    • Date of issue


    • Accepted
      08 May 2007
    • Received
      28 Nov 2006
    Sociedade Brasileira de Física Caixa Postal 66328, 05389-970 São Paulo SP - Brazil - São Paulo - SP - Brazil