Animation of Quantum Tunneling: A code in Python

Last Updated on August 18, 2024 by Max

In this article, we explore the phenomenon of quantum tunneling through a detailed Python simulation. Quantum tunneling is a fundamental concept in quantum mechanics where particles pass through potential barriers, defying classical physics.

The code provided demonstrates the process of quantum tunneling using a 2D Gaussian wave packet and a potential barrier, showcasing the wave-particle duality.

This animation visualizes how the wave packet evolves over time, offering an intuitive understanding of quantum mechanics through computational simulation.

Table of Contents

Animation Result

The video below shows the simulation results. The detailed method of the simulation and the simulation parameters can be found in Ref. [1]

Figure 1: Animation of quantum tunneling.

Animation Code in Python

Wave Packet Simulation


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

# Constants
hbar = 1.0545718e-34  # Reduced Planck constant (Joule second)
m_e = 9.10938356e-31  # Mass of electron (kg)
wavelength = 20e-9  # de Broglie wavelength (meters)
k0 = 2 * np.pi / wavelength  # Wave number

# Grid parameters
Lx, Ly = 1e-6, 1e-6  # Simulation domain size (meters)
Nx, Ny = 512, 512  # Number of grid points
dx, dy = Lx / Nx, Ly / Ny  # Grid spacing
x = np.linspace(-Lx / 2, Lx / 2, Nx)
y = np.linspace(-Ly / 2, Ly / 2, Ny)
X, Y = np.meshgrid(x, y)

# Time parameters
dt = 3e-14  # Time step (seconds)
Nt = 501  # Number of time steps

# Wave packet parameters
sigma_x = sigma_y = 60e-9  # Width of the Gaussian packet (meters)
x0, y0 = -2e-7, 0  # Initial position of the packet (meters)
kx0, ky0 = k0, 0  # Initial wave vector

# Barrier parameters
barrier_width = 1e-8  # Width of the square barrier (meters)
barrier_height = 1.05 * hbar**2 * k0**2 / (2 * m_e)  # Height of the barrier (Joules)

# Initial Gaussian wave packet
psi0 = np.exp(-(X - x0) ** 2 / (2 * sigma_x ** 2)) * np.exp(-(Y - y0) ** 2 / (2 * sigma_y ** 2))
psi0 = psi0.astype(np.complex128)
psi0 *= np.exp(1j * (kx0 * X + ky0 * Y))

# Normalize the initial wave packet
psi0 /= np.sqrt(np.sum(np.abs(psi0) ** 2))

# Potential energy (Square barrier in the center)
V = np.zeros_like(X)
V[np.abs(X) < barrier_width / 2] = barrier_height

# Fourier space components
kx = np.fft.fftfreq(Nx, dx) * 2 * np.pi
ky = np.fft.fftfreq(Ny, dy) * 2 * np.pi
KX, KY = np.meshgrid(kx, ky)
K2 = KX ** 2 + KY ** 2

# Split-step Fourier method
psi = psi0.copy()

# Prepare for animation
fig, ax = plt.subplots(figsize=(16, 9))

def update(t):
    global psi
    # (a) 1/2 Evolution for the potential energy in real space
    psi *= np.exp(-1j * V * dt / (2 * hbar))
    
    # (b) Forward transform
    psi_k = np.fft.fft2(psi)
    
    # (c) Full evolution for the kinetic energy in Fourier space
    psi_k *= np.exp(-1j * hbar * K2 * dt / (2 * m_e))
    
    # (d) Inverse Fourier transform
    psi = np.fft.ifft2(psi_k)
    
    # (e) 1/2 Evolution for the potential energy in real space
    psi *= np.exp(-1j * V * dt / (2 * hbar))

    # Normalize the wave function
    psi /= np.sqrt(np.sum(np.abs(psi) ** 2) * dx * dy)

    # Visualization of wave packet evolution
    ax.clear()
    ax.set_xlabel(r'$x$ (m)', fontsize=24)
    ax.set_ylabel(r'$y$ (m)', fontsize=24)
    ax.tick_params(which="major", axis="both", direction="in", top=True, right=True, length=5, width=1, labelsize=20)
    normalized_psi = np.abs(psi) ** 2 / np.max(np.abs(psi) ** 2)
    normalized_V = V / np.max(V)
    combined = normalized_psi + 0.5 * normalized_V
    im = ax.imshow(combined, extent=(-Lx / 2, Lx / 2, -Ly / 2, Ly / 2), cmap='hot', animated=True)
    ax.set_title(f'Time step: {t}', fontsize=20)

    return [im]

