Solution of the 1d Schr\"odinger Equation for a Symmetric Well

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


I. 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 1-6 . 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 interesting to have certain classes of problems to be solved in a Quantum Mechanical course.
Nowadays, several numerical techniques have 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 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 . Because of that, we suggest a mathematical potential well, which can be expanded in even exponent power series inside the well, becoming an excellent exercise to treat it perturbatively or using the WKB method or other well-known 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 II we present the mathematical potential well. In section III we apply the potential well to the 1d-Schrödinger equation and analyses some particular aspects for both analytical and WKB approximation. In section IV we present the Stationary Perturbation Theory to calculate the energy levels and some wave functions for comparison. In section V we discuss about the Numerov's numerical approach to solve the 1d-Schrödinger equation. Finally, in section VII we draw conclusions and perspectives.

II. 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 where k has m −1 units, r = 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 equals π. As a consequence k = π/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 |z| 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.
In particular, we are interested in those curves with minimum values given by 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 = |x|. However, as we already mentioned, Eq. (1) is spherically symmetric, which means that if we substitute r = |x| by −r = −|x|, the potential remains the same. As a consequence we can drop the absolute value of |x| and write only x to have a 1d version of this mathematical potential FIG. 2. 1d plot of the potential.
This 1d potential is plotted in figure 2. As one can see, in the symmetrical interval kx ∈ [−π, π], it takes a minimum value at V 0 /3 and diverges at ±π. This is the potential we will use to solve the 1d Schrödinger equation. To assess the singularities of this potential we firstly consider the series expansion of the cotangent function around the origin given by where B 2n are the Bernoulli numbers: B 0 = 1, B 1 = − 1 2 , B 2 = 1 6 , B 4 = − 1 30 , B 6 = 1 42 , B 8 = −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 where the singularity is now removed. If we use the value of B 0 = 1 and divide both sides by (kx) 2 with x = 0, we can rewrite this series as This series converges for all the values in the interval kx ∈] − π, π[, as one can check taking the x → 0 limit in Eq. (6). Thus, the 1d potential can be expressed as The first terms of the infinite series are given by 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.

III. 1D SCHÖDINGER EQUATION
The 1d Schrödinger equation is given by 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 ξ = kx such that Eq. (9) can be rewritten as where φ(ξ) ≡ ψ(ξ/k), υ 0 = mV0 2 k 2 , and ε = mE 2 k 2 . We want to solve this equation for ξ with the boundary condition φ(±π) = 0, where the potential diverges. Note that υ 0 is dimensionless. V 0 can assume any positive value in dimensions of 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, then we plug the above identity in Eq. (10) to obtain As we can note, if we choose a Spherical Bessel Function j 0 (ξ) = sin(ξ)/ξ we obtain, In this sense we must impose υ 0 = 1 and 2ε = 1 to obtain a common solution for Spherical Bessel Function of the first kind of the zeroth-order. Consequently we have valid for a wave function φ(ξ) = A sin(ξ)/ξ. We shall demonstrate a posteriori that this is actually the ground state wave function with eigenvalue E = 0.5 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 (±π) → ∞ implies φ(±π) = 0. In the pedagogical exercise of onedimensional square well potential with perfectly rigid walls in a symmetrical interval −a ≤ x ≤ a where V (x) = 0 and V (±a) = ∞, the WKB method gives us an exact solution of energy levels [3][4][5] . Clearly the potential (7) has rigid walls in −π ≤ kx ≤ π 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 where σ(x) = 1 p(x)dx, and p(x) is the momentum given by The linear combination of solutions (15) is also a solution and Due to the symmetry of the well −π ≤ kx ≤ π, it admits even wave functions and odd wave functions As already mentioned, the potential (7) has rigid walls in kx = ±π and the wave function must be zero at these points. For even wave function, φ even (π) = 0, we have while for the odd wave function, φ odd (π) = 0, we have so we assign It follows that then it is straightforward to obtain Now one can write explicitly the potential (7) to , such that the momentum (16) can be written as . and the integral (21) can now be calculated it follows that (23) where we put the label n in E n to indicate that the energy depends on n. Note that for n → ∞, 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 ≈ 0.5244 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 2 k 2 /m. For n = 8 one obtains E 8 = 8.406 2 k 2 /m and 8.333 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 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) where there are some difficulties to calculate it when the potential is a function expanded in a power series.
One must employ the WKB approach calculating the turning points correctly for each energy, which is a tough task once these points depend on the energy, and the condition (21) does not apply anymore. However, 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.

