{
"cells": [
{
"attachments": {},
"cell_type": "markdown",
"metadata": {},
"source": [
"# Dynamics of particle in a box with a delta potential\n",
"[Jay Foley and Lane Tolley, University of North Carolina Charlotte](https://foleylab.github.io/)\n",
"\n",
"#### Objectives\n",
"- To demonstrate the use of the numerical methods, including the Euler method and the Runge-Kutta method, to solve the time-dependent Schrodinger equation\n",
"- To illustrate the dynamics of a particle in a box with a delta potential\n",
"\n",
"#### Learning Outcomes\n",
"By the end of this workbook, students should be able to\n",
"- Express the time-dependent Schrodinger equation using matrices for the Hamiltonian operator and vectors for the wavefunction\n",
"- Utilize matrix-vector operations in `numpy` to compte time-derivatives of the wavefunction\n",
"- Utilize numerical integration methods, particularly the Euler method and the Runge-Kutta method, to evolve the wavefunction in time.\n",
"\n",
"\n",
"#### Summary\n",
"We will investigate the particle in a box with a delta potential that was explored in the notebook on the [Linear Variational Method](https://escip.io/notebooks/linear_variational_method.html), but in the context of the time-dependence induced by the delta potential on an initial state, rather than through the application of the variational method to identify stationary states of this Hamiltonian. We will choose an initial state to be the ground state of the ordinary particle in a box system, so again we will find it useful to build the Hamiltonian matrix in the basis of\n",
"the particle in a box states, and to express the wavefunction in terms of the amplitudes of each of these basis states.\n",
"\n",
"\n"
]
},
{
"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",
"import math\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 box size and delta function\n",
"x0 = 5\n",
"L = 10\n",
"Amp = np.sqrt(2/L)\n",
"\n",
"### grid of x values for the functions\n",
"x = np.linspace(0, L, 500)\n",
"\n",
"### parameters for plot\n",
"ax.set_xlim((0, L))\n",
"ax.set_ylim((-2*Amp, 2*Amp))\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,)\n",
"plt.show()"
]
},
{
"cell_type": "markdown",
"metadata": {},
"source": [
"# Approach\n",
"In general, the dynamics of a quantum system are determined by the time-dependent Schrodinger equation (TDSE),\n",
"$$ i\\hbar \\frac{d}{dt}\\Psi(x, t) = \\hat{H} \\Psi(x, t). \\tag{1} $$\n",
"\n",
"Just as [before](https://escip.io/notebooks/linear_variational_method.html) can write the Hamiltonian operator for this system between $x=0$ and $x=10$ as follows:\n",
"\n",
"$$\\hat{H} = -\\frac{\\hbar^2}{2m} \\frac{d^2}{dx^2} + \\delta(x-5). \\tag{2} $$\n",
"\n",
"Also as before, we can build a matrix representation of this Hamiltonian $ {\\bf H}$ by introducing a basis. We will use the \n",
"energy eigenfunctions of the ordinary particle in a box as this basis, where these basis functions are given by\n",
"\n",
"$$ \\psi_n(x) = \\sqrt{ \\frac{2}{10} } {\\rm sin}\\left( \\frac{n \\pi x}{10} \\right), \\tag{3} $$\n",
"\n",
"and the elements of $ {\\bf H}$ will be given by \n",
"\\begin{equation}\n",
"H_{nm} = \\int_0^{10} \\psi^*_n(x) \\hat{H} \\psi_m(x) dx. \\tag{4} \n",
"\\end{equation}\n",
"Because the basis set comprised of all $\\psi_n(x)$ is complete, we can represent our wavefunction at any instant in time as a linear combination of these basis functions:\n",
"\n",
"$$ \\Psi(x, t) = \\sum_n c_n(t) \\psi_n(x) \\tag{5} $$\n",
"where we are denoting that the basis functions are time-independent, and the time-dependence is contained only in\n",
"the expansion coefficients $c_n(t)$. Thus we can map the time-dependence of the collection of these coefficients, which we will call ${\\bf c}(t)$ to the time-dependence of the wavefunction. \n",
"This suggests that the TDSE can be recast as\n",
"$$ i\\hbar \\frac{d}{dt} {\\bf c}(t) = {\\bf H} {\\bf c}(t). \\tag{6} $$"
]
},
{
"cell_type": "markdown",
"metadata": {},
"source": [
"#### Numerically Solving the TDSE Part 1: Euler Method\n",
"The key idea for numerically solving the TDSE is that if we know the coefficients at one instant in time,\n",
"${\\bf c}(t_0)$ and wish to find the coefficients at a later time $t_0 + \\Delta t$, we can approximate \n",
"${\\bf c}(t_0 + \\Delta t)$ using information about the time derivative of the coefficients given by the TDSE. The simplest approximation is known as the Euler update, and has the form\n",
"\n",
"$$ {\\bf c}(t_0 + \\Delta t) = {\\bf c}(t_0) + \\frac{d}{dt} {\\bf c}(t_0) \\Delta t, \\tag{7} $$\n",
"\n",
"where from the TDSE, we can see\n",
"\n",
"$$ \\frac{d}{dt} {\\bf c}(t_0) = -\\frac{i}{\\hbar} {\\bf H} {\\bf c}(t_0). \\tag{8} $$"
]
},
{
"cell_type": "markdown",
"metadata": {},
"source": [
"#### Numerically Solving the TDSE Part 2: Runge-Kutta Method\n",
"A more numerically stable approach for updating the wavefunction can be obtained by applying the 4$^{\\rm th}$-order\n",
"Runge-Kutta Method. In this approach, we approximate the ${\\bf c}$ vector at a future time in the following way:\n",
"\n",
"$$ {\\bf c}(t_0 + \\Delta t) = {\\bf c}(t_0) + \\frac{1}{6}\\left(k_1 + 2k_2 + 2k_3 + k_4 \\right)\\Delta t \\tag{9} $$\n",
"\n",
"where\n",
"\n",
"\\begin{align}\n",
"k_1 &= -\\frac{i}{\\hbar} {\\bf H}(t_0) {\\bf c}(t_0) \\\\\n",
"k_2 &= -\\frac{i}{\\hbar} {\\bf H}\\left(t_0 + \\frac{\\Delta t}{2}\\right) \\left( {\\bf c}(t_0) + k_1 \\frac{\\Delta t}{2} \\right) \\\\ \\tag{10}\n",
"k_3 &= -\\frac{i}{\\hbar} {\\bf H}\\left(t_0 + \\frac{\\Delta t}{2}\\right) \\left( {\\bf c}(t_0) + k_2 \\frac{\\Delta t}{2} \\right) \\\\\n",
"k_4 &= -\\frac{i}{\\hbar} {\\bf H}\\left(t_0 + \\Delta t\\right) \\left( {\\bf c}(t_0) + k_3 \\Delta t \\right). \\\\\n",
"\\end{align}\n",
"If the Hamiltonian does not depend upon time, as is the case for the particle in a box with a delta potential, then ${\\bf H} = {\\bf H}(t_0) = {\\bf H}\\left(t_0 + \\frac{\\Delta t}{2}\\right) = {\\bf H}\\left(t_0 + \\Delta t\\right)$.\n",
"\n",
"Each $k_n$ term can be seen as a particular time derivative that results from the action of the Hamiltonian on an \n",
"estimate of the ${\\bf c}$ vector at some instant in time. We can see that the $k_1$ is exactly equal to Eq. (8), \n",
"that is it results from the action of ${\\bf H}$ on ${\\bf c}$ at the current instant in time $t_0$. \n",
"In the case that ${\\bf H}$ is time-independent, then $k_2$ is simply defined as the action of ${\\bf H}$ on a \n",
"${\\bf c}$ vector that has been updated by half an Euler step, i.e. equivalent to Eq. (7) with an update of $\\frac{\\Delta t}{2}$. More can be read about the Runge-Kutta method [here](https://en.wikipedia.org/wiki/Runge%E2%80%93Kutta_methods)."
]
},
{
"cell_type": "markdown",
"metadata": {},
"source": [
"### Questions Part 1:\n",
"\n",
"1. We will include the helper functions from before that will build the Hamiltonian matrix in the PIB basis. Add a function called `compute_c_dot(H, c)` that will take the Hamiltonian matrix and a given c vector and return the time-derivative of the c vector using the TDSE using Eq. (8).\n",
"\n",
"2. Write a function called `euler_update(H, c, dt)` that takes the Hamiltonian matrix `H`, a c vector `c` and a time increment `dt` and returns the c vector at a future time using Eq. (7). It can call `compute_c_dot(H, c)` to get the time-derivative of `c`.\n",
"\n",
"3. Write a function called `rk4_update(H, c, dt)` that takes `H`, a c vector `c` and a time increment `dt` and returns the c vector at a future time using Eq. (9).\n",
"\n",
"4. If we use only the kinetic energy matrix, then we know analytically time-dependence of any arbitrary state is given by $\\Psi(x,t) = \\sum_n c_n(t_0){\\rm e}^{-iE_n t/\\hbar} \\psi_n(x) {\\rm e}^{-iE_n t/\\hbar}$. Write a test to confirm that evolving the initial state ${\\bf c}(t_0) = (1, 0, ..., 0)$ to ${\\bf c}(0.1)$ using `rk4_update()` with the kinetic energy matrix matches the analytical result. "
]
},
{
"cell_type": "code",
"execution_count": 2,
"metadata": {},
"outputs": [],
"source": [
"def energy_eigenvalue(n, L, m):\n",
" \"\"\" Helper function to take the quantum number n of the particle in a box, the length \n",
" of the box L, and the mass of the particle m and return the energy eigenvalue in atomic units.\n",
" Both the length and mass 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",
" m : float\n",
" the mass of the particle in atomic units\n",
" \"\"\"\n",
" \n",
" return n ** 2 * np.pi ** 2 / ( 2 * m * L ** 2)\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\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)\n",
"\n",
"def Hamiltonian_matrix_element(n, m):\n",
" \"\"\" Helper function to take two indices (n, m) and return the Hamiltonian matrix element H_{n,m}, \n",
" which is in Dirac's Bra-Ket notation. This will call Kinetic_matrix_element(n, m) and\n",
" Potential_matrix_element(n, m).\n",
" \n",
" Arguments\n",
" ---------\n",
" n : int\n",
" the index corresponding to the bra state\n",
" m : int\n",
" the index corresponding to the ket state \n",
" \n",
" Returns\n",
" -------\n",
" H_mn : float\n",
" the matrix element corresponding to = + \n",
" \n",
" Example\n",
" -------\n",
" >>> H_12 = Hamiltonian_matrix_element(1,2)\n",
" >>> H_33 = Hamiltonian_matrix_element(3,3)\n",
" \n",
" Notes\n",
" -----\n",
" if n == m, you need to include kinetic and potential contributions of .\n",
" e.g. = + = + = E_n + = E_n + \n",
" where E_n is the energy eigenvalue of the ordinary particle in a box state \\psi_n(x)\n",
" \n",
" if n != m, you need to include only the potential contribution .\n",
" \n",
" \"\"\"\n",
" H_nm = Kinetic_matrix_element(n, m) + Potential_matrix_element(n, m)\n",
" return H_nm\n",
"\n",
"def Kinetic_matrix_element(n, m):\n",
" \"\"\" Function to take two indices (n, m) and return the Kinetic Energy matrix element T_{n,m}, \n",
" which is in Dirac's Bra-Ket notation. \n",
" \n",
" Arguments\n",
" ---------\n",
" n : int\n",
" the index corresponding to the bra state\n",
" m : int\n",
" the index corresponding to the ket state \n",
" \n",
" Returns\n",
" -------\n",
" T_nm : float\n",
" the matrix element corresponding to = E_m \n",
" \n",
" Example\n",
" -------\n",
" >>> T_12 = Kinetic_matrix_element(1,2) -> 0\n",
" >>> T_33 = Kinetic_matrix_element(3,3) -> E_3\n",
" \n",
" \"\"\"\n",
" T_nm = energy_eigenvalue(m, 10, 1) * (n == m)\n",
" return T_nm\n",
"\n",
"def Potential_matrix_element(n, m):\n",
" \"\"\" Function to take two indices (n, m) and return the Hamiltonian matrix element V_{n,m}, \n",
" which is in Dirac's Bra-Ket notation and V = \\delta(x-5)\n",
" \n",
" Arguments\n",
" ---------\n",
" n : int\n",
" the index corresponding to the bra state\n",
" m : int\n",
" the index corresponding to the ket state \n",
" \n",
" Returns\n",
" -------\n",
" V_nm : float\n",
" the matrix element corresponding to = \\psi_n(5) * psi_m(5)\n",
" \n",
" Example\n",
" -------\n",
" >>> V_12 = Potential_matrix_element(1,2)\n",
" >>> V_33 = Potential_matrix_element(3,3)\n",
" \"\"\"\n",
" V_nm = energy_eigenfunction(n, 10, 5) * energy_eigenfunction(m, 10, 5)\n",
" return V_nm\n",
"\n",
"def build_matrices(n_basis):\n",
" \"\"\" Function that will build and return the Hamiltonian, Kinetic, and Potential energy matrices\n",
" \n",
" Arguments\n",
" ---------\n",
" n_basis : int\n",
" the number of basis functions used to build the matrices\n",
" \n",
" Returns\n",
" -------\n",
" H : n_basis x n_basis array of floats\n",
" the total Hamiltonian matrix\n",
" \n",
" T : n_basis x n_basis array of floats\n",
" the kinetic energy matrix\n",
" \n",
" V : n_basis x n_basis array of floats\n",
" the potential energy matrix \n",
" \n",
" Example\n",
" -------\n",
" >>> H, T, V = build_matrices(3)\n",
" \"\"\"\n",
" H = np.zeros((n_basis, n_basis))\n",
" T = np.zeros((n_basis, n_basis))\n",
" V = np.zeros((n_basis, n_basis))\n",
" \n",
" for i in range(n_basis):\n",
" n = i + 1 # i starts from 0 but n should start at 1\n",
" for j in range(n_basis):\n",
" m = j + 1 # j starts at 0 but m should start at 1\n",
" H[i, j] = Hamiltonian_matrix_element(n, m)\n",
" T[i, j] = Kinetic_matrix_element(n, m)\n",
" V[i, j] = Potential_matrix_element(n, m)\n",
" return H, T, V\n",
"\n",
"def compute_c_dot(H, c):\n",
" \"\"\" Function that will compute d/dt c(t_0) using -i/hbar H * c \n",
" \n",
" Arguments\n",
" ---------\n",
" H : n_basis x n_basis numpy array of floats\n",
" the Hamiltonian matrix in atomic units\n",
" \n",
" c : 1 x n_basis numpy array of complex floats\n",
" the vector of coefficients that defines the time-dependent wavefunction\n",
" \n",
" Returns\n",
" -------\n",
" c_dot : 1 x n_basis numpy array of complex floats\n",
" the vector of time-derivatives of the coefficients\n",
" \n",
" \"\"\"\n",
" # <== define imaginary unit as ci = 0+1j\n",
" ci = 0 + 1j\n",
" # <== compute matrix-vector product of H and c to get c_dot\n",
" c_dot = (-ci) * np.dot(H, c)\n",
" return c_dot\n",
"\n",
"def euler_update(H, c, dt):\n",
" \"\"\" Function that will take c(t0) and H and return c(t0 + dt)\n",
" \n",
" Arguments\n",
" ---------\n",
" H : n_basis x n_basis numpy array of floats\n",
" the Hamiltonian matrix in atomic units\n",
" \n",
" c : 1 x n_basis numpy array of complex floats\n",
" the vector of coefficients that defines the time-dependent wavefunction\n",
" \n",
" dt : float\n",
" the increment in time in atomic units\n",
" \n",
" \n",
" Returns\n",
" -------\n",
" c_new : 1 x n_basis numpy array of complex floats\n",
" the vector of coefficients at time t0 + dt\n",
" \n",
" \"\"\"\n",
" # <== call compute_c_dot(H, c) to compute c_dot\n",
" c_dot = compute_c_dot(H, c)\n",
" # <== compute c_new using c, c_dot, and dt\n",
" c_new = c + c_dot * dt\n",
" return c_new\n",
"\n",
"def rk4_update(H, c, dt):\n",
" \"\"\" Function that will take c(t0) and H and return c(t0 + dt)\n",
" \n",
" Arguments\n",
" ---------\n",
" H : n_basis x n_basis numpy array of floats\n",
" the Hamiltonian matrix in atomic units\n",
" \n",
" c : 1 x n_basis numpy array of complex floats\n",
" the vector of coefficients that defines the time-dependent wavefunction\n",
" \n",
" dt : float\n",
" the increment in time in atomic units\n",
" \n",
" \n",
" Returns\n",
" -------\n",
" c_new : 1 x n_basis numpy array of complex floats\n",
" the vector of coefficients at time t0 + dt\n",
" \n",
" \"\"\"\n",
" # <== define k_1, k_2, k_3, k_4\n",
" ci = 0 + 1j\n",
" k_1 = -ci * np.dot(H, c) \n",
" k_2 = -ci * np.dot(H, (c + k_1 * dt / 2))\n",
" k_3 = -ci * np.dot(H, (c + k_2 * dt / 2))\n",
" k_4 = -ci * np.dot(H, (c + k_3 * dt))\n",
" # <== define c_new as 1/6 * (k_1 + 2 * k_2 + 2 * k_3 + k_4) * dt\n",
" c_new = c + (1 / 6) * (k_1 + 2 * k_2 + 2 * k_3 + k_4) * dt\n",
" return c_new"
]
},
{
"cell_type": "markdown",
"metadata": {},
"source": [
"# Test for rk4_update function \n",
"Assume initial condition $\\Psi(x,t_0) = \\psi_1(x)e^{-iE_1 t_0/\\hbar}$ and $t_0 = 0$\n",
"and update using RK4 to $t_f = 1$ atomic units with a time-step of $\\Delta t = 0.01$. We will \n",
"feed the kinetic energy matrix to the RK4 method so that the time evolution should correspond give us\n",
"$\\Psi(x,t_f) = \\psi_1(x)e^{-iE_1 t_f/\\hbar}$ and $t_f = 1$.\n"
]
},
{
"cell_type": "code",
"execution_count": 3,
"metadata": {},
"outputs": [
{
"name": "stdout",
"output_type": "stream",
"text": [
"PASSED!\n"
]
}
],
"source": [
"# imaginary unit\n",
"ci = 0+1j\n",
"\n",
"# E1 \n",
"E1 = energy_eigenvalue(1, 10, 1)\n",
"\n",
"# dt\n",
"dt = 0.01\n",
"\n",
"# tf\n",
"tf = 1\n",
"\n",
"# define initial c vector for RK4 update\n",
"c_rk4 = np.array([1, 0], dtype=complex)\n",
"\n",
"# number of basis functions\n",
"n_basis = 2\n",
"\n",
"# build H, T, and V\n",
"Hamiltonian_matrix, Kinetic_matrix, Potential_matrix = build_matrices(n_basis)\n",
"\n",
"# Define analytical result for c(t_f)\n",
"c_analytical = np.array([np.exp(-ci * E1 * tf ), 0+0j])\n",
"\n",
"# update c_rk4 100 times using the RK4 method\n",
"for i in range(1, 101):\n",
" c_rk4 = rk4_update(Kinetic_matrix, c_rk4, dt)\n",
"\n",
"\n",
"assert np.allclose(c_analytical, c_rk4)\n",
"print(\"PASSED!\")\n"
]
},
{
"cell_type": "markdown",
"metadata": {},
"source": [
"# Time-dependence with the delta potential\n",
"Now we will propage the same initial state using the full Hamiltonian including the delta potential. We will use a larger basis this time to capture the dynamics imparted by the potential.\n",
"The following lines will build the Hamiltonian matrix, define an initial c-vector, and update it using the RK4 method for a number of time-steps storing the time-dependent c-vector in a n_basis x n_time numpy array of complex floats"
]
},
{
"cell_type": "code",
"execution_count": 4,
"metadata": {},
"outputs": [],
"source": [
"# number of basis functions\n",
"n_basis = 70\n",
"\n",
"# build H, T, and V\n",
"Hamiltonian_matrix, Kinetic_matrix, Potential_matrix = build_matrices(n_basis)\n",
"\n",
"# number of time steps\n",
"N_time = 500\n",
"\n",
"# create the c_vec array as a n_basis x N_time numpy array of complex floats\n",
"c_vec = np.zeros((n_basis, N_time), dtype=complex)\n",
"\n",
"# define a time increment of 0.005 atomic units\n",
"dt = 0.01\n",
"\n",
"# initialize the c(0) vector as (1+0j, 0, 0, ...) -> initial state = \\psi_1(x)\n",
"c_vec[0,0] = 1+0j\n",
"\n",
"# loop over time-steps and update the c-vector at each time\n",
"for i in range(1, N_time):\n",
" c_vec[:,i] = rk4_update(Hamiltonian_matrix, c_vec[:,i-1], dt)\n",
"\n"
]
},
{
"cell_type": "markdown",
"metadata": {},
"source": [
"The following lines will animate the wavefunction as shown in Eq. (5). Now, the elements of ${\\bf c}(t)$ are stored\n",
"in `c_vec`, and the $\\psi_n(x)$ functions can be obtained by calling `energy_eigenfunction(n, L, x)`. "
]
},
{
"cell_type": "code",
"execution_count": 5,
"metadata": {},
"outputs": [
{
"data": {
"text/html": [
"\n",
"\n",
"\n",
"\n",
"\n",
"\n",
"