Abstracts
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 onedimensional potential, providing accurate solutions not only for bound states but also for scattering and resonant states, as demonstrated here for a few examples.
onedimensional 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
ARTIGOS GERAIS
General recursive solution for onedimensional 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^{1}1 Email: morelhao@if.usp.br.; André V. Perrotta
Instituto de Física, Universidade de São Paulo, São Paulo, SP, Brazil
ABSTRACT
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 onedimensional potential, providing accurate solutions not only for bound states but also for scattering and resonant states, as demonstrated here for a few examples.
Keywords: onedimensional quantum potentials, numerical solutions, computational methods, Schrödinger equation.
RESUMO
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.
Palavraschave: 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 onedimensional problems can be quite difficult to solve analytically and often approximative, such as in perturbation theory, and/or numerical methods are required [16]. Semiconductor heterostructures are good examples of onedimensional systems of technological interest where numerical solutions of the Schrödinger equation (SE) are often employed for selfconsistent studies of electronic devices [711]. 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 onedimensional problems when summarized in the radial SE; they consist another class of important examples in which numerical methods are also often employed [1215].
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 [1921], 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 noexperts in numerical methods.
2. Theory
For a given potential energy function U(x), it is always possible to define the interval X = [x_{0}, x_{N}] 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 d_{j}(x) = 1 for x Î [x_{j},x_{j}_{+1}] and 0 otherwise. The under arrow indicates that this format will be used here for lefthand solutions, as better explained below.
Incremental distances dx_{j} = x_{j}_{+1}x_{j}, are small enough to guarantee a good numerical solution when taking U(x) = U(x_{j})d_{j}(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 A_{j} and B_{j} amplitudes recursively, initial values at the boundaries of the interval X have to be provided. For instance, normalized solutions for incidence from the lefthand side, i.e. incidence from x < x_{0}, must have B_{N} = 0 and A_{0} = 1. Conservation of the probability current imposes that
at every step of the potential, i.e. at every x_{j}. The prime symbol indicates the first derivative in x. Lefthand solutions are obtained by applying these conditions, Eqs. (4), first at x_{N} where B_{N} = 0 and
and then at x_{N}_{1} where
and
and then at x_{N}_{2} where
and A_{N}_{2} is analogous to Eq. (7) when replacing N by N1.
This is a recursive procedure, repeatable since j = N down to j = 1 as summarized by
where
and
All amplitudes are proportional to A_{0}, and within the condition assumed when defining the interval X, R_{N}_{+1} = 0.
Solutions for incident particles from the righthand side, i.e. incidence from x > x_{N}, 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 x_{1} where C_{1} = 0 up to x_{N}_{+1} where D_{N}_{+1} = 1 (for normalized solutions). It leads to the recursive equationsof the righthand solutions where
and
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 nontunnelable 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 lefthand and righthand scattering solutions must coexist in the confinement region, i.e. _{j}(x) = _{j}(x), which takes us back to Eqs. (11). Since B_{j} = A_{j}R_{j}_{+1} and C_{j}_{+1} = D_{j}_{+1 }
_{j}, as given by Eqs. (9) and (13), all amplitudes cancel each other out in Eq. (11) so that _{j}R_{j}_{+1} = exp(2ik_{j}dx_{j}). 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 ascan be defined with j running over all values where U(x_{j}) < 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 onedimensional SE at a given local minimum of U(x).
Eigenfunctions Y_{n}, 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(x_{h}) < _{n}, is taken as the matching point, _{h} = _{h} and hence D_{h} = A_{h}(1+R_{h}_{+1})/(1+_{h}_{1}). By providing an arbitrary value to A_{h}, such as A_{h} = 1, the nonnormalized eigenfunction
can be calculated over the entire interval X. Since S = Y_{n }(x_{j})^{2}dx_{j} 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 Vshaped well the middle. By scanning the energy values E, Eq. (3), in a chosen range from E_{min} to E_{max} in N_{E} steps of width dE so that E_{n} = E_{min} + (n1)dE and = E_{max}, a T(E) curve is obtained as the one shown in Fig. 2(a). By storing the coefficients A_{j}(E_{n}) and B_{j}(E_{n}) when calculating the T values, highresolution plots of the wave functions in the interval X also as a function of energy are already available, i.e.
given that dx_{j} << l_{j} = 2p/k_{j}. Hence, amplitudes and attenuations of the waves through the entire barrier can be visualized. For this particular potential, it shows that the mode undergoing resonanttunneling, at E = 0.406eV, has two maxima in the well as can be seen in Fig. 2(b).
Optionally one can compute only the R_{j}'s and T_{j}'s, Eq. (10), to obtain T = T_{j}^{2} and either stop the calculation there if only the T(E) curve is desired or compute the A_{j}'s and B_{j}'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}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 freeparticle wave packet in t = 0, centered at x_{0}, and described by a single mode of energy E_{0} so that k_{0} = j. It is normalized since òY(x,0)^{2}dx = 1, and its full width at half maximum (FWHM), W, can be chosen by providing s_{x} = W/@ 0.6 W.
Numerovtype solutions of the timedependent 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
where
and s_{k} = s_{x}^{1}. For sake of normalization, the wave vectors k_{n} must be equally spaced, i.e. k_{n+1}k_{n} = dk, so that
Since the c_{n}'s out of the range k_{n}k_{0} = ±3.5s_{k} have negligible amplitude, the N_{E} modes with wave vectors k_{1}, k_{2},..., are taken in this range to compose the wave packet, whose spectrum of energy is given by E_{n} = (k_{n}/j)^{2}. Energy difference between adjacent modes,
increases as n runs toward N_{E}. Hence, the maximum interval of time by which the evolution of the wave packet can be studied is
when the energy range E_{0 }± DE comprises all E_{n}'s.
To be more illustrative, let construct an initial guassian wave packet within a previously chosen range of energy, for instance, with E_{n} in the narrow range E_{0 }± 0.058 eV around the energy E_{0} = 0.406 eV of the resonant mode in Fig. 2. It leads to s_{k} @ 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 lefthand and righthand 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 E_{n}_{+1}E_{n} = dE. But, most of these unilateral solutions do not exist under confinement. Hence, prior to calculate A_{j}'s, B_{j}'s, and T_{j}'s from Eqs. (9) and (10a), we first compute the R_{j}'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 LennardJones 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 Y_{n}(x), Eq. (16), are calculated. For n = 1, 2, and 3, the eigenfunctions are shown in Fig. 4(b) where the onedimensional 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 H_{2} 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 nonconstancy 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 (N_{E} = 1000); x_{0} = 0.002 nm, x_{N} = 0.2 nm, and dx_{j} = 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 phaselag, errors of the Numerovtype 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 nonnull 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) = e^{2}/[x+e] with e ® 0. In this case, the odd solutions with respect to the singularity provide the wellknown 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 E_{n} Î [E_{min},E_{max}], we first compute the R_{j}'s and _{j}'s in the entire interval X that contains the wells, i.e. from x_{0} to x_{N}. Then, by selecting a position x_{b} somewhere in between the wells as indicated in Fig. 5(b), two f(E_{n}) values are obtained: one for x_{j} Î [x_{0},x_{b}] and another for x_{j} Î [x_{b}_{+1},x_{N}], and always obeying that U(x_{j}) < E_{n}. 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 spincoupled electrons) in Josephsonjunction 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) = U_{X}(x)+U_{Y}(y), in which the wave functions are obtained by the product of independent onedimensional 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 onedimensional 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.
Acknowledgments
During the time dedicated to this work, financial supports were received from CNPq (301617/19953 and 304916/20064), FAPESP (2002/103875), and CAPESGRICES (140/05).
Recebido em 28/11/2006; Aceito em 8/5/2007
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 Y_{n}(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
Publication Dates

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

Accepted
08 May 2007 
Received
28 Nov 2006