Periodic Pitfalls

I recently hacked up a simulation of one of the simplest quantum dynamics problems – propagating a free particle in time. In the process, I came to appreciate something about the use of Fourier transforms in quantum dynamics that I did not previously understand. You can see how the wave function of the particle spreads out over time, like it is supposed to, but eventually develops a strange wave pattern that should not exist:

Gaussian wave packet at time zero.
Gaussian wave packet at time zero.
Gaussian wave packet at end of simulation.

The free particle is a surprisingly interesting case in quantum mechanics because it has no stationary states; the solutions to the time-independent Schrodinger equation for the free particle are not normalizable, and thus do not correspond to physically realizable states. In practice, we obtain a wave mechanics description of the free particle that matches our classical intuition of a billiard ball flying through space by using a wave packet, an expansion in plane waves:

 \displaystyle \Psi(x) = \int dk X(k) e^{ikx}

Often in quantum mechanics one will see these types of expressions written with k transformed using the de Broglie relation p = \hbar k, but we will ignore those details here.

A particularly useful type of wave packet is the Gaussian wave packet:

\displaystyle \Psi(x) = e^{-\zeta x^2}

Because the momentum distribution corresponding to this wave packet is also a Gaussian, the Gaussian wave packet matches our classical ideas of a particle moving through space with position and momentum that are both known to within some uncertainty.

How do we know that the electrons emitted from some source are described by a Gaussian wave packet, and not any other type of function? Well, we can measure the momentum distribution of the particles when they hit the detector. There  is experimental and theoretical evidence that thermionic electrons obey a Gaussian momentum distribution.  The momentum operator commutes with the free-particle Hamiltonian, and we thus know that the momentum-space wave function of the particle (and thus the position-space wave function as well), was always Gaussian along its trajectory to the detector.

We propagate this Gaussian in time by using the standard rearrangement of the time-dependent Schrodinger equation for a time-independent Hamiltonian:

\displaystyle i \hbar \frac{d \Psi(x,t)}{dt} = \left(\hat{H}\Psi\right)(x,t_0) \Longrightarrow \Psi(x,t) = \left(e^{-\frac{i}{\hbar}\hat{H}t}\Psi\right)(x,t_0)

The simplicity of our Hamiltonian allows easy application of the propagator operator on the right-hand side using the plane-wave expansion of our wave function:

\begin{aligned}\Psi(x,t) &= \left(e^{-\frac{i}{\hbar} \frac{\hat{p}^2}{2m}t}\Psi\right)(x ,t_0) \\ &= \int dk\ X(k,t) e^{i\left(kx -\frac{\hbar k^2 t}{2m}\right)} \end{aligned}

The really cool part is that the momentum-space wave function can be computed easily by sampling values of the position space wave function along our interval and using the standard discrete Fourier transform numerical method to calculate the coefficients of the plane waves. We can see that the momentum space wave function calculated this way matches the exact result, and is constant in time as it should be:



However, having never worked in signal processing, I had never really pondered the nature of the object yielded by the DFT.  While the analytic FT yields a perfect representation of an aperiodic function, the DFT performed over any finite interval yields a representation with a period equal to the length of the interval.  

This means that when we evolve the wave packet in time, what we really do is evolve an infinite number of evenly-spaced wave packets. This comes with the territory of using the inverse Fourier transform to generate our position-space wave function at each time step – the result is intrinsically periodic:

\displaystyle \Psi(x+L,t) = \sum_{k=0}^{N-1} X_k(t) e^{2\pi i k\frac{(x+L)-x_0}{L}} =\sum_{k=0}^{N-1} X_k(t) e^{2\pi i k\frac{x-x_0}{L}}

We don’t get the choice of ignoring these other images. This can cause problems when the wave packet develops a significant amplitude near the edges of the interval,  as the wave packet images outside of our interval begin to “bleed” into it, corrupting our wave function.

After noticing this effect, I made a quick proof-of-concept demonstration. We begin  by initializing a Gaussian wave packet in the center of our interval and propagating it in time:

\displaystyle \Psi(x,t) = e^{-\frac{i}{\hbar} \frac{\hat{p}^2}{2m}t}e^{-\zeta x^2}

If we apply the Fourier transform and inverse Fourier transform in the course of propagating for one step, we can see that if we plot our new wave function on a larger interval, it has become periodic:


Running the simulation out further causes interference between the images to develop:


In practice, this effect is avoid by increasing the size of the box. In some cases, such as when modeling a truly periodic chemical system such as a crystal, the periodicity of the Fourier representation confers a great advantage.

The code used to make all of these plots is on GitHub.