Time evolution of a quantum particle in a one-dimensional box (infinite potential well)

1. Demonstration

The video shows a python simulation of the time evolution of the wave function associated with a quantum particle in a one dimensional box (also referred to as the 'infinite square well'). The theory and implementation is described in the sections that follow.

2. Methodology

2.1 Theoretical Background

It is assumed that the reader is already familiar with how the problem of the one dimensional box (or the infinite well) is solved in Quantum Mechanics. If not, see reference [1] and/or [2] at the end of this article. This section presents a brief background on the equations and the theoretical jargon that will be directly used in writing the code.

The infinite well potential is given by,

Solving the time independent Schrӧdinger's equation for this potential, we obtain time independent separable solutions of the form,

The general time dependent solution can be written as a linear combination of the separable solutions along with their time dependent part,

... (1)

For a given the value of 𝚿 at t = 0 (the initial state) the complex coefficients c_n can be determined by exploiting the the property of orthogonality of the separable solutions,

Using these coefficients and the separable states, the general time dependent solution can be generated using equation (1). Now according to the Statistical Interpretation of Quantum Mechanics, the square of the absolute value of 𝚿 will give the probability density at that point. This is what we are going to plot.

2.2 Python implementation

Note: You will need to install the python libraries (Numpy, Scipy and Matplotlib) if you haven't already. A prior experience in these libraries will also be helpful. More info and examples can be found on their respective websites. See references [3], [4] and [5] for links.

Prerequisites:

Basics of numpy
Numerical integration using scipy (scipy.integrate.quad() )
Plotting and animation in matplotlib. See reference [6] for a brief overview on animation.

We start by importing the required modules,

import numpy as np
import scipy.integrate as integrate
from matplotlib import pyplot
import matplotlib.animation as animation

Next, we define the required constants. The values of ħ and m have been taken according to convenience.

L = 20; h_cross = 1; m = 1/2;

We define the initial wave function at time t = 0 as a travelling (towards the right) gaussian wave packet given by,

The values of 𝜎 and k0 control the standard deviation and the mean momentum of the wave packet respectively. Giving appropriate values to these constants, we implement this in a python function,

def psi0(x):
    if(x<=L and x>=0):
        return 0.4*np.exp(-((x-5)**2)/20)*np.exp(complex(0,2*x));
    else:
        return 0;

A has been set equal to 0.4 arbitrarily. If you want, you can also use a value of A such that the packet is normalised.

We now define a function that will return the time independent stationary states (𝜓ₙ) of the infinite potential well,

def psi_n(x,n):
    if(x>=0 and x<=L):
        return np.sqrt(2/L)*np.sin(n*np.pi*x/L);
    else:
        return 0;

Define a numpy array to store the complex coefficients cₙ, and following that, a function to determine these coefficients based on the initial wave function. Scipy's quad integration function doesn't support the complex data type so the real and imaginary parts have been extracted and integrated separately.

#coefficients (first 30 states)
c = np.zeros(30,dtype = np.complex_);

for i in range(30):
    n = i+1;
    I_Real = lambda x: psi_n(x,n)*np.real(psi0(x));
    I_Imag = lambda x: psi_n(x,n)*np.imag(psi0(x));
    c[i] = complex(integrate.quad(I_Real,0,L)[0],integrate.quad(I_Imag,0,L)[0]);

The coefficients of only the first 30 states have been calculated. The higher energy states have less contribution towards the net wave, so it is a fair approximation. For more accurate results you can go up to 50 or 100 states at the expense of performance, you will notice some small scale differences but the large scale net behavior of the wave function will remain the same.

Now we define a function to return the net wave function as a function of position and time (equation (1)),

def PSI(x,t):
    val = 0
    for i in range(30):
        n = i+1;
        val+=c[i]*psi_n(x,n)*np.exp(complex(0,-(n**2)*(np.pi**2)*h_cross*t/(2*m*(L**2))))
    
    return val;

Thats it! The square of the absolute of the value returned by this function will give the probability density at that point and time. The following code calls this function over and over and draws an animated plot.

fig = pyplot.figure()
ax = pyplot.axes(xlim=(0, L), ylim=(0.0, 0.4))


x = np.linspace(0,L,100);
PSI_SQ = np.zeros(100);

#plots the probability density
def plot_probability(frame_number):
    t = frame_number/15;
    for i in range(100):
        PSI_SQ[i] = abs(PSI(x[i],t))**2; #probability density
    p = pyplot.fill_between(x, PSI_SQ,color='orange');
    return p,
        
anim = animation.FuncAnimation(fig, plot_probability,frames=1000, interval=30, blit=True)

animation.FuncAnimation calls the plot_probability() function every 30 milliseconds (interval) each time passing the frame number (iterating from 1 to 1000) as an argument to the function. The plot_probability function then plots the the probability density at time t = frame_number/15. You can use a higher or lower value instead of 15 to slow down or speed up the animation. Detailed info about how FuncAnimation works can be found in the reference [7].

You can try experimenting with the code by changing the parameters (for instance the deviation of the gaussian) and observe their effects. A completely different initial wave function can also be used, for instance, the the second half of the demonstration the initial wave function is shaped like a triangle, describe by the following code,

def psi0(x):
    if(x<=10 and x>=0):
        return 0.05*x;
    elif(x<=20 and x>=10):
        return 1-0.05*x;
    else:
        return 0;

3. Final code

import numpy as np
import scipy.integrate as integrate
from matplotlib import pyplot
import matplotlib.animation as animation

L = 20; h_cross = 1; m = 1/2;


#Initial wave function of the particle (not normalised)
def psi0(x):
    if(x<=L and x>=0):
        return 0.4*np.exp(-((x-5)**2)/20)*np.exp(complex(0,2*x));
    else:
        return 0;
    

#stationary states (time independent part)
def psi_n(x,n):
    if(x>=0 and x<=L):
        return np.sqrt(2/L)*np.sin(n*np.pi*x/L);
    else:
        return 0;

    
#coefficients (first 30 states)
c = np.zeros(30,dtype = np.complex_);

#calculate co-efficients using fourier's trick
for i in range(30):
    n = i+1;
    I_Real = lambda x: psi_n(x,n)*np.real(psi0(x));
    I_Imag = lambda x: psi_n(x,n)*np.imag(psi0(x));
    c[i] = complex(integrate.quad(I_Real,0,L)[0],integrate.quad(I_Imag,0,L)[0]);
    
    
#Wave function at time t    
def PSI(x,t):
    val = 0
    for i in range(30):
        n = i+1;
        val+=c[i]*psi_n(x,n)*np.exp(complex(0,-(n**2)*(np.pi**2)*h_cross*t/(2*m*(L**2))))
    
    return val;



#plot and animate

fig = pyplot.figure()
ax = pyplot.axes(xlim=(0, L), ylim=(0.0, 0.4))


x = np.linspace(0,L,100);
PSI_SQ = np.zeros(100);

#plots the probability density
def plot_probability(frame_number):
    t = frame_number/15;
    for i in range(100):
        PSI_SQ[i] = abs(PSI(x[i],t))**2; #probability density
    p = pyplot.fill_between(x, PSI_SQ,color='orange');
    return p,
        
anim = animation.FuncAnimation(fig, plot_probability,frames=1000, interval=30, blit=True)