ani = animation.FuncAnimation(fig, update, frames=Nt, repeat=False)

# Save the animation
ani.save('quantum_tunneling_simulation.mp4', writer='ffmpeg', fps=20)

plt.show()

Code Overview

Here are the steps in the animation code -

Initialization of Constants: The code starts by defining key physical constants, including the reduced Planck constant (hbar), electron mass (m_e), and the de Broglie wavelength (wavelength). These are essential for calculating the quantum properties of the system, such as the wave number (k0).

Grid and Simulation Setup: The simulation domain size (Lx, Ly) is set to 1 micron in both dimensions. The grid resolution is defined by the number of points (Nx, Ny), and the spatial grid is created using np.linspace to divide the domain into evenly spaced intervals. A meshgrid (X, Y) is generated to facilitate 2D calculations.

Time and Wave Packet Parameters: The time evolution of the system is controlled by the time step (dt) and the total number of time steps (Nt). The initial Gaussian wave packet is defined with specific width parameters (sigma_x, sigma_y), and positioned at (x0, y0). The wave vector components (kx0, ky0) determine the packet's initial momentum.

Barrier Setup: A potential barrier is introduced in the center of the simulation domain. The barrier is characterized by its width (barrier_width) and height (barrier_height), which are calculated using quantum mechanical principles.

Initial Wave Packet Creation: The initial wave packet is constructed as a complex Gaussian function, incorporating both amplitude and phase factors. The packet is then normalized to ensure the total probability is unity.

Fourier Space Preparation: Fourier space components (kx, ky) are calculated to enable the momentum space evolution of the wave packet. The meshgrid (KX, KY) is used to compute the squared wave numbers (K2), necessary for kinetic energy calculations.

Wave Packet Evolution: The wave packet evolves over time using the split-step Fourier method. This involves alternating between real space, where the potential energy operator is applied, and Fourier space, where the kinetic energy operator is used. The method ensures accurate simulation of quantum dynamics.

Normalization: After each time step, the wave packet is renormalized to maintain the total probability, preventing numerical errors from accumulating.

Animation Setup: The update function is defined to handle the visualization for each time step. It updates the wave packet, applies the necessary transformations, and prepares the data for plotting.

Animation Creation: The FuncAnimation function from matplotlib.animation is used to create the animation. It repeatedly calls the update function for each time step, generating a sequence of frames that depict the evolution of the wave packet.
Saving and Displaying the Animation: The animation is saved as an MP4 file using ani.save(), and then displayed with plt.show() for real-time visualization of the quantum tunneling process.

Conclusion

In conclusion, this Python code provides a detailed simulation of quantum tunneling, a fundamental phenomenon in quantum mechanics. By visualizing the evolution of a 2D Gaussian wave packet as it interacts with a potential barrier, the animation offers an intuitive understanding of how particles can penetrate barriers that classical physics deems impassable. The step-by-step construction and execution of the code demonstrate the application of the split-step Fourier method, making this a valuable tool for exploring and teaching quantum concepts.

References

[1] Max, "Quantum Tunneling: Theoretical Insights and Numerical Simulation", MatterWaveX.Com, August 13, (2024).

Max

I am a science enthusiast and writer, specializing in matter-wave optics and related technologies. My goal is to promote awareness and understanding of these advanced fields among students and the general public.

Animation Result

Animation Code in Python

Code Overview

Conclusion

References

Leave a Comment Cancel Reply

Related Posts