TDSE Solution - Barrier Penetration and Tunneling for an Initially Gaussian Wavepacket

Program 1

Barrier Penetration and Tunneling for an Initially Gaussian Wavepacket

Python 3.11 Ready

import numpy as np
import matplotlib.pyplot as plt
from scipy.integrate import trapezoid

# Parameters
L, x0, k, hbar, m, a, N, V0 = 40.0, 4.0, 5.0, 1, 1, 2.0,150, 14
xs, h = np.linspace(0, L, N, retstep=True)
us = np.exp(-(xs-x0)**2/a) * np.exp(1.j*k*xs)
Vs = xs[1:-1]*0.0
Vs[35:42] = V0
Vs1 = np.zeros(N)
Vs1[35:42] = 0.5

dt = 0.2
mu = (0. + 1.j)*dt/(2*h**2)
t, tf = 0, 10.0
E = (hbar**2/2.0/m)*(k**2 + 0.5/a)
print('wavelength energy:', E)
print('Barrier height:', V0)

# Evolution matrices
A1 = np.eye(N-2)*(1.0 + 2.0*mu + 1.j*dt/2*Vs) - mu*np.eye(N-2, k=1) - mu*np.eye(N-2, k=-1)
A2 = np.eye(N-2)*(1.0 - 2.0*mu + 1.j*dt/2*Vs) + mu*np.eye(N-2, k=1) + mu*np.eye(N-2, k=-1)
B = np.dot(np.linalg.inv(A1), A2)

for frame in range(int(tf/dt)):
    us[1:-1] = np.dot(B, us[1:-1])
    t += dt
    normus = us / np.sqrt(trapezoid(us * np.conj(us), xs))
    plt.clf()  
    plt.title(f'Energy={E:.2f}, Barrier Height={V0}')
    plt.xlim(xs[0], xs[-1])
    plt.ylim(-0.6, 0.6)
    # Plot potential
    plt.plot(xs, Vs1, label='V(x)')
    # Plot real part
    plt.plot(xs, np.real(normus), label='Real part')
    # Plot probability density
    plt.plot(xs, abs(normus)**2, label='|psi|Â²')

plt.legend()
    plt.pause(0.01)

Click "Run Code" to execute script and render figures.

Mathematical Problem Formulation

Simulating Time Evolution of a Gaussian Wave Packet in a Potential Well

Consider a one-dimensional system where a quantum particle is evolving in time under the influence of a potential barrier. The time evolution of the wave function \( \psi(x, t) \) is governed by the time-dependent SchrÃ¶dinger equation:

\[ i\hbar \frac{\partial \psi(x,t)}{\partial t} = \left( -\frac{\hbar^2}{2m} \frac{\partial^2}{\partial x^2} + V(x) \right) \psi(x,t) \]

take: \( m = 1.0 \) and \( \hbar = 1.0 \)

The corresponding initial wave function is:

\[ \psi(x, 0) = A \exp\left( -\frac{(x - x_0)^2}{2\sigma^2} \right) e^{ik_0 x} \]

where \( A \) is a normalization constant.

The energy of the wave packet is given by:

\[ E = \frac{\hbar^2}{2m} \left( k_0^2 + \frac{1}{2\sigma^2} \right) \]

Crank-Nicolson Code to simulate the time evolution of the wave packet using finite difference methods, and plot the real and imaginary parts of the wave function \( \psi(x, t) \) as well as the probability density \( |\psi(x,t)|^2 \) at various time steps.

Program 2

Barrier Penetration and Tunneling for an Initially Gaussian Wavepacket

Python 3.11 Ready

import numpy as np
import matplotlib.pyplot as plt
from scipy.integrate import trapezoid

# Parameters
L, x0, k, hbar, m, a, N, V0 = 40.0, 4.0, 5.0, 1, 1, 2.0,150, 14
xs, h = np.linspace(0, L, N, retstep=True)
us = np.exp(-(xs-x0)**2/a) * np.exp(1.j*k*xs)
Vs = xs[1:-1]*0.0
Vs[35:50] = V0
Vs1 = np.zeros(N)
Vs1[35:50] = 0.5

