In quantum mechanics, the wavefunction of a particle can appear in places where a classical particle would never be found. For instance, a portion of the wavefunction can appear on the other side of a barrier which is impassible to a classical particle. This occurs in quantum mechanics due to the uncertainty in the position of the wavefunction.
The simulation below shows how the probability distribution of a Gaussian wavefunction changes in time according to a Hamiltonian with the potential
where a is the height of the barrier in the centre and b is its width. These two parameters can be adjusted in order to observe how they determine the behaviour of the wavefunction.
The units used here are dimensionless, with ℏ=m=1.
Elise Kenny
2016
barrierWidth = 0.12;
barrierHeight = 2.1;
xmax = 10;
dt = 1.0e-2; % Time step size
N = 512; % Simulation spatial resolution
x = (-N/2:N/2-1)*2*xmax/N; % Position axis
p = (-N/2:N/2-1)*2*pi/xmax; % Momentum axis
% Calculate potential
potential = 0.25*x.^2 + barrierHeight*exp(-(x/barrierWidth).^2);
% Calculate initial wavefunction
wavefunction = sqrt(sqrt(1/(2*pi)))*exp(-0.25*(x+4).^2);
% Calculate the propogators
kprop = exp(-1i*dt*p.^2);
xprop = exp(-1i*dt*potential);
% Setup plot
figure;
[ax, hp, hw] = plotyy(x, potential, x, abs(wavefunction).^2);
axis(ax(1), [-xmax xmax 0 25]); set(ax(1), 'Ytick', 0:5:25);
axis(ax(2), [-xmax xmax 0 1]); set(ax(2), 'Ytick', 0:0.1:1);
legend('Potential', '|Wavefunction|^2');
xlabel('Position');
% Evolve the system and show the animation
while ishandle(hw)
% Evolve the wavefunction
kphi = fftshift(fft(fftshift(wavefunction))).*kprop;
wavefunction = fftshift(ifft(fftshift(kphi))).*xprop;
set(hw, 'ydata', abs(wavefunction).^2);
pause(0.05);
drawnow;
end
import numpy as np
from numpy import fft
import matplotlib.pyplot as plt
import matplotlib.animation as animation
barrierWidth = 0.12;
barrierHeight = 2.1;
animLength = 1000;
xmax = 10;
dt = 1.0e-2; # Time step size
N = 512; # Simulation spatial resolution
x = np.arange(-N/2, N/2)*2.0*xmax/N; # Position axis
p = np.arange(-N/2, N/2)*2.0*np.pi/xmax; # Momentum axis
# Calculate potential
potential = 0.25*x**2 + barrierHeight*np.exp(-(x/barrierWidth)**2.0);
# Calculate initial wave-function
wavefunction = np.sqrt(np.sqrt(1/(2.0*np.pi)))*np.exp(-0.25*(x+4)**2.0);
# Calculate the propagators
kprop = np.exp(-1j*dt*p**2);
xprop = np.exp(-1j*dt*potential);
# Setup plot
fig, ax1 = plt.subplots()
line1, = ax1.plot(x, potential, 'b')
ax1.set_ylim(0.0, 25.0)
ax1.set_xlim(-xmax, xmax)
ax1.set_xlabel("Position")
ax2 = ax1.twinx()
line2, = ax2.plot(x, np.abs(wavefunction)**2.0, 'r')
ax2.set_ylim(0.0, 1.0)
fig.legend((line1, line2), ('Potential', '|Wavefunction|^2'))
def animate(step_number):
global wavefunction
# Evolve the wave-function
kphi = fft.fftshift(fft.fft(fft.fftshift(wavefunction)))*kprop;
wavefunction = fft.fftshift(fft.ifft(fft.fftshift(kphi)))*xprop;
# Update the plot
line2.set_data(x, np.abs(wavefunction)**2.0)
return line2,
ani = animation.FuncAnimation(fig, animate, xrange(1, animLength),
interval=25, blit=True, repeat=False)
#ani.save('output.mp4')
plt.show()
#include <fstream>
#include <vector>
#include <complex>
#include <cmath>
#include <fftw3.h> // For Fourier transforms, link with -lfftw3
using std::complex;
const double barrierWidth = 0.12; // Width of barrier
const double barrierHeight = 2.1; // Height of barrier
const unsigned animLength = 1000; // Number of output frames (files)
const double xmax = 10.0; // Half width of the simulation
const double dt = 1.0e-2; // Time step size
const unsigned N = 512; // Simulation spatial resolution
// Implementation of Matlab's fftshift function
void fftshift(std::vector<complex<double>>& vec) {
std::vector<complex<double>>::iterator mid = vec.begin() + vec.size()/2;
std::swap_ranges(vec.begin(), mid, mid);
}
int main(int argc, char** argv) {
std::vector<complex<double>> kprop(N);
std::vector<complex<double>> xprop(N);
std::vector<complex<double>> wavefunction(N);
std::vector<complex<double>> pwavefunction(N);
std::complex<double> cJ(0.0, 1.0);
// Create a plan for the Fourier transforms
fftw_plan forward, backward;
forward = fftw_plan_dft_1d(N, (double(*)[2]) wavefunction.data(),
(double(*)[2]) pwavefunction.data(), FFTW_FORWARD, FFTW_ESTIMATE);
backward = fftw_plan_dft_1d(N, (double(*)[2]) pwavefunction.data(),
(double(*)[2]) wavefunction.data(), FFTW_BACKWARD, FFTW_ESTIMATE);
std::ofstream fpp("potential.dat");
fpp << "# x\tPotential\n";
for (unsigned i = 0; i < N; ++i) {
double x = 2*i*xmax/N - xmax;
double p = 2*M_PI*i/xmax - N*M_PI/xmax;
double potential = 0.25*x*x + barrierHeight*exp(-pow(x/barrierWidth, 2.0));
// Calculate the propagators
kprop[i] = exp(-cJ*dt*p*p)/pow(N, 0.5);
xprop[i] = exp(-cJ*dt*potential)/pow(N, 0.5);
// Calculate initial wave-function
wavefunction[i] = sqrt(sqrt(1.0/(2.0*M_PI)))*exp(-0.25*pow(x+4.0, 2.0));
// Write the potential to file
fpp << x << '\t' << potential << '\n';
}
fpp.close();
fftshift(kprop);
for (unsigned i = 0; i < animLength; ++i) {
// Evolve the wave-function
fftw_execute(forward);
for (unsigned k = 0; k < N; ++k) pwavefunction[k] *= kprop[k];
fftw_execute(backward);
for (unsigned k = 0; k < N; ++k) wavefunction[k] *= xprop[k];
// Write the wave-function to file
char buf[128];
snprintf(buf, 128, "output_%06d.dat", i);
std::ofstream fp(buf);
fp << "# x\t|Wavefunction|^2\n";
for (unsigned k = 0; k < N; ++k) {
double x = 2*k*xmax/N - xmax;
fp << x << '\t' << pow(abs(wavefunction[k]), 2.0) << '\n';
}
}
// Cleanup the plans
fftw_destroy_plan(forward);
fftw_destroy_plan(backward);
return 0;
}
set terminal gif animate delay 10
set output 'output.gif'
set xrange [-10:10]
set x2range [-10:10]
set yrange [0.0:25.0]
set y2range [0.0:512*512.0]
set xlabel "Position"
do for [i=0:1] {
plot "potential.dat" with lines title "Potential" axes x1y1, \
sprintf("output_%06d.dat", i) with lines \
title sprintf("|Wavefunction|^2 (at %d)", i) axes x2y2
}