Dynamics of particle in a box with a delta potential

Dynamics of particle in a box with a delta potential#

Jay Foley and Lane Tolley, University of North Carolina Charlotte

Objectives#

To demonstrate the use of the numerical methods, including the Euler method and the Runge-Kutta method, to solve the time-dependent Schrodinger equation
To illustrate the dynamics of a particle in a box with a delta potential

Learning Outcomes#

By the end of this workbook, students should be able to

Express the time-dependent Schrodinger equation using matrices for the Hamiltonian operator and vectors for the wavefunction
Utilize matrix-vector operations in numpy to compte time-derivatives of the wavefunction
Utilize numerical integration methods, particularly the Euler method and the Runge-Kutta method, to evolve the wavefunction in time.

Summary#

We will investigate the particle in a box with a delta potential that was explored in the notebook on the Linear Variational Method, 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 the particle in a box states, and to express the wavefunction in terms of the amplitudes of each of these basis states.

Step 0: Import libraries and initialize a plot object for the animation#

import numpy as np
import matplotlib.pyplot as plt
import math

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 box size and delta function
x0 = 5
L = 10
Amp = np.sqrt(2/L)

### grid of x values for the functions
x = np.linspace(0, L, 500)

### parameters for plot
ax.set_xlim((0, L))
ax.set_ylim((-2*Amp, 2*Amp))

line, = ax.plot([], [], lw=2)

# initialization function: plot the background of each frame
def init():
    line.set_data([], [])
    return (line,)
plt.show()

Approach#

In general, the dynamics of a quantum system are determined by the time-dependent Schrodinger equation (TDSE), $$ i\hbar \frac{d}{dt}\Psi(x, t) = \hat{H} \Psi(x, t). \tag{1} $$

Just as before can write the Hamiltonian operator for this system between $x=0$ and $x=10$ as follows:

\[\hat{H} = -\frac{\hbar^2}{2m} \frac{d^2}{dx^2} + \delta(x-5). \tag{2} \]

Also as before, we can build a matrix representation of this Hamiltonian $ {\bf H}$ by introducing a basis. We will use the energy eigenfunctions of the ordinary particle in a box as this basis, where these basis functions are given by

\[ \psi_n(x) = \sqrt{ \frac{2}{10} } {\rm sin}\left( \frac{n \pi x}{10} \right), \tag{3} \]

and the elements of $ {\bf H}$ will be given by

(19)#\[\begin{equation} H_{nm} = \int_0^{10} \psi^*_n(x) \hat{H} \psi_m(x) dx. \tag{4} \end{equation}\]

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:

\[ \Psi(x, t) = \sum_n c_n(t) \psi_n(x) \tag{5} \]

where we are denoting that the basis functions are time-independent, and the time-dependence is contained only in 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.
This suggests that the TDSE can be recast as $$ i\hbar \frac{d}{dt} {\bf c}(t) = {\bf H} {\bf c}(t). \tag{6} $$

Numerically Solving the TDSE Part 1: Euler Method#

The key idea for numerically solving the TDSE is that if we know the coefficients at one instant in time, ${\bf c}(t_0)$ and wish to find the coefficients at a later time $t_0 + \Delta t$, we can approximate ${\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

\[ {\bf c}(t_0 + \Delta t) = {\bf c}(t_0) + \frac{d}{dt} {\bf c}(t_0) \Delta t, \tag{7} \]

where from the TDSE, we can see

\[ \frac{d}{dt} {\bf c}(t_0) = -\frac{i}{\hbar} {\bf H} {\bf c}(t_0). \tag{8} \]

Numerically Solving the TDSE Part 2: Runge-Kutta Method#

A more numerically stable approach for updating the wavefunction can be obtained by applying the 4$^{\rm th}$-order Runge-Kutta Method. In this approach, we approximate the ${\bf c}$ vector at a future time in the following way:

\[ {\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} \]

where

(20)#\[\begin{align} k_1 &= -\frac{i}{\hbar} {\bf H}(t_0) {\bf c}(t_0) \\ 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} 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) \\ k_4 &= -\frac{i}{\hbar} {\bf H}\left(t_0 + \Delta t\right) \left( {\bf c}(t_0) + k_3 \Delta t \right). \\ \end{align}\]

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)$.

Each $k_n$ term can be seen as a particular time derivative that results from the action of the Hamiltonian on an estimate of the ${\bf c}$ vector at some instant in time. We can see that the $k_1$ is exactly equal to Eq. (8), that is it results from the action of ${\bf H}$ on ${\bf c}$ at the current instant in time $t_0$.
In the case that ${\bf H}$ is time-independent, then $k_2$ is simply defined as the action of ${\bf H}$ on a ${\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.

Questions Part 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).
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.
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).
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.