dt = 0.1
mu = (0. + 1.j)*dt/(2*h**2)
t, tf = 0, 10.0
E = (hbar**2/2.0/m)*(k**2 + 0.5/a)
print('wavelength energy:', E)
print('Barrier height:', V0)

# Evolution matrices
A1 = np.eye(N-2)*(1.0 + 2.0*mu + 1.j*dt/2*Vs) - mu*np.eye(N-2, k=1) - mu*np.eye(N-2, k=-1)
A2 = np.eye(N-2)*(1.0 - 2.0*mu + 1.j*dt/2*Vs) + mu*np.eye(N-2, k=1) + mu*np.eye(N-2, k=-1)
B = np.dot(np.linalg.inv(A1), A2)

for frame in range(int(tf/dt)):
    us[1:-1] = np.dot(B, us[1:-1])
    t += dt
    normus = us / np.sqrt(trapezoid(us * np.conj(us), xs))
    plt.clf()  
    plt.title(f'Energy={E:.2f}, Barrier Height={V0}')
    plt.xlim(xs[0], xs[-1])
    plt.ylim(-0.6, 0.6)
    # Plot potential
    plt.plot(xs, Vs1, label='V(x)')
    # Plot real part
    plt.plot(xs, np.real(normus), label='Real part')
    # Plot probability density
    plt.plot(xs, abs(normus)**2, label='|psi|Â²')

plt.legend()
    plt.pause(0.001)

Click "Run Code" to execute script and render figures.

Program 3

Barrier Penetration and Tunneling for an Initially Gaussian Wavepacket

Python 3.11 Ready

import numpy as np
import matplotlib.pyplot as plt
from scipy.integrate import trapezoid

# Parameters
L, x0, k, hbar, m, a, N, V0 = 40.0, 4.0, 5.8, 1, 1, 2.0,150, 10
xs, h = np.linspace(0, L, N, retstep=True)
us = np.exp(-(xs-x0)**2/a) * np.exp(1.j*k*xs)
Vs = xs[1:-1]*0.0
Vs[35:42] = V0
Vs1 = np.zeros(N)
Vs1[35:42] = 0.5

dt = 0.1
mu = (0. + 1.j)*dt/(2*h**2)
t, tf = 0, 10.0
E = (hbar**2/2.0/m)*(k**2 + 0.5/a)
print('wavelength energy:', E)
print('Barrier height:', V0)

# Evolution matrices
A1 = np.eye(N-2)*(1.0 + 2.0*mu + 1.j*dt/2*Vs) - mu*np.eye(N-2, k=1) - mu*np.eye(N-2, k=-1)
A2 = np.eye(N-2)*(1.0 - 2.0*mu + 1.j*dt/2*Vs) + mu*np.eye(N-2, k=1) + mu*np.eye(N-2, k=-1)
B = np.dot(np.linalg.inv(A1), A2)

for frame in range(int(tf/dt)):
    us[1:-1] = np.dot(B, us[1:-1])
    t += dt
    normus = us / np.sqrt(trapezoid(us * np.conj(us), xs))
    plt.clf()  
    plt.title(f'Energy={E:.2f}, Barrier Height={V0}')
    plt.xlim(xs[0], xs[-1])
    plt.ylim(-0.6, 0.6)
    # Plot potential
    plt.plot(xs, Vs1, label='V(x)')
    # Plot real part
    plt.plot(xs, np.real(normus), label='Real part')
    # Plot probability density
    plt.plot(xs, abs(normus)**2, label='|psi|Â²')

plt.legend()
    plt.pause(0.001)

Click "Run Code" to execute script and render figures.

Mathematical Problem Formulation

Barrier penetration and tunneling for E > V0

Ch. 19.5 Barrier Penetration and Tunneling for an Initially Gaussian Wavepacket

Barrier Penetration and Tunneling for an Initially Gaussian Wavepacket

Mathematical Problem Formulation

Simulating Time Evolution of a Gaussian Wave Packet in a Potential Well

Barrier Penetration and Tunneling for an Initially Gaussian Wavepacket

Barrier Penetration and Tunneling for an Initially Gaussian Wavepacket

Mathematical Problem Formulation