Time evolution of a Gaussian wavepacket#
Prof. Jay Foley, University of North Carolina Charlotte
To demonstrate the expansion of a Gaussian wavepacket as a superposition of particle-in-a-box energy eigenstates
To demonstrate the time-dependence of a Gaussian wavepacket traveling between two infinitely tall potential walls
To calculate the momentum uncertainty of the Gaussian wavepacket as a function of time
By the end of this workbook, students should be able to
Expand an arbitrary wavefunction in the energy eigenstate basis, and specifically apply this to the expansion of a Gaussian wavepacket in the basis of the particle-in-a-box energy eigenstates
numpyto perform numeric integration
Use the known time dependence of energy eigenstates of the particle-in-a-box to simulate the time evolution of the Gaussian wavepacket travelling between two infinitely tall ptential walls
Use the capabilities of
matplotlibto animate the motion of the Gaussian wavepacket
Compute time-dependent expectation values of the Gaussian wavepacket, and specifically track the momentum uncertainty of this state as a function of time.
We will demonstrate the time-dependence of an electron initially described by a Gaussian wavepacket that is travelling between two infinitely tall potential walls located at \(x=0\) and \(x = L\). The initial state of this wavepacket will be given by
Such a wavefunction describes a particle initially with a mean momentum and position of \(\hbar k_0\) and \(x_0\), respectively, and \(\sigma\) relates to the spread (or position uncertainty) of the particle. In the code that follows, we will have \(x_0 = 200\), \(k_0 = 0.4\), and \(\sigma = 15\) atomic units. We will specify that the location of the left-most infinite potential is at \(L = 500\) atomic units. We will use the fact that we can we can expand this wavefunction in any complete basis, and we will use the basis of the energy eigenfunctions of the particle in a box with \(L = 500\) atomic units.
Hence, we will expand the wavefunction as follows:
where we have used the known time-dependence of each energy eigenstate. The energy eigenvalues have the form
Since we are working in atomic units, \(\hbar\) and \(m\) will have the value of 1.
The coefficients for the above expansion can be determined simply by
where we have set \(L = 500\). We will use the
trapz(fx, x) to perform a trapezoid rule approximation to
Step 0: Import libraries and initialize a plot object for the animation#
import numpy as np import matplotlib.pyplot as plt from matplotlib import animation, rc from IPython.display import HTML # First set up the figure, the axis, and the plot element we want to animate fig, ax = plt.subplots() plt.close() ### parameters for the gaussian wavepacket and the PIB functions x0 = 200 sig = 15 k0 = 0.4 L = 500 ### grid of x values for the functions x = np.linspace(0,500,500) ### parameters for plot ax.set_xlim((0, L)) ax.set_ylim((-0.03, 0.03)) line, = ax.plot(, , lw=2) # initialization function: plot the background of each frame def init(): line.set_data(, ) return (line,)
Step 1: Define helper functions that will compute the PIB basis functions and energy eigenvalues#
It will be helpful to have helper functions that will evaluate the energy eigenfunctions and eigenvalues of the particle in a box given the quantum number \(n\), the length of the box \(L\), and the grid of \(x\) values for the eigenfunction.
def energy_eigenvalue(n, L): """ Helper function to take the quantum number n of the particle in a box and the length of the box L and return the energy eigenvalue in atomic units. This assumes the mass of the particle is 1 in atomic units (i.e. it assumes the particle is an electron). The length should be in atomic units. Arguments --------- n : int the quantum state of the particle in a box L : float the length of the box in atomic units """ m = 1 return n ** 2 * np.pi ** 2 / ( 2 * m * L ** 2) def energy_eigenfunction_phase(n, L, t): """ Helper function to take the quantum number n of the particle in a box and the length of the box L, and the time value and return the time-dependent phase. This assumes the mass of the particle is 1 in atomic units (i.e. it assumes the particle is an electron). The length and time should be in atomic units. Arguments --------- n : int the quantum state of the particle in a box L : float the length of the box in atomic units t : float the time in atomic units """ ci = 0+1j En = energy_eigenvalue(n, L) return np.exp(-ci * En * t) def energy_eigenfunction(n, L, x): """ Helper function to take the quantum number n of the particle in a box, the length of the box L, and x-coordinate value(s) (single value or list) and return the corresponding energy eigenstate value(s) of the ordinary particle in a box. L and x should be in atomic units, and the maximum value of x should be L. Arguments --------- n : int the quantum state of the particle in a box L : float the length of the box in atomic units x : float (or numpy array of floats) the position variable for the energy eigenstate in atomic units """ return np.sqrt(2 / L) * np.sin(n * np.pi * x / L)
Step 2: Define a helper function that will compute values of the Gaussian wavepacket#
We will define a function to take the \(\sigma\), \(x_0\), and \(k_0\) parameters and \(x\) value(s) that define Eq. (1) and return \(\Psi(x,t_0)\).
def gaussian_wavepacket(x0, sig, k0, x): """ Helper function to take the the spread (sig), mean wavenumber (k0), mean position value (x0) and grid of x-values (x) and return the corresponding Gaussian wavepacket values evaluated on that grid Arguments --------- x0 : float mean position of the wavepacket in atomic units sig : float spread of the wavepacket in atomic units k0 : float the mean wavenumber in atomic units of the wavepacket, related to mean momentum by \hbar k0 x : float (or numpy array of floats) the position variable for the energy eigenstate in atomic units """ ci = 0+1j # <== the imaginary unit pre = 1 / (sig * np.sqrt(2 * np.pi)) # <== prefactor gauss = np.exp(-0.5 * ((x - x0) / sig) ** 2) # <== gaussian envelope plane = np.exp(ci * k0 * x) # <== plane-wave component return pre * gauss * plane
Step 3: Set up plot and function parameters and form the Gaussian Wavepacket at time t=t0#
### Compute the Gaussian wavepacket psi_gp = gaussian_wavepacket(x0, sig, k0, x) plt.plot(x, psi_gp)
/opt/hostedtoolcache/Python/3.8.16/x64/lib/python3.8/site-packages/matplotlib/cbook/__init__.py:1335: ComplexWarning: Casting complex values to real discards the imaginary part return np.asarray(x, float)
[<matplotlib.lines.Line2D at 0x7fdaf93c4c40>]
Step 3: Determine the expansion coefficients for the Gaussian Wavepacket in the PIB basis#
### Create an array of quantum numbers for the PIB states... 1 to 100 should be adequate n_array = np.linspace(1,100,100, dtype=int) ### Each quantum number n has a corresponding (complex) coefficient c_n cn_array = np.zeros(len(n_array),dtype=complex) for i in range(0,100): ### multiply psi_n by gaussian wavepacket integrand = np.conj(energy_eigenfunction(n_array[i], L, x) ) * psi_gp ### get coefficient c_n from the integral of psi_n * Psi_gp cn_array[i] = np.trapz(integrand, x)
Step 4: Put it all together and animate!#
N_time = 300 y = np.zeros(len(x),dtype=complex) t = np.zeros(N_time) # animation function. This is called sequentially def animate(i): y = np.zeros(len(x),dtype=complex) for j in range(0,100): ft = energy_eigenfunction_phase(n_array[j], L, i * 10) fx = energy_eigenfunction(n_array[j], L, x) y = y + cn_array[j] * fx * ft t[i] = i*10 line.set_data(x, np.real(y)) return (line,) anim = animation.FuncAnimation(fig, animate, init_func=init, frames=N_time, interval=100, blit=True) # Note: below is the part which makes it work on Colab rc('animation', html='jshtml') anim