Solving the Time-Dependent Schrödinger Equationa

David G. Robertson†
Department of Physics, Otterbein University, Westerville, OH 43081
(Dated: October 10, 2011)

Abstract
Methods for solving the time-dependent Schrödinger equation in one dimension are discussed.
Possible simulation projects include the study of scattering by barriers and wells, the analysis of
time evolution by expansion in energy eigenstates, and tests of time-dependent perturbation theory.

Keywords: quantum mechanics, Schrödinger equation, time dependence, scattering, perturbation

theory, Crank-Nicholson method

a
Work supported by the National Science Foundation under grant CCLI DUE 0618252.
†
Email: drobertson@[Link]

1
 CONTENTS

I. Module Overview 2

II. What You Will Need 2

III. Quantum Mechanics Background 3

IV. Solving the Time-Dependent Schrödinger Equation 5

V. Simulation Projects 12

References 16

Glossary 18

I. MODULE OVERVIEW

This module discusses methods for solving the time-dependent Schrödinger equation in
one space dimension. This is the central problem in quantum mechanics, and courses in
quantum physics typically devote considerable time to developing solutions in analytically
tractable cases (e.g., square wells, the harmonic oscillator). The ability to determine the
wavefunction numerically allows exploration of much wider class of potentials, however, as
well as the study of e.g., scattering with actual wavepackets rather than plane waves. One
can also test the predictions of time-dependent perturbation theory. Such exercises allow
the development of deeper intuition regarding the properties of quantum systems.

II. WHAT YOU WILL NEED

The minimal physics background required includes quantum mechanics at the level of
a typical sophomore-level course on modern physics, specifically the basic properties of
wavefunctions and the Schrödinger equation. Students will ideally have worked through the
free particle in one dimension and be familiar with wave packets, dispersion, and so on. An
upper-level course in quantum mechanics will allow a richer exploration of problems with
the tools developed here, for example the study of time-dependent perturbations.

2
 On the mathematical side, facility with differential and integral calculus is essential, as
is a basic familiarity with differential equations and complex numbers. Students who have
passed through a sophomore-level modern physics course and the standard calculus sequence
through the sophomore year should have the necessary background.
The computing resources needed are actually fairly minimal, as the “size” of the compu-
tational problem is modest. The natural framework for scientific computing is a high-level
language like Fortran, C or C++, and this is an excellent project for students to sharpen
their programming skills. However, the needed calculations can be done using Matlab or
Mathematica. Support for complex arithmetic is essential.
Some facility for generating plots will also be necessary, for example gnuplot or (less
ideally) Excel. Some suggestions in this regard are given with the simulation projects.

III. QUANTUM MECHANICS BACKGROUND

The central problem of quantum mechanics is to solve the Schrödinger equation, which
determines the time evolutions of the wavefunction for a system [1]. For a single non-
relativistic particle in one dimension this takes the form
∂Ψ(x, t) h̄2 ∂ 2 Ψ(x, t)
ih̄ =− + V (x, t)Ψ(x, t), (3.1)
∂t 2m ∂x2
where m is the particle mass and V (x, t) is the potential energy. Given an initial wavefunc-
tion Ψ(x, 0), the Schrödinger equation determines Ψ for all later (and earlier) times.
The physical meaning of Ψ is that it gives the probability density for finding the particle
at different locations, if its position is measured [3]. Specifically, the probability dP to find
the particle in a small range from x to x + dx at time t is given by

dP = |Ψ(x, t)|2 dx. (3.2)

The probability to find the particle in some finite interval, say between x = a and x = b, is
then the sum of the probabilities for each infinitesimal interval between a and b:
Z b
P = |Ψ(x, t)|2 dx. (3.3)
a

A critical requirement for the consistency of the framework is that the total probability to
find the particle somewhere be unity; hence we must require that
Z ∞
|Ψ(x, t)|2 dx = 1. (3.4)
−∞

3
Wavefunctions that satisfy this condition are said to be “normalized.” An important feature
of the time evolution defind by the Schrödinger equation is that if the wavefunction is
normalized at one time, then it will be normalized for all other times as well. This property
is known as “unitarity.”
The physical importance of the normalization requirement cannot be over-emphasized;
it is this condition that puts the “quantum” in quantum mechanics. Wavefunctions that do
not satisfy this fall into two classes: functions for which the normalization integral is finite
but not equal to one, and functions for which the normalization integral is infinite. Given a
function in the first class, we can easily produce a normalized wavefunction. Assume that
we have found a solution of the Schrödinger equation Ψ for which
Z ∞
|Ψ(x, t)|2 dx = A, (3.5)
−∞