def energy_eigenvalue(n, L, m):
    """ Helper function to take the quantum number n of the particle in a box, the length 
        of the box L, and the mass of the particle m and return the energy eigenvalue in atomic units.
        Both the length and mass 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 : float
        the mass of the particle in atomic units
    """
        
    return n ** 2 * np.pi ** 2 / ( 2 * m * L ** 2)

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
        
    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)

def Hamiltonian_matrix_element(n, m):
    """ Helper function to take two indices (n, m) and return the Hamiltonian matrix element H_{n,m}, 
        which is <n|H|m> in Dirac's Bra-Ket notation.  This will call Kinetic_matrix_element(n, m) and
        Potential_matrix_element(n, m).
    
    Arguments
    ---------
    n : int
        the index corresponding to the bra state
    m : int
        the index corresponding to the ket state 
        
    Returns
    -------
    H_mn : float
         the matrix element corresponding to <n|H|m> = <n|T|m> + <n|V|m>
         
    Example
    -------
    >>> H_12 = Hamiltonian_matrix_element(1,2)
    >>> H_33 = Hamiltonian_matrix_element(3,3)
    
    Notes
    -----
    if n == m, you need to include kinetic and potential contributions of <n|H|n>.
    e.g. <n|H|n> = <n|T|n> + <n|V|n> = <n|E_n|n> + <n|V|n> = E_n<n|n> + <n|V|n> = E_n + <n|V|n>
    where E_n is the energy eigenvalue of the ordinary particle in a box state \psi_n(x)
    
    if n != m, you need to include only the potential contribution <n|V|m>.
    
    """
    H_nm = Kinetic_matrix_element(n, m) + Potential_matrix_element(n, m)
    return H_nm

def Kinetic_matrix_element(n, m):
    """ Function to take two indices (n, m) and return the Kinetic Energy matrix element T_{n,m}, 
        which is <n|T|m> in Dirac's Bra-Ket notation.  
    
    Arguments
    ---------
    n : int
        the index corresponding to the bra state
    m : int
        the index corresponding to the ket state 
        
    Returns
    -------
    T_nm : float
         the matrix element corresponding to <n|T|m> = E_m <n|m>
         
    Example
    -------
    >>> T_12 = Kinetic_matrix_element(1,2) -> 0
    >>> T_33 = Kinetic_matrix_element(3,3) -> E_3
    
    """
    T_nm = energy_eigenvalue(m, 10, 1) * (n == m)
    return T_nm

def Potential_matrix_element(n, m):
    """ Function to take two indices (n, m) and return the Hamiltonian matrix element V_{n,m}, 
        which is <n|V|m> in Dirac's Bra-Ket notation and V = \delta(x-5)
    
    Arguments
    ---------
    n : int
        the index corresponding to the bra state
    m : int
        the index corresponding to the ket state 
        
    Returns
    -------
    V_nm : float
         the matrix element corresponding to <n|V|m> = \psi_n(5) * psi_m(5)
         
    Example
    -------
    >>> V_12 = Potential_matrix_element(1,2)
    >>> V_33 = Potential_matrix_element(3,3)
    """
    V_nm = energy_eigenfunction(n, 10, 5) * energy_eigenfunction(m, 10, 5)
    return V_nm

def build_matrices(n_basis):
    """ Function that will build and return the Hamiltonian, Kinetic, and Potential energy matrices
    
    Arguments
    ---------
    n_basis : int
        the number of basis functions used to build the matrices
        
    Returns
    -------
    H : n_basis x n_basis array of floats
         the total Hamiltonian matrix
         
    T : n_basis x n_basis array of floats
         the kinetic energy matrix
         
    V : n_basis x n_basis array of floats
         the potential energy matrix 
         
    Example
    -------
    >>> H, T, V = build_matrices(3)
    """
    H = np.zeros((n_basis, n_basis))
    T = np.zeros((n_basis, n_basis))
    V = np.zeros((n_basis, n_basis))
    
    for i in range(n_basis):
        n = i + 1 # i starts from 0 but n should start at 1
        for j in range(n_basis):
            m = j + 1 # j starts at 0 but m should start at 1
            H[i, j] = Hamiltonian_matrix_element(n, m)
            T[i, j] = Kinetic_matrix_element(n, m)
            V[i, j] = Potential_matrix_element(n, m)
    return H, T, V

def compute_c_dot(H, c):
    """ Function that will compute d/dt c(t_0) using -i/hbar H * c 
    
    Arguments
    ---------
    H : n_basis x n_basis numpy array of floats
        the Hamiltonian matrix in atomic units
        
    c : 1 x n_basis numpy array of complex floats
        the vector of coefficients that defines the time-dependent wavefunction
    
    Returns
    -------
    c_dot : 1 x n_basis numpy array of complex floats
        the vector of time-derivatives of the coefficients
    
    """
    # <== define imaginary unit as ci = 0+1j
    ci = 0 + 1j
    # <== compute matrix-vector product of H and c to get c_dot
    c_dot = (-ci) * np.dot(H, c)
    return c_dot

