Ah, quantum mechanics, inscrutable and stupendous way of doing physics and yet, no-one knows why it works. I have always been amazed by the wonders of quantum world and always looked for the ways of grasping its significance. One really illustrative example of the quantum phenomena is called a finite square well, and here I present my solution of this popular problem solved in the numerical computation package of Python.
Finite square well, what does it mean? Actually it’s quite simple to comprehend – a finite square well is a one-dimensional function V(x) which has a constant value V0 everywhere except where |x| < L, when it drops to zero (Figure). V(x) is called the potential function and it determines behavior of the quantum particle. Don’t be afraid now- even though quantum physics was designed to accommodate subatomic particles, you don’t need to know particle physics to understand this post. Quantum theory works well even if you imagine normal everyday objects like footballs, cars or rabbits instead of protons and electrons.
The idea is that particle is bounded within the region -L<x<L, and does not have enough energy to leave it. This is called a bound state. Classical interpretation would be a tiny ball in the box bouncing back and forth from two sides, forever. What sense does it make? Well, in school grade physics you learned that higher the position of the body, bigger is its potential energy. So, top of the mountain”contains” higher potential than the mountain’s bottom. If there was car between two mountains it would be forever stuck in the valley between if it’s engine cannot make more (kinetic) energy than it is potential at the top of the mountains.
Here comes the Schrödinger
How to describe behavior of such particle? In many textbooks, online classes, websites, etc. you’ll find that you can not treat quantum particle as an ordinary physical body (even though it can look like an ordinary body, have mass like ordinary body and so on), simply because you can not know for sure its performances like position, speed or momentum. Instead, you attach some transcendental function which is space dependent and use it to find a probability of particle’s position, momentum etc at the given time. To get nice and clean expression for the -function, we introduce Her Majesty – the Schrödinger equation:
Note: this is time-independent form of the Schrödinger equation, and it will give us a time-independent solutions of . But don’t worry, time-independent solutions are very insightful and – now comes the best part – complete solution is just a linear combination of time-independent solutions!
To make long story short – to describe the behavior of the particle, just give me the Schrödinger equation and the potential function and I will tell you all you want to know about the particle.
Do it the Python way
As you see, problem includes differential equation and it’s solution may be tricky if you are a fan of “pen and paper” approach. But, there are people out there who made computers do all this dirty work, and some of these tools are available in our beloved Python. The one I will use is a routine called odeint() (which stands for ordinary differential equation-something) and is included in standard scipy.integrate package. The odeint() works in a two-state-space representation of : state one is function the way we want it and state two is a first derivative of .
The thing is – this way we can get infinitely many solutions. But not all of them are meaningful. We need to find only those solutions where does not diverge at the edges of x. If it diverges, it disobeys the condition of square-integrability: and we deal with the physical falseness.
So here’s what we do: we run value of particle’s energy E from 0 to some value, say Vo, and calculate for each E. Value of (at the far outside of the finite well – where b>>L), is put into the separate list whose roots (places where it crosses zero) are to be found. Only those energies E for which is zero are taken in consideration and represent valid state of the particle. Here we see first of many weird things happening in the quantum world – particle can not have any value for kinetic energy, only some of discrete numbers.
Equation solving code
As said earlier, to solve the differential equation we use scipy.integrate.odeint() routine. To properly work with it, we need to rewrite the Schrödinger equation in state-space representation, where first state is and second state . Since the odeint() works stepwise, it will take states from the previous step in each computation. The state-space function will finally look like this:
def SE(psi, x): """ Returns derivatives for the 1D schrodinger eq. Requires global value E to be set somewhere. State0 is first derivative of the wave function psi, and state1 is its second derivative. """ state0 = psi state1 = 2.0*(V(x) - E)*psi return array([state0, state1])
(Note: I’ve put m so that , just for simplicity). This function prepares state-space only. Following function does the real job, calls the state-space function and calculates the wave function for a given energy E. It also returns value of at the location x=b to test divergence.
def Wave_function(energy): """ Calculates wave function psi for the given value of energy E and returns value at point b """ global psi global E E = energy psi = odeint(SE, psi0, x) return psi[-1,0]
Of course, for the first step there is no previous step – that’s why we introduce variable , which holds initial conditions: and
Finding allowed states
Once we have the values of for a range of energies we need to distinguish those states where diverges. This is done using the Brent method to find a zero of the function f on the sign changing interval [a , b]. The method is available as a Python’s routine scipy.optimize.brentq(). For an array of energies e we will check the sign of the . When the sign changes, it means that crosses zero point and now we can call brentq() to find the exact zero-point. The function is embedded as follows:
def find_all_zeroes(x,y): """ Gives all zeroes in y = f(x) """ all_zeroes =  s = sign(y) for i in range(len(y)-1): if s[i]+s[i+1] == 0: zero = brentq(Wave_function, x[i], x[i+1]) all_zeroes.append(zero) return all_zeroes
The function will return an array all_zeroes which contains all zero points of the . How nicely the function works, you can see on this image:
In this example, I put V0 to 20, L to 1, and m so that . We obtain only 4 four allowed bound states with energies 0.92, 3.65, 8.09, and 14.
Now there’s a little trick. You remember that fantastic little function which solves differential equations, odeint(), needs some initial conditions to work, right?. These are the values of at x=0. What we get from it is the function only right from 0! What happens on the left side? It depends on how we set the initial conditions. For if we chose and we get something called symmetrical case, so left side is just mirrored right side. But if we reverse initial conditions so that and , we get antisymmetric case. Why is that so? I’ll tell you just that symmetrical states correspond to the cosine based functions, while antisymmetric correspond to the sine based functions. More about it will come a bit later.
Let’s see how the states of look like. For the symmetrical case (, ):
And for the antisymmetrical case (, ):
This is nice and correct way to solve the problem of the finite square well. But, I don’t like this boring pictures, so I decided to cheat a little bit :). You see, the odeint() requires the initial conditions at x = 0, assuming x goes from 0 to b. But what if I give it values -b < x< b? Normally, we expect to be o. But if we put initial conditions to be zero, we won’t get any non-zero result. That’s how the odeint() works, I can not prevent it. So here is where I cheat: I put x to be from -b to b, but initial conditions are not zero. Instead, they are some arbitrarily small value. Of course, odeint() will mess up values, but, interestingly, the shape of the function will stay correct! Additionally, it will calculate both symmetrical and antisymmetrical states. Great!
So, let’s do it again. Energies of allowed states are already calculated, it’s the first graph you see in this article. Bound states are here, presented in symmetrical and antisymmetrical case:
As said earlier, values are messed up, since the integral of squared is not 1. But it doesn’t matter cause is normalized anyway later. Few things about these states.
- Particle exists only in discrete states and is associated with the corresponding energy in each state. The states are orthonormal. One state can not be represented as a linear combination of other states. So, the particle can not be in the state which is not representable as a combination of given states.
- Probability of finding particle in space is integral of . We see that in the first state (ground state) most likely to find the particle is in the center of the well. But if you give the particle a bit of energy and it jumps into the second state (excited state) you will find it most likely between center and edges.
- What happens with the particle behind the walls of the well? We see that is not zero left from -L and right from L. It means that there is a non-zero probability of finding a particle outside the well! Remember the analogy with the car between the mountains? It would mean that car has somehow ended behind the mountain although it can not climb it. The phenomenon is called quantum tunnelling and it may seem crazy in case of a car and a mountain, but it is quite common and normal thing in subatomic particles.
Is all this right?
Did I do something wrong? Who guarantees that my functions odeint() and brentq() are used properly? Are these solutions of Schrödinger equation correct? To check it, I run a numerical computation of analytical model, provided in “Introduction to Quantum Mechanics” by D. Griffiths, page 62. In that model, we need to find solution of the equation where and . This will give us energies of allowed symmetrical states. For antisymmetrical states we should use . Python method that does the job is presented here:
def find_analytic_energies(en): """ Calculates Energy values for the finite square well using analytical model (Griffiths, Introduction to Quantum Mechanics, page 62.) """ z = sqrt(2*en) z0 = sqrt(2*Vo) z_zeroes =  f_sym = lambda z: tan(z)-sqrt((z0/z)**2-1) # Formula 2.138, symmetrical case f_asym = lambda z: -1/tan(z)-sqrt((z0/z)**2-1) # Formula 2.138, antisymmetrical case # first find the zeroes for the symmetrical case s = sign(f_sym(z)) for i in range(len(s)-1): # find zeroes of this crazy function if s[i]+s[i+1] == 0: zero = brentq(f_sym, z[i], z[i+1]) z_zeroes.append(zero) print "Energies from the analytical model are: " print "Symmetrical case)" for i in range(0, len(z_zeroes),2): # discard z=(2n-1)pi/2 solutions cause that's where tan(z) is discontinuous print "%.4f" %(z_zeroes[i]**2/2) # Now for the asymmetrical z_zeroes =  s = sign(f_asym(z)) for i in range(len(s)-1): # find zeroes of this crazy function if s[i]+s[i+1] == 0: zero = brentq(f_asym, z[i], z[i+1]) z_zeroes.append(zero) print "(Antisymmetrical case)" for i in range(0, len(z_zeroes),2): # discard z=npi solutions cause that's where cot(z) is discontinuous print "%.4f" %(z_zeroes[i]**2/2)
As a result, I get output in Python console:
Energies from the analyitical model are: (Symmetrical case) 0.9179 8.0922 19.9726 (Antisymmetrical case) 3.6462 14.0022
and that corresponds completely to programs computation. We can be sure that program works for as long as E is smaller than , and is not unreasonablly large. Finally, the source code for the whole program is here:
# -*- coding: utf-8 -*- """ Created on Mon Sep 03 21:02:59 2014 @author: Pero 1D Schrödinger Equation in a finite square well. Program calculates bound states and energies for a finite potential well. It will find eigenvalues in a given range of energies and plot wave function for each state. For a given energy vector e, program will calculate 1D wave function using the Schrödinger equation in a finite square well defined by the potential V(x). If the wave function diverges on x-axis, the energy e represents an unstable state and will be discarded. If the wave function converges on x-axis, energy e is taken as an eigenvalue of the Hamiltonian (i.e. it is alowed energy and wave function represents allowed state). Program uses differential equation solver &amp;amp;quot;odeint&amp;amp;quot; to calculate Sch. equation and optimization tool &amp;amp;quot;brentq&amp;amp;quot; to find the root of the function. Both tools are included in the Scipy module. The following functions are provided: - V(x) is a potential function of the square well. For a given x it returns the value of the potential - SE(psi, x) creates the state vector for a Schrödinger differential equation. Arguments are: psi - previous state of the wave function x - the x-axis - Wave_function(energy) calculates wave function using SE and &amp;amp;quot;odeint&amp;amp;quot;. It returns the wave-function at the value b far outside of the square well, so we can estimate convergence of the wave function. - find_all_zeroes(x,y) finds the x values where y(x) = 0 using &amp;amp;quot;brentq&amp;amp;quot; tool. Values of m and L are taken so that h-bar^2/m*L^2 is 1. v2 adds feature of computational solution of analytical model from the usual textbooks. As a result, energies computed by the program are printed and compared with those gained by the previous program. """ from pylab import * from scipy.integrate import odeint from scipy.optimize import brentq def V(x): """ Potential function in the finite square well. Width is L and value is global variable Vo """ L = 1 if abs(x) > L: return 0 else: return Vo def SE(psi, x): """ Returns derivatives for the 1D schrodinger eq. Requires global value E to be set somewhere. State0 is first derivative of the wave function psi, and state1 is its second derivative. """ state0 = psi state1 = 2.0*(V(x) - E)*psi return array([state0, state1]) def Wave_function(energy): """ Calculates wave function psi for the given value of energy E and returns value at point b """ global psi global E E = energy psi = odeint(SE, psi0, x) return psi[-1,0] def find_all_zeroes(x,y): """ Gives all zeroes in y = Psi(x) """ all_zeroes =  s = sign(y) for i in range(len(y)-1): if s[i]+s[i+1] == 0: zero = brentq(Wave_function, x[i], x[i+1]) all_zeroes.append(zero) return all_zeroes def find_analytic_energies(en): """ Calculates Energy values for the finite square well using analytical model (Griffiths, Introduction to Quantum Mechanics, 1st edition, page 62.) """ z = sqrt(2*en) z0 = sqrt(2*Vo) z_zeroes =  f_sym = lambda z: tan(z)-sqrt((z0/z)**2-1) # Formula 2.138, symmetrical case f_asym = lambda z: -1/tan(z)-sqrt((z0/z)**2-1) # Formula 2.138, antisymmetrical case # first find the zeroes for the symmetrical case s = sign(f_sym(z)) for i in range(len(s)-1): # find zeroes of this crazy function if s[i]+s[i+1] == 0: zero = brentq(f_sym, z[i], z[i+1]) z_zeroes.append(zero) print "Energies from the analyitical model are: " print "Symmetrical case)" for i in range(0, len(z_zeroes),2): # discard z=(2n-1)pi/2 solutions cause that's where tan(z) is discontinous print "%.4f" %(z_zeroes[i]**2/2) # Now for the asymmetrical z_zeroes =  s = sign(f_asym(z)) for i in range(len(s)-1): # find zeroes of this crazy function if s[i]+s[i+1] == 0: zero = brentq(f_asym, z[i], z[i+1]) z_zeroes.append(zero) print "(Antisymmetrical case)" for i in range(0, len(z_zeroes),2): # discard z=npi solutions cause that's where ctg(z) is discontinous print "%.4f" %(z_zeroes[i]**2/2) N = 1000 # number of points to take psi = np.zeros([N,2]) # Wave function values and its derivative (psi and psi') psi0 = array([0,1]) # Wave function initial states Vo = 20 E = 0.0 # global variable Energy needed for Sch.Eq, changed in function "Wave function" b = 2 # point outside of well where we need to check if the function diverges x = linspace(-b, b, N) # x-axis def main(): # main program en = linspace(0, Vo, 100) # vector of energies where we look for the stable states psi_b =  # vector of wave function at x = b for all of the energies in en for e1 in en: psi_b.append(Wave_function(e1)) # for each energy e1 find the the psi(x) at x = b E_zeroes = find_all_zeroes(en, psi_b) # now find the energies where psi(b) = 0 # Print energies for the bound states print "Energies for the bound states are: " for E in E_zeroes: print "%.2f" %E # Print energies of each bound state from the analytical model find_analytic_energies(en) # Plot wave function values at b vs energy vector figure() plot(en/Vo,psi_b) title('Values of the $\Psi(b)$ vs. Energy') xlabel('Energy, $E/V_0$') ylabel('$\Psi(x = b)$', rotation='horizontal') for E in E_zeroes: plot(E/Vo, , 'go') annotate("E = %.2f"%E, xy = (E/Vo, 0), xytext=(E/Vo, 30)) grid() # Plot the wavefunctions for first 4 eigenstates figure(2) for E in E_zeroes[0:4]: Wave_function(E) plot(x, psi[:,0], label="E = %.2f"%E) legend(loc="upper right") title('Wave function') xlabel('x, $x/L$') ylabel('$\Psi(x)$', rotation='horizontal', fontsize = 15) grid() if __name__ == "__main__": main()