with A a finite number. Then the rescaled wavefunction

1
Ψ0 (x, t) = √ Ψ(x, t) (3.6)
A

will be properly normalized, i.e., will satisfy

Z ∞
|Ψ0 (x, t)|2 dx = 1. (3.7)
−∞

Note that because the Schrödinger equation is linear in the wavefunction, the rescaled wave-
function remains a solution.
Functions for which the normalization integral is infinite, on the other hand, cannot be
normalized at all. Such wavefunctions are simply unacceptable for describing real physical
systems, and must be discarded despite being solutions to the Schrödinger equation [4].
The problem of finding solutions to the Schrödinger equation is usually approached using
the technique of separation of variables. This leads to the “time independent” Schrödinger
equation, which determines the energy eigenstates ψ(x):

h̄2 d2 ψ(x)
− + V (x)ψ(x) = Eψ(x), (3.8)
2m dx2

Here E is a constant to be determined by solving the equation [5]. It has the form of an
eigenvalue equation,
Ĥψ = Eψ, (3.9)

4
where Ĥ is the Hamiltonian operator

h̄2 d2
Ĥ = − + V (x) (3.10)
2m dx2

and E is the eigenvalue. It should be emphasized that the TISE determines both ψ and E,
that is, we find that each solution ψ works only with some particular value of E. We can
label the solutions by an index n, so that each function ψn has a corresponding eigenvalue
En . We must also be sure that our wavefunctions are normalizable; this means that ψn
themselves must be normalizable. It is convenient to require that
Z ∞
|ψn (x)|2 dx = 1. (3.11)
−∞

The solution to the Schrödinger equation can then be obtained by expanding the initial
wavefunction in terms of the stationary states,
X
Ψ(x, 0) = cn ψn (x) (3.12)
n

for some coefficients cn . The full solution is then

X
Ψ(x, t) = cn ψn (x)e−iEn t/h̄ (3.13)
n

An alternative to this procedure is to develop a numerical approach to solving the

Schrödinger equation directly, without expanding in stationary states. It is to this task
that we now turn.

IV. SOLVING THE TIME-DEPENDENT SCHRÖDINGER EQUATION

In this section we will develop techniques for solving the full Schrödinger equation nu-
merically. The first step is to introduce a grid of space points, separated by some distance
∆x, on which we will determine the wavefunction. We will also discretize in time, that is,
evaluate the wavefunction only for a discrete set of times separated by some ∆t. Hence we
replace the function Ψ(x, t) with the discrete values Ψni , where i labels the space point and
n labels the time value. Specifically, the grid points are located at

xi = xmin + i∆x (4.1)

5
and the discrete times are given by (taking the initial time to be t = 0)

tn = n∆t. (4.2)

The next step is to approximate the derivatives in the Schrödinger equation as differences.
These should approximate the corresponding derivatives well as long as ∆x and ∆t are “small
enough.” More precisely, the change in Ψ from one grid point to the next or one time step
to the next should be small compared to Ψ itself.
Now, an obvious approximation to the time derivative is

∂Ψ(x, t) Ψn+1 − Ψni

≈ i (4.3)
dt ∆t

which becomes exact as ∆t → 0. To develop a discrete approximation to the second deriva-

tive with respect to x, imagine that we Taylor expand Ψ(x, t) about the point x (I will
suppress the dependence on t here for clarity):

∂Ψ 1 2 ∂ 2 Ψ
Ψ(x + ∆x) = Ψ(x) + ∆x + ∆x + ··· (4.4)
∂x 2 ∂x2

where all derivatives are evaluated at the point x. Eq. (4.4) implies that

∂ 2Ψ
Ψ(x + ∆x) + Ψ(x − ∆x) = 2Ψ(x) + ∆x2 2
+ O(∆x4 ). (4.5)
∂x

If we now drop the higher order terms as unimportant for sufficiently small ∆x, then this
can be interpreted as an approximate formula for the second derivative:

∂ 2Ψ Ψ(x + ∆x) + Ψ(x − ∆x) − 2Ψ(x)

≈ . (4.6)
∂x2 ∆x2