def euler_update(H, c, dt):
    """ Function that will take c(t0) and H and return c(t0 + dt)
    
    Arguments
    ---------
    H : n_basis x n_basis numpy array of floats
        the Hamiltonian matrix in atomic units
        
    c : 1 x n_basis numpy array of complex floats
        the vector of coefficients that defines the time-dependent wavefunction
        
    dt : float
        the increment in time in atomic units
    
    
    Returns
    -------
    c_new : 1 x n_basis numpy array of complex floats
        the vector of coefficients at time t0 + dt
    
    """
    # <== call compute_c_dot(H, c) to compute c_dot
    c_dot = compute_c_dot(H, c)
    # <== compute c_new using c, c_dot, and dt
    c_new = c + c_dot * dt
    return c_new

def rk4_update(H, c, dt):
    """ Function that will take c(t0) and H and return c(t0 + dt)
    
    Arguments
    ---------
    H : n_basis x n_basis numpy array of floats
        the Hamiltonian matrix in atomic units
        
    c : 1 x n_basis numpy array of complex floats
        the vector of coefficients that defines the time-dependent wavefunction
        
    dt : float
        the increment in time in atomic units
    
    
    Returns
    -------
    c_new : 1 x n_basis numpy array of complex floats
        the vector of coefficients at time t0 + dt
    
    """
    # <== define k_1, k_2, k_3, k_4
    ci = 0 + 1j
    k_1 = -ci * np.dot(H, c) 
    k_2 = -ci * np.dot(H, (c + k_1 * dt / 2))
    k_3 = -ci * np.dot(H, (c + k_2 * dt / 2))
    k_4 = -ci * np.dot(H, (c + k_3 * dt))
    # <== define c_new as 1/6 * (k_1 + 2 * k_2 + 2 * k_3 + k_4) * dt
    c_new = c + (1 / 6) * (k_1 + 2 * k_2 + 2 * k_3 + k_4) * dt
    return c_new

Test for rk4_update function#

Assume initial condition $\Psi(x,t_0) = \psi_1(x)e^{-iE_1 t_0/\hbar}$ and $t_0 = 0$ and update using RK4 to $t_f = 1$ atomic units with a time-step of $\Delta t = 0.01$. We will feed the kinetic energy matrix to the RK4 method so that the time evolution should correspond give us $\Psi(x,t_f) = \psi_1(x)e^{-iE_1 t_f/\hbar}$ and $t_f = 1$.

# imaginary unit
ci = 0+1j

# E1 
E1 = energy_eigenvalue(1, 10, 1)

# dt
dt = 0.01

# tf
tf = 1

# define initial c vector for RK4 update
c_rk4 = np.array([1, 0], dtype=complex)

# number of basis functions
n_basis = 2

# build H, T, and V
Hamiltonian_matrix, Kinetic_matrix, Potential_matrix = build_matrices(n_basis)

# Define analytical result for c(t_f)
c_analytical = np.array([np.exp(-ci * E1 * tf  ), 0+0j])

# update c_rk4 100 times using the RK4 method
for i in range(1, 101):
    c_rk4 = rk4_update(Kinetic_matrix, c_rk4, dt)


assert np.allclose(c_analytical, c_rk4)
print("PASSED!")

PASSED!

Time-dependence with the delta potential#

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. 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

# number of basis functions
n_basis = 70

# build H, T, and V
Hamiltonian_matrix, Kinetic_matrix, Potential_matrix = build_matrices(n_basis)

# number of time steps
N_time = 500

# create the c_vec array as a n_basis x N_time numpy array of complex floats
c_vec = np.zeros((n_basis, N_time), dtype=complex)

# define a time increment of 0.005 atomic units
dt = 0.01

# initialize the c(0) vector as (1+0j, 0, 0, ...) -> initial state = \psi_1(x)
c_vec[0,0] = 1+0j

# loop over time-steps and update the c-vector at each time
for i in range(1, N_time):
    c_vec[:,i] = rk4_update(Hamiltonian_matrix, c_vec[:,i-1], dt)

The following lines will animate the wavefunction as shown in Eq. (5). Now, the elements of ${\bf c}(t)$ are stored in c_vec, and the $\psi_n(x)$ functions can be obtained by calling energy_eigenfunction(n, L, x).

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, n_basis):
        fx = energy_eigenfunction(j+1, 10, x)
        y = y + c_vec[j,i] * fx
    t[i] = i
    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