![\"Open](\"https://colab.research.google.com/assets/colab-badge.svg\"){kind=icon}
"
]
},
{
"cell_type": "markdown",
"metadata": {},
"source": [
"# Time evolution of a Gaussian wavepacket\n",
"[Prof. Jay Foley, University of North Carolina Charlotte](https://foleylab.github.io/)\n",
"\n",
"#### Objectives\n",
"- To demonstrate the expansion of a Gaussian wavepacket as a superposition of particle-in-a-box energy eigenstates\n",
"- To demonstrate the time-dependence of a Gaussian wavepacket traveling between two infinitely tall potential walls\n",
"- To calculate the momentum uncertainty of the Gaussian wavepacket as a function of time\n",
"\n",
"#### Learning Outcomes\n",
"By the end of this workbook, students should be able to\n",
"- 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\n",
"- Utilize the `trapz` function of `numpy` to perform numeric integration\n",
"- 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\n",
"- Use the capabilities of `matplotlib` to animate the motion of the Gaussian wavepacket\n",
"- Compute time-dependent expectation values of the Gaussian wavepacket, and specifically track the momentum uncertainty of this state as a function of time."
]
},
{
"cell_type": "markdown",
"metadata": {
"colab_type": "text",
"id": "DqCDR6mqeCqx"
},
"source": [
"#### Summary\n",
"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\n",
"\\begin{equation}\n",
"\\Psi(x,t_0) = \\frac{1}{\\sigma \\sqrt(2\\pi)} {\\rm exp}\\left(\\frac{-1}{2}\\left(\\frac{x-x_0}{\\sigma}\\right)^2\\right) {\\rm exp}(i k_0 x). \\tag{1}\n",
"\\end{equation}\n",
"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)\n",
"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. \n",
"\n",
"Hence, we will expand the wavefunction as follows:\n",
"\\begin{equation}\n",
"\\Psi(x,t) = \\sum_{n=1}^N c_n \\sqrt{\\frac{2}{L}} {\\rm sin}\\left( \\frac{n \\pi x}{L} \\right)\n",
"{\\rm exp}\\left(-\\frac{i E_n t}{\\hbar}\\right), \\tag{2}\n",
"\\end{equation}\n",
"where we have used the known time-dependence of each energy eigenstate. The energy eigenvalues have the form\n",
"\\begin{equation}\n",
"E_n = \\frac{n^2 \\pi^2 \\hbar^2}{2 m L^2}. \\tag{3}\n",
"\\end{equation}\n",
"Since we are working in atomic units, $\\hbar$ and $m$ will have the value of 1.\n",
"\n",
"The coefficients for the above expansion can be determined simply by\n",
"\\begin{equation}\n",
"c_n = \\sqrt{\\frac{2}{500}} \\int_0^{500} {\\rm sin}\\left(\\frac{n \\pi x}{500}\\right) \\cdot \\Psi(x,t_0) \\: dx \\tag{4}\n",
"\\end{equation}\n",
"where we have set $L = 500$. We will use the `numpy` function `trapz(fx, x)` to perform a trapezoid rule approximation to \n",
"this integral."
]
},
{
"cell_type": "markdown",
"metadata": {},
"source": [
"### Step 0: Import libraries and initialize a plot object for the animation"
]
},
{
"cell_type": "code",
"execution_count": 1,
"metadata": {},
"outputs": [],
"source": [
"import numpy as np\n",
"import matplotlib.pyplot as plt\n",
"\n",
"from matplotlib import animation, rc\n",
"from IPython.display import HTML\n",
"\n",
"\n",
"# First set up the figure, the axis, and the plot element we want to animate\n",
"fig, ax = plt.subplots()\n",
"plt.close()\n",
"\n",
"### parameters for the gaussian wavepacket and the PIB functions\n",
"x0 = 200\n",
"sig = 15\n",
"k0 = 0.4\n",
"L = 500\n",
"\n",
"### grid of x values for the functions\n",
"x = np.linspace(0,500,500)\n",
"\n",
"### parameters for plot\n",
"ax.set_xlim((0, L))\n",
"ax.set_ylim((-0.03, 0.03))\n",
"\n",
"line, = ax.plot([], [], lw=2)\n",
"\n",
"# initialization function: plot the background of each frame\n",
"def init():\n",
" line.set_data([], [])\n",
" return (line,)"
]
},
{
"cell_type": "markdown",
"metadata": {
"colab_type": "text",
"id": "xW874_pZD6eg"
},
"source": [
"#### Step 1: Define helper functions that will compute the PIB basis functions and energy eigenvalues\n",
"It will be helpful to have helper functions that will evaluate the energy eigenfunctions and eigenvalues of the \n",
"particle in a box given the quantum number $n$, the length of the box $L$, and the grid of $x$ values for the eigenfunction."
]
},
{
"cell_type": "code",
"execution_count": 2,
"metadata": {},
"outputs": [],
"source": [
"def energy_eigenvalue(n, L):\n",
" \"\"\" Helper function to take the quantum number n of the particle in a box and the length \n",
" of the box L and return the energy eigenvalue in atomic units. This assumes the mass of the particle\n",
" is 1 in atomic units (i.e. it assumes the particle is an electron).\n",
" The length should be in atomic units.\n",
" \n",
" Arguments\n",
" ---------\n",
" n : int\n",
" the quantum state of the particle in a box\n",
" \n",
" L : float\n",
" the length of the box in atomic units\n",
" \"\"\"\n",
" m = 1 \n",
" return n ** 2 * np.pi ** 2 / ( 2 * m * L ** 2)\n",
"\n",
"def energy_eigenfunction_phase(n, L, t):\n",
" \"\"\" Helper function to take the quantum number n of the particle in a box and the length \n",
" of the box L, and the time value and return the time-dependent phase. This assumes the mass of the particle\n",
" is 1 in atomic units (i.e. it assumes the particle is an electron).\n",
" The length and time should be in atomic units.\n",
" \n",
" Arguments\n",
" ---------\n",
" n : int\n",
" the quantum state of the particle in a box\n",
" \n",
" L : float\n",
" the length of the box in atomic units\n",
" \n",
" t : float\n",
" the time in atomic units\n",
" \"\"\"\n",
" ci = 0+1j\n",
" En = energy_eigenvalue(n, L)\n",
" return np.exp(-ci * En * t)\n",
" \n",
"\n",
"def energy_eigenfunction(n, L, x):\n",
" \"\"\" Helper function to take the quantum number n of the particle in a box, the length \n",
" of the box L, and x-coordinate value(s) (single value or list) and return the corresponding\n",
" energy eigenstate value(s) of the ordinary particle in a box. L and x should be in atomic units,\n",
" and the maximum value of x should be L.\n",
" \n",
" Arguments\n",
" ---------\n",
" n : int\n",
" the quantum state of the particle in a box \n",
" L : float\n",
" the length of the box in atomic units\n",
" x : float (or numpy array of floats)\n",
" the position variable for the energy eigenstate in atomic units\n",
" \"\"\"\n",
" return np.sqrt(2 / L) * np.sin(n * np.pi * x / L)"
]
},
{
"cell_type": "markdown",
"metadata": {},
"source": [
"#### Step 2: Define a helper function that will compute values of the Gaussian wavepacket \n",
"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)$."
]
},
{
"cell_type": "code",
"execution_count": 3,
"metadata": {},
"outputs": [],
"source": [
"def gaussian_wavepacket(x0, sig, k0, x):\n",
" \"\"\" Helper function to take the the spread (sig), mean wavenumber (k0), mean position value (x0) \n",
" and grid of x-values (x) and return the corresponding Gaussian wavepacket values evaluated on that grid\n",
" \n",
" Arguments\n",
" ---------\n",
" x0 : float\n",
" mean position of the wavepacket in atomic units \n",
" sig : float\n",
" spread of the wavepacket in atomic units\n",
" k0 : float\n",
" the mean wavenumber in atomic units of the wavepacket, related to\n",
" mean momentum by \\hbar k0\n",
" x : float (or numpy array of floats)\n",
" the position variable for the energy eigenstate in atomic units\n",
" \"\"\"\n",
" ci = 0+1j # <== the imaginary unit\n",
" pre = 1 / (sig * np.sqrt(2 * np.pi)) # <== prefactor\n",
" gauss = np.exp(-0.5 * ((x - x0) / sig) ** 2) # <== gaussian envelope\n",
" plane = np.exp(ci * k0 * x) # <== plane-wave component\n",
" return pre * gauss * plane\n",
" "
]
},
{
"cell_type": "markdown",
"metadata": {},
"source": [
"#### Step 3: Set up plot and function parameters and form the Gaussian Wavepacket at time t=t0"
]
},
{
"cell_type": "code",
"execution_count": 4,
"metadata": {
"colab": {},
"colab_type": "code",
"id": "QLRBwgFqdr83"
},
"outputs": [
{
"name": "stderr",
"output_type": "stream",
"text": [
"/opt/anaconda3/envs/escip/lib/python3.10/site-packages/matplotlib/cbook/__init__.py:1369: ComplexWarning: Casting complex values to real discards the imaginary part\n",
" return np.asarray(x, float)\n"
]
},
{
"data": {
"text/plain": [
"[