This difference formula is analogous to, e.g.,

∂Ψ Ψ(x + ∆x) − Ψ(x)

≈ (4.7)
∂x ∆x

which gives a discrete approximation to the first derivative. Actually, it is more analogous
to the “symmetric” first derivative:

∂Ψ Ψ(x + ∆x) − Ψ(x − ∆x)

≈ , (4.8)
∂x 2∆x

showing that there is actually quite a bit of freedom in how these approximations are con-
structed. Both eqs. (4.7) and (4.8) approach the exact derivative as ∆x → 0 and so either is

6
a perfectly valid approximation for finite ∆x. The latter is often preferred, however, because
it is more accurate for a given value of ∆x.
Switching to the grid notation where Ψ(xi ± ∆x, t) = Ψni±1 , the discretized version of the
Schrödinger equation thus takes the form

h̄2 Ψni+1 − 2Ψni + Ψni−1

 n+1
Ψi − Ψni
  
ih̄ =− + Vin Ψni (4.9)
∆t 2m ∆x2

Now, imagine that we know the wavefunction at any one time; eq. (4.9) could then be used
to evolve it forward to the next time step. On has only to solve for Ψn+1
i , which will be
given in terms of Ψ at the current time step and the potential function. The new Ψ then
becomes the basis for a further step, and so on. Hence given some initial wavefunction we
can calculate Ψ at any later (or earlier) time. Note that in general the value of Ψn+1 at
some grid point will depend on the value of Ψ at that point and its nearest neighbors – this
is the effect of the space derivatives in the Hamiltonian.
Unfortunately, this approach is totally useless because it is numerically unstable [2]. The
basic reason for this can be seen as follows. Let us adopt a matrix notation and write eq.
(4.9) as
Ψn+1 − Ψn
 
ih̄ = ĤΨn (4.10)
∆t
Here Ψn represents a vector of values (across the space grid) and Ĥ is a (hermitian) matrix
defined by the requirement that it reproduce the right-hand side of eq. (4.9). Solving for
Ψn+1 then gives  
n+1 i∆t
Ψ = 1− Ĥ Ψn (4.11)
h̄
Now, imagine that we expand Ψn in terms of the eigenvectors of Ĥ, that is the vectors ϕα
satisfying:
Ĥϕα = α ϕα (4.12)

where the eigenvalues α are real (since Ĥ is hermitian). We write

X
Ψn = cnα ϕα (4.13)
α

and in the exact solution to the Schrödinger equation we would have

X
Ψn+1 = cnα ϕα e−iα ∆t/h̄ (4.14)
α

7
The approximate solution eq. (4.11) amounts to replacing the exponential by the first-order
expression
iα ∆t
e−iα ∆t/h̄ ≈ 1 − (4.15)
h̄
But this is a complex number with magnitude greater than one:
1/2
2 ∆t2

iα ∆t
1− = 1+ α 2 ≥1 (4.16)
h̄ h̄

Hence the solution will be unstable; all modes grow without bound as we iterate. In fact it
is “unconditionally unstable” – it is unstable no matter how small we make ∆t. It is just
doomed.
To fix this we need a new discrete equation. Fortunately we have a great deal of freedom
in this regard – the main restriction is that it must reduce to the correct Schrödinger equation
in the limits ∆t → 0 and ∆x → 0. For example, we could replace the right hand side of eq.
(4.9) with its value at the future time step, that is,

Ψn+1 − Ψni h̄2 Ψn+1 n+1

i+1 − 2Ψi + Ψn+1
   
i i−1
ih̄ =− 2
+ Vin+1 Ψn+1
i (4.17)
∆t 2m ∆x

or, in matrix notation,

Ψn+1 − Ψn
 
ih̄ = ĤΨn+1 (4.18)
∆t
This clearly has the same continuum limit as the original equation – it is essentially what
we obtain by taking the “backwards” difference for the time derivative. It is also implicit,
in the sense that Ψn+1 appears on both sides of the equation, and we must solve this system
of linear equations to determine it. Still, this could be done (details to follow) so this is
another possible discretization we can use.
But does it help? Well, yes and no. The evolution is now stable; solving for Ψn+1 we find

1
Ψn+1 = i∆t
Ψn (4.19)
1+ h̄
Ĥ

Now the energy eigenstates are multiplied at each time step by a complex number of mag-
nitude !1/2
1 1
iα ∆t
= 2α ∆t2
(4.20)
1+ h̄ 1+ h̄2

which is less than one. However, the “time evolution operator” in eq. (4.19) – the matrix
that multiplies Ψn to give Ψn+1 – is not unitary! (Recall that Ĥ is a hermitian matrix.) This

8
means that the norm of the wavefunction will not be constant; in fact it will decrease mono-
tonically with time. Since conservation of probability is a central requirement of quantum
theory, we should not use this discretization either.
However, a useful approach is now close at hand. What if we instead replace the right
hand side of eq. (4.9) with the average of those terms at timestep n and timestep n + 1:
 n+1
− Ψn

Ψ 1 n+1 n

ih̄ = ĤΨ + ĤΨ (4.21)
∆t 2
This clearly leads to the same continuum limit (∆x → 0 and ∆t → 0), and so is a valid
approximation. It leads to the solution
! 
n+1 1 i∆t
Ψ = i∆t
1− Ĥ Ψn (4.22)
1+ 2h̄
Ĥ 2h̄

Now the evolution operator is clearly unitary, so that probability will be conserved in this
discrete evolution. In addition, the complex number multiplying each eigenmode has exactly
unit modulus, so the evolution should be stable. This version of the discretization is known
as the Crank-Nicholson scheme, and is the one we shall adopt.
Let us write it out in detail so that we can see the computational problem(s) to be solved.
We begin by putting eq. (4.22) in the form
   
i∆t n+1 i∆t
1+ Ĥ Ψ = 1− Ĥ Ψn (4.23)
2h̄ 2h̄

Now recall that Ĥ acts according to the difference formula we derived earlier; that is
h̄2 ∂ 2
Ĥ = − + V (x) (4.24)
2m ∂x2
with
∂ 2Ψ Ψi+1 − 2Ψi + Ψi−1
2
→ (4.25)
∂x ∆x2
With this substitution and some algebra, you can convince yourself that eq. (4.23) can be
written in the form
Ψn+1 n+1 n+1
i+1 + Ψi−1 + Ai Ψi = Bi (4.26)

where
4im∆x2 2m∆x2 n+1
Ai = −2 + − Vi (4.27)
h̄∆t h̄2
4im∆x2 2m∆x2 n
 
n n n
Bi = −Ψi+1 − Ψi−1 + Ψi 2 + + Vi (4.28)
h̄∆t h̄2

9
Eq. (4.26) is a set of linear equations that determine Ψn+1 in terms of Ψn .
To solve these, observe that the system is tri-diagonal, meaning that if we write it as a
matrix equation M Ψn+1 = B:
    
 A1 1 0 0 0 . . .   Ψn+1
1 B
  1
   n+1   
 1 A2 1 0 0 . . .   Ψ2   B2 
  =  (4.29)
  n+1  
 0 1 A3 1 0 . . .   Ψ3   B3 
 
 . . .. .. .. . .   ..   .. 
.. .. . . . . . .
then the only nonzero elements in M are on the central diagonal and the two neighboring
ones. This reflects that the derivatives couple grid points to their nearest neighbors only.
This allows a straightforward solution by row reduction, as follows.
Let us first convert eq. (4.29) to the form
    
n+1
 1 U 1 0 0 0 . . . Ψ
 1   1 R
 0 1 U2 0 0 . . .   Ψn+1
    
 R 
  2  =  2  (4.30)
 0 0 1 U3 0 . . .   Ψn+1   R3 
    
3
. . . . .  .   . 
.. .. .. .. .. . . . .. ..

by means of the following manipulations. Multiply the first row of eq. (4.29) by 1/A1 ≡ U1 ,
and further define R1 ≡ B1 /A1 = B1 U1 . This gives
    
n+1
 1 U1 0 0 0 . . .   Ψ1   R1 
. . .   Ψn+1
    
 1 A2 1 0 0  B 
  2  =  2  (4.31)
 0 1 A3 1 0 . . .   Ψn+1   B3 
    
3
. . . . .  .   . 
.. .. .. .. .. .. .. ..
.
Next, subtract the first row from the second to obtain
    
n+1
 1 U 1 0 0 0 . . . Ψ
 1   R 1

  n+1  
 0 A2 − U1 1 0 0 . . .   Ψ2   B2 − R1 
 
  =  (4.32)
  n+1  
0 1 A3 1 0 . . .   Ψ3   B3 
 
. .. .. .. .. . .   ..   .. 
.. . . . . . . .
Dividing the second row by A2 − U1 then gives
    
n+1
1 U1 0 0 0 . . .   Ψ1   R1 
1 U2 0 0 . . .   Ψn+1
    
0  R 
  2  =  2  (4.33)
0 1 A3 1 0 . . .   Ψn+1   B3 
    
3
. .. .. .. .. . .   ..   .. 
.. . . . . . . .

10
where
1
U2 ≡ (4.34)
A2 − U1
and
A2 − R1
R2 ≡ = (A2 − R1 )U2 (4.35)
A2 − U1
The procedure generalizes, so that eq. (4.30) is obtained with

1 1
U1 = , Ui = (i > 1) (4.36)
D1 Ai − Ui−1

R1 = B1 U1 , Ri = (Bi − Ri−1 )Ui (i > 1) (4.37)

We can now solve equations (4.30) in reverse order, starting from the last row, which reads
simply
Ψn+1
N = RN (4.38)

The next-to-last row reads

Ψn+1 n+1
N −1 + UN −1 ΨN = RN −1 (4.39)

so that
Ψn+1 n+1
N −1 = RN −1 − UN −1 ΨN (4.40)

and so on. In general we have

Ψn+1
i = Ri − Ui Ψn+1
i+1 (4.41)

The following sequence of operations summarizes the above steps (written in C):

U[1] = 1.0/A[1];
for (i=1; i<N; i++)
U[i] = 1.0/(A[i] - U[i-1]);

R[1] = B[1]*U[1];
for (i=1; i<N; i++)
R[i] = (B[i] - R[i-1])*U[i];

psi[N] = R[N];
for (i=N-1; i>0; i--)
psi[i] = R[i] - U[i]*psi[i+1];

11
As presented, this is actually somewhat wasteful of space – not all quantities need to be
stored in separate arrays, for example. For the present problem this inefficiency is rela-
tively harmless, although for more computationally demanding problems it might be worth
eliminating.

V. SIMULATION PROJECTS

The projects suggested here are all based on the above scheme. Once the basic code is
in place and working (Project 1), the others should be very straightforward – they amount
mainly to modifying the potential and/or initial conditions. I recommend that you imple-
ment these calculations in a programming language like C, C++ or Fortran, but Matlab or
Mathematica could be employed if necessary.
A good way of visualizing the solutions will significantly enhance the value of the projects.
Once can simply write data files containing the wavefunction (real and imaginary parts
and/or magnitude, as desired) at different times. A plotting program like gnuplot can then
be used to view the action.
Even better is to produce animations, e.g., of scattering processes. A simple and effective
way to do this with gnuplot in a Unix environment makes use of gnuplot’s facility for
plotting data from the command line. One has merely to make the code print the necessary
gnuplot commands to stdout and then run it, redirecting (or “piping”) the output directly
to gnuplot. For example, consider the following code snippet:

printf("plot ’-’ with lines notitle\n");

for (i=1; i<N; i++)
printf ("%lf %lf\n", x[i], psi[i]);
printf("e\n");

The first printf issues a gnuplot command indicating that data will be entered via the
standard input. (It also indicates that the data points will be connected with lines, and
that the title will be omitted.) If you entered this on the gnuplot command line, it would
subsequently expect you to enter pairs of (x, y) data via the keyboard, separated by spaces
or tabs and followed by carriage returns. To signal the end of the data entry you would
type “e,” at which point your plot would displayed. The subsequent code above prints just

12
these items to the standard output. If this code is part of a loop over time steps, then at
each time step it generates a new plot with the current values. If the code is then run as,
for example,

./[Link] | gnuplot

then the result is an animation of the time evolution.

A detailed example of this is given in the sample code tdse-gnuplot.c.

1. The Free Particle

The first project should be a simulation of a free particle (V = 0) with a gaussian initial
wavefunction. This problem can be solved exactly, so the correctness (or otherwise)
of the code can be checked easily.

For complete details see any text on quantum mechanics [1]. The normalized initial
wavefunction is given by

(x − x0 )2
 
1
Ψ(x, 0) = 1/4 √ exp − exp [ip0 x/h̄] (5.1)
π σ 2σ 2

This represents a particle localized to within ∆x ∼ σ of x0 , and with average momen-

tum p0 . The solution to the Schrödinger equation is then (for V = 0):

(x − (x0 + p0 t/m))2
 
1
Ψ(x, t) = p exp − exp [i(p0 x − Et)/h̄]
π 1/4 σ(1 + ih̄t/mσ 2 ) 2σ 2 (1 + ih̄t/mσ 2 )
(5.2)
where E = p20 /2m.

In building the simulation you will need an array to hold the values of the wavefunction
at each spatial grid point. Say there are N of these. Recall that to update each grid
point we need the values on the neighboring grid points; obviously something will need
to be done for the first and last point, which are missing one neighbor each. Here you
should just specify Ψ = 0 for all time; this should introduce no problems as long as the
wavepacket is well away from the edges. Alternatively, this is equivalent to putting
infinitely high potential “walls” at these locations.

You will also want to adopt a “reasonable” system of units so that the numerical values
of quantities like x, p0 , etc. are neither too large nor too small. One convenient choice

13
 is to use nanometers (1 nm = 10−9 m) for distance, femtoseconds (1 fs = 10−15 s) for
time, and electron volts (1 eV = 1.6 × 10−19 J) for energy. In these units we have

h̄ = 0.6582 eV · fs (5.3)

c = 300 nm/fs (5.4)

A useful combination of these is h̄c = 197.5 eV · nm. The mass of a particle is conve-
niently expressed in terms of its rest energy, e.g., for an electron me c2 = 0.511 MeV
or 511,000 eV. Hence its mass is me = 511, 000/3002 = 5.68 in these units.

A standard test of a program for time evolution is to evolve the solution forward for a
while, then evolve backwards by the same amount. The solution should return to its
initial value. Note that this does not directly test whether the time evolution is being
simulated accurately, only that it is reversible. Compare your result also to the exact
solution and see what sorts of values for ∆x and ∆t are useful. In general we expect
the difference approximations to be valid when the change in Ψ from one grid point
to the next, or one time step to the next, is small compared to Ψ itself.

Observe the spreading of the wavepacket with time. Try other forms for the initial
wavefunction, for example the Lorentzian function

A2
Ψ(x, 0) = eip0 x (5.5)
(x − x0 )2 + γ 2

where A4 = 2γ 3 /π, or or a step function or “tent” function multiplied by exp(ip0 x).

Add code that will calculate the expectation value of the energy at any time, and
verify that this quantity is conserved.

Verify that the evolution is unitary, that is, that the wavefunction remains properly
normaized as it evolves in time.

Experiment with either of the unsuitable approaches to the discretization, and note
what happens.

2. Scattering

An extremely instructive application of the simulation is to study the scattering of a

wavepacket from a potential barrier or well. Modify your code to do this.

14
 Examine the behavior for infinite and finite barriers. The infinite barrier can be simu-
lated by the requirement that Ψ = 0 at the location of the barrier. For a finite barrier,
observe the penetration of the wavefunction into the classically forbidden region. Ac-
cording to the semi-classical approximation, when a wavepacket with average energy
E encounters a barrier of height V > E, it persists in the classically forbidden region
for a time
h̄
T ∼p (5.6)
E(V − E)
Check this result using your simulation. What happens as V → E?

Next, consider scattering off of barriers and wells of finite width. Use your simulation
to calculate the fraction of the initial probability that is transmitted and reflected,
and compare to the textbook results derived using plane waves.

For the case of the finite well, note the formation of a metastable state that it trapped
(for a while) in the well. The initial wavepacket separates into a “prompt” component
that emerges immediately, and a delayed component that leaks out more slowly.

3. Energy Eigenstates

Study the behavior of a wavepacket in a harmonic oscillator potential,

1
V (x) = mω02 x2 (5.7)
2
Begin by setting the initial wavefunction to be a harmonic oscillator stationary state
[1], and confirm that the probability density |Ψ|2 is time-independent for such states.

Next, study linear combinations of a few stationary states. In this case you can
observe the oscillation of the wave packet in the potential well, with a frequency you
can calculate. Compare the calculation to what your simulation shows.

Verify again that the expectation value of the energy is constant, and has the expected
values.

You can also study the evolution of gaussian or other wavepackets in this potential.

4. Periodic Perturbations

Consider again the harmonic oscillator and add to the Hamiltonian a time-dependent
interaction of the form
V (x, t) = F x cos ωt (5.8)

15
 where F is a constant. In classical theory this would correspond to a sinusoidal (in
time) external driving force. In quantum theory it can be thought of as inducing
transitions between the stationary states of the unperturbed harmonic oscillator.

Let’s assume the particle starts out at t = 0 in its ground state (n = 0). The according
to lowest-order perturbation theory it will be found in the first excited state (n = 1)
at time t with probability
2
F2

sin[(ω − ω0 )t/2]
P (t) = (5.9)
2h̄mω0 ω − ω0
p
where ω0 ≡ k/m is the natural frequency of the oscillator [1]. (The probability for
transitions to other states is zero in lowest order.) This is actually an approximation
to the full perturbation theory result, valid when we are “near resonance,” i.e., when

|ω0 − ω|  ω0 + ω (5.10)

The full expression may be found in the standard texts.

Note that the transition probability itself oscillates in time with a period 2π/|ω − ω0 |,
and has its maximum value at t = π/|ω − ω0 |, where

F2
Pmax = (5.11)
2h̄mω0 (ω − ω0 )2

This should be significantly less than one or perturbation theory is invalid.

Modify your harmonic oscillator code from Project 3 to include such a contribution
and study the transitions. Compare the “exact” transition probability (i.e., computed
using your solution of the full Schrödinger equation) to the result of perturbation
theory. Examine the transition probability as a function of time and observe its
periodic behavior.

[1] For additional details on background, see a standard introductory text on quantum mechanics,
e.g., D.J. Griffiths, Introduction to Quantum Mechanics, 2nd ed. (Prentice Hall, 2005); R.L.
Liboff, Introductory Quantum Mechanics, 4th ed. (Addison Wesley, 2003).

16
[2] For example, S.E. Koonin, Computational Physics (Benjamin/Cummings, 1985); T. Pang,
An Introduction to Computational Physics (Cambridge University Press, 2006); L.D. Fosdick,
E.R. Jessup, C.J.C. Schauble, and G. Domik, Introduction to High-Performance Scientific
Computing (MIT Press, 1996).
[3] It also encodes the probabilities for any other measurement, for example of the momentum.
[4] There is an important case where non-normalizable wavefunctions are useful despite their un-
physicality: plane wave solutions for a free particle (V = 0), corresponding to a particle with
a definite momentum. In this case the lack of normalizability is connected to the failure of
such states to properly respect the uncertainty principle: a particle with a definite p would
have ∆p = 0 and hence ∆x = ∞. However, normalizable states can be constructed as linear
combinations of these plane waves, and, indeed, if sufficient care is exercised, the plane waves
themselves can often be used directly to obtain physical results (e.g., transition and reflection
probabilities for potential scattering).
[5] E can be shown to be a real number, assuming V is real.

17
 GLOSSARY

classically forbidden region Region of space in which the total energy is less than the
potential energy; forbidden in classical theory, but quan-
tum wavefunctions can extend into such regions.

Crank-Nicholson method Discretization of the Schrödinger equation that is unitary

and stable.

energy eigenstate Eigenstate of the Hamiltonian operator for a system. A

stationary state.

expectation value Average value of many measurements on identically pre-

pared systems.

Hamiltonian In quantum theory, the operator corresponding to the total

energy of a system.

harmonic oscillator Classical or quantum system in which the potential energy

is proportional to x2 .

normalization Requirement applied to wavefunctions; its physical content

is that the probability of finding the particle somewhere
should be unity.

perturbation theory Approximation scheme based on considering small pertur-

bations of a solvable problem.

resonance Condition when an oscillatory system is driven at close to

its natural frequency.

round-off error Error in floating-point computations on a computer intro-

duced due to the discrete representation of real numbers.

row reduction Matrix technique for solving systems of linear equations.

Schrödinger equation The central equation of quantum theory, which determines

the time development of the wavefunction.

18
separation of variables Technique for separating a partial differential equation into
ordinary differential equations. The basic assumption is a
product form for the solutions.

stationary state Eigenstate of the Hamiltonian operator for a system. An

energy eigenstate.

transition probability The probability that a quantum system initially in some

state will be found in another state at some later time.

unitarity Property of quantum systems that insures the wavefunc-

tion normalization remains constant in time; reflects con-
servation of probability.

wavefunction The entity that describes the quantum state of a system.

Its modulus squared gives the probability density for po-
sition.

wavepacket A superposition of plane waves to produce a wavefunction

that is localized in space.

19