IV. 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 2,4-6 . 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 9 . 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, as we already set V 0 = 2 k 2 m , we write the unperturbed Hamiltonian as where H 0 is the unperturbed Hamiltonian, which is the Hamiltonian for an 1d Harmonic Oscillator with Using the transformation ξ = kx we now write the 1d Schödinger equation as where we have ω = mΩ k 2 ≡ 2 45 and ε = mE 2 k 2 dimensionless constants. Now, it is convenient to define the creation and annihilation operators respectively 6,9 , where we define the operator p = −i d dξ . These operators obey the rules below also we have a |n = √ n |n − 1 and a † |n = √ n + 1 |n + 1 .
(29) The unperturbed Hamiltonian of the one-dimensional harmonic oscillator yields energy values with ε (0) n ≡ mE 0 2 k 2 , where the Stationary Perturbation Theory will add corrections to (30).
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 , where we set the perturbing term as 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) 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 ξ 6 , we will simplify the calculations to second-order of correction to the energy using just only the term ξ 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 The terms of W ξ 4 = 2ξ 4 945 that contribute to secondorder are 1 2ω 2 [a 4 + a † 4 + (4N + 6)a 2 + (4N − 2)a † 2 ], so that it follows n + 2|W ξ 4 |n = 1 84 (4n + 6) (n + 1)(n + 2) n + 4|W ξ 4 |n = 1 84 (n + 1)(n + 2)(n + 3)(n + 4), The values of the matrix elements n| W ξ 4 |n ′ necessary to compute the energy corrections as well as the respective wave functions are listed in Appendix A. In what follows, when we calculate the second-order correction for the ground state energy we find with ε (i) = mE (i) 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 The eigenfunctions of the unperturbed onedimensional harmonic oscillator are where u = mω x ≡ 2 45 1/4 ξ and H n (u) are the Hermite polynomials, Thus, we can calculate the ground state wave function using the equation (42), it follows that we can now put the values of unperturbed functions of (43) into above result and obtain In the same way we can obtain the first excited wave function as fallows 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.

V. NUMERICAL METHOD
For simplicity, we will write the 1D Schrödinger Equation exactly like we did in Eq. (10), where υ 0 ≡ mV0 ℏ 2 k 2 and ε ≡ mE ℏ 2 k 2 . Unlike what we did in the last sections, we will use from now on x as dimensionless variable instead of ξ for convenience.
The Numerov Method 7 can always be applied if the differential equation is in the format Then we want to solve the equation (47) when the potential V (x) is given by we will rewrite it as and if we make g(x) = 2 [ε − υ(x)] then in this case the function s(x) appearing in Eq. (48) is zero. With this in mind, x k = −π + hk, and k = 0, 1, 2, · · · , N . Also assume that at x k we have ψ k = ψ(x k ). Since Eq. (51) is an homogeneous linear differential equation. Therefore, if ϕ 1 (x) is solution to the equation, then is as well ϕ 2 (x) = Aϕ 1 (x), where A is some real constant.
It is worth to mention the special case when A = −1 then ϕ 2 (x) = −ϕ 1 (x). This fact will be important because in the Numerov Method we need two consecutive values for the solution, ψ k and ψ k+1 . So let us assume we start with Ψ 0 = 0 then we can guess Ψ 1 = ±ǫ.
From the Boundary Conditions, at x 0 = −π we must have Ψ 0 = 0, in other words, this point (−π, 0) has to be a value for the solution. Guessing a positive or negative value for Ψ 1 will only drive us to +Ψ(x) or −Ψ(x) then it doesn't matter if we guess a positive or negative variation value in the neighborhood of ψ 0 , since we will get either Ψ or −Ψ but both are solutions. Assume that the amplitudes of some solution ψ when x is discretized in the interval [−π, π] be given by ψ = {ψ 0 , ψ 1 , ψ 2 , · · · , ψ N }. Accordingly, Ψ = Aψ = {Aψ 0 , Aψ 1 , Aψ 2 , · · · , Aψ 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ψ 1 = ǫ, for instance, ǫ = 0.0001, or in other words: A = 10 −4 /ψ 1 . That is reasonable-since Ψ 0 = 0, and because the solution is assumed to be a smooth function-then, we will not expect huge variations from one value to the ones in some neighborhood. If we can calculate the next value Ψ 2 based on Ψ 0 and Ψ 1 , we can repeat the process and find all other values of Ψ k . We can use the Numerov's Method to find Ψ k+1 = f (Ψ k , Ψ k−1 ), which is derived in the Appendix B. This method has order 4 of convergence 7 .
The selection goes as following, for those values of E, starting with Ψ 0 = 0 and Ψ 1 = ǫ which finish with Ψ N ≈ 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 ±π or Ψ(π) = Ψ(−π) = 0.

VI. RESULTS
Unlike what was done in section III 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 with n = 0, 1, 2, · · · (52) where x1 ≡ x1(E) and x2 ≡ x2(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 using Newton-Raphson algorithm to find roots, varying the energy with 0.0001 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 III we could assign an initial value below these values to find the energy which minimizes this new function for a given n. The resulting energies can be visualized in the first column in Table I. 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 3,4 . Now we call attention to the first level of Numerov's algorithm calculation. The value of 0.500 2 k 2 /m is exactly the ground state value obtained analytically in section III. 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 |ψ| 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 III 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  I. 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 I.

VII. 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 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.