Program 1

Barrier Penetration and Tunneling for an Initially Gaussian Wavepacket

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.

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)

Output 1

Program 2

Barrier Penetration and Tunneling for an Initially Gaussian Wavepacket

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)

Output 2

Program 3

Barrier Penetration and Tunneling for an Initially Gaussian Wavepacket

Barrier penetration and tunneling for E > V0

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)

Output 3

TDSE Solution

Barrier Penetration and Tunneling for an Initially Gaussian Wavepacket

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