Computational Guided Inquiry for PChem (Neshyba, 2014)
Exploring 1st-order perturbation theory
¶OBJECTIVE: Obtain 1st-order solutions to a perturbed particle-in-a-box
SELF-ASSESSMENTS:
- Define the terms zeroth-order eigenfunction, corrected eigenfunction, and ...
Introduction¶
In 1st-order perturbation theory, one constructs new, "corrected" eigenfunctions and eigenenergies based on known solutions to Schroedinger's time-independent equation for a similar, usually simpler system. We call the simpler system the zeroth-order system, $H_o$, and the system of interest is called the perturbed system. They are related according to
$ H = H_o + V(x) \qquad$ (1)
where $V(x)$ is called the perturbation potential. If $H_o$ is a particle in a box, we know its eigenfunctions are simple sinusoids, which we'll denote $\psi_{0,n}$, and its eigenenergies are proportional to the square of the quantum number $n$, and denoted $E_{0,n}$. See http://en.wikipedia.org/wiki/Particle_in_a_box for details. A perturbation to this system can take on lots of different forms. Here's one, for example,
$V(x) = V_o(\dfrac{x}{L})^p \qquad$ (2)
whose form ensures that the perturbation potential begins at zero at the left-hand-side of the box, and increases to maximum value of $V_o$ at the right-hand-side; the parameter $p$ controls how the potential rises (e.g., if $p=2$, the potential rises quadratically). First-order perturbation theory provides a way to calculate approximate eigenfunctions and energies of the perturbed problem (i.e., $H$). These new eigenenergies are given by
$ E_{1,n} = E_{0,n}+<\psi_{0,n} | V | \psi_{0,n}> \qquad$ (3)
where the subscript "1" is meant to indicate this is a first-order perturbation theory result. The new eigenfunctions are given by
$ \psi_n = \psi_{0,n} + \sum{c_k \psi_{0,k}} \qquad$ (4)
where the summation excludes n=k, and where the coefficients are given by
$ c_k = \dfrac{<\psi_{0,n} | V | \psi_{0,k}>}{E_{0,n}-E_{0,k}} \qquad$ (5)
The notation $<\psi | ... | \psi>$ indicates an integral (often referred to as a matrix element), in this case given by
$ <\psi_{0,n} | V | \psi_{0,k}> \quad = \quad \int{\psi_{0,n}^*V\psi_{0,k} dx} \qquad$ (6)
See http://en.wikipedia.org/wiki/Perturbation_theory for details.
Resources
import numpy as np
from numpy import linspace, sqrt, pi, sin, trapz, zeros, size
%matplotlib inline
import matplotlib.pyplot as plt
from matplotlib.pyplot import plot, title, show, xlabel, ylabel, legend
This code defines some constants
hbar = 1
m = 1
This code calculates zero-order wavefunctions
# Define parameters of the system
L = 1.0
x = linspace(0,L,101)
# This inline function returns the nth eigenfunction of a particle-in-a-box
def psi0(n):
A = sqrt(2./L)
kn = n*pi/L
psi = A*sin(kn*x)
return psi
# This inline function returns the nth eigenvalue of a particle-in-a-box
def E0(n):
kn = n*pi/L
E = (hbar*kn)**2/(2*m)
return E
# Testing the inline functions
n = 1
plot(x,psi0(n+1))
title('zeroth-order eigenfunction of POB')
show()
plot(x,psi0(n)**2)
title('zeroth-order probability density of POB')
show()
print "E = ", E0(n)
print "Normalization = ", trapz(psi0(n)**2,x)
Pause for Analysis #1: Sketch these graphs (both $\psi$ and $\psi^2$) in your comp book. Then modify the code to produce eigenstates 2 and 3, and sketch these too.
Next, we create an inline function that generates a perturbation to the hamiltonian.
# This inline function defines a pertubration running from 0 to V0, proportional to x^2
def pert2(V0):
V = (x/L)**2*V0
return V
# Graph it to see what it looks like
V = pert2(.1)
plot(x,V)
ylabel('V(x) (au)')
xlabel('x')
title('Perturbation potential')
Pause for Analysis #2: Sketch this function in your comp book for a few different values of $V_o$ (superimposed on the same graph) so you have a sense of how it affects the perturbation.
The code below uses a trapezoial rule to calculate $<\psi$|V|$\psi>$, the 1st-order perturbation energy correction to the energy of perturbed $\psi$.
# Choose the ground-state zero-order wavefunction
n = 1
# Calculate the perturbation
V = pert2(0.1)
# Get 1st order energy correction
E1n = trapz(psi0(n)**2*V,x)
print "Old energy, correction, and new energy:", E0(n), E1n, E0(n)+E1n
Pause for Analysis #3: Run this cell a few times, for increasing perturbation strength (e.g., $V_o$ = 0.1 to 1). Record any trends you see in your comp book.
Next, to test your skill: create a new inline function, called "pert3", that defines a perturbation running from 0 to $V_o$, but is proportional to x^3. Test it by graphing it, then evaluate the the 1st-order pertubation energy correction.
# This inline function defines a pertubration running from 0 to V0, proportional to x^3
# Graph it to see what it looks like
# Choose the ground-state zero-order wavefunction
# Calculate the perturbation
# Get 1st order energy correction
A little python interlude ...
A little python interlude ... In Python, the "range" command lets you construct a list of numbers. Play with the code below to see how this works.
nlist = range(1,5)
print nlist
Also check out how you can select out some iterations.
for k in nlist:
if k==2:
pass # this does nothing
else:
print k
... and now continuing with our project.
The code below repeats the above energy analysis, but adds on a calculation of the 1st-order corrected wavefunction.
# Choose a zero-order wavefunction
n = 2
# Calculate a perturbation potential, V(x)
V = pert2(0.1)
# Get the 1st order energy correction, <psi|V|psi>
E1n = trapz(psi0(n)**2*V,x)
print "Old energy, correction, and new energy:", E0(n), E1n, E0(n)+E1n
# Generate a list of zeroth-order wavefunctions to sum over
kmax = 8
nlist = range(1,kmax)
cklist = zeros(size(nlist))
print nlist
# Loop to get a 1st-order corrected wavefunction (this is not normalized)
psi_perturbed = psi0(n) # Start off with the zeroth order wavefunction, with coefficient = 1
ik = 0
for k in nlist:
if k==n:
pass
else:
ck = trapz(psi0(n)*V*psi0(k))/(E0(n)-E0(k)) # Calculate the coefficient for the kth 0-order eigenfunction
cklist[ik] = ck
psi_perturbed += ck*psi0(k) # Add that contribution to the corrected eigenfunction, with appropriate coefficient
print k, ck**2 # Print some results
ik += 1
print "cklist =", cklist
# Normalize
norm1 = trapz(psi_perturbed**2,x) # Get <psi|psi> for the corrected eigenfunction
psi_perturbed = psi_perturbed/norm1**.5 # Normalize this |psi>
norm2 = trapz(psi_perturbed**2,x) # Just checking that the new <psi|psi> equals 1
print "Old, new normalization:", norm1, norm2
# Graph the zero-order and corrected eigenfunctions side-by-side
plot(x,psi0(n))
plot(x,psi_perturbed)
title("probability amplitudes")
legend(['Old','New'])
show()
# Also the probabilities
plot(x,psi0(n)**2)
plot(x,psi_perturbed**2)
title("Probability densities")
legend(['Old','New'])
show()
