Last Updated on December 10, 2024 by Max
Numerical simulations of quantum systems provide a robust framework for exploring particle dynamics in varying potential landscapes. We simulate the motion of a quantum particle in a two-dimensional harmonic potential by solving the Time-Dependent Schrödinger Equation (TDSE) using the Split-Step Fourier Method (SSFM), an efficient scheme for evolving the wave function in real and Fourier spaces.
This study provides key insights into the quantum dynamics of harmonic confinement, relevant to systems like trapped ions and optical lattices, which are foundational to quantum computing and precision measurement technologies.
Numerical Methods
A detailed explanation of solving the TDSE has been provided earlier [1]. For the sake of completeness, a concise summary of the method is presented below.
Time-Dependent Schrödinger Equation (TDSE)
The TDSE in two dimensions governs the evolution of the wave function \( \psi(x, y, t) \) and is expressed as:
\[i\hbar \frac{\partial \psi}{\partial t} = -\frac{\hbar^2}{2m_e} \nabla^2 \psi + V(x, y) \psi,\]
where \( \hbar \) is the reduced Planck constant, \( m_e \) is the particle’s mass, and \( V(x, y) \) is the potential energy.
Initial Wave Packet
The initial wave function is a Gaussian packet defined as,
\[\psi(x, y, t=0) = \exp\left[-\frac{(x – x_0)^2}{2\sigma_x^2} – \frac{(y – y_0)^2}{2\sigma_y^2}\right] e^{i(k_{x0}x + k_{y0}y)},\]
where \( \sigma_x, \sigma_y \) are the widths, \( (x_0, y_0) \) is the initial position, and \( (k_{x0}, k_{y0}) \) is the initial wave vector.
Wave Function Normalization
The wave function is normalized at each step to ensure conservation of total probability,
\[\psi(x, y, t) \rightarrow \frac{\psi(x, y, t)}{\sqrt{\int |\psi(x, y, t)|^2 dx dy}}.\]
Split-Step Fourier Method
The TDSE is solved iteratively using the Split-Step Fourier Method (SSFM), which alternates between real and Fourier space:
1. Real-Space Evolution (Potential Energy):
\[ \psi(x, y, t + \frac{\Delta t}{2}) = \psi(x, y, t) e^{-\frac{i V(x, y) \Delta t}{2\hbar}}. \]
2. Fourier Transform to Momentum Space:
\[ \tilde{\psi}(k_x, k_y, t) = \mathcal{F}[\psi(x, y, t)]. \]
3. Momentum-Space Evolution (Kinetic Energy):
\[ \tilde{\psi}(k_x, k_y, t + \Delta t) = \tilde{\psi}(k_x, k_y, t) e^{-\frac{i \hbar (k_x^2 + k_y^2) \Delta t}{2m_e}}. \]
4. Inverse Fourier Transform to Real Space:
\[ \psi(x, y, t + \Delta t) = \mathcal{F}^{-1}[\tilde{\psi}(k_x, k_y, t + \Delta t)]. \]
5. Final Real-Space Evolution (Potential Energy):
\[ \psi(x, y, t + \Delta t) = \psi(x, y, t + \frac{\Delta t}{2}) e^{-\frac{i V(x, y) \Delta t}{2\hbar}}. \]
This method efficiently combines the effects of potential and kinetic energies, providing accurate simulations of quantum systems.
| Parameter | Value |
|---|---|
| Reduced Planck constant (\( \hbar \)) | \( 1.0545718 \times 10^{-34} \, \text{J·s} \) |
| Mass of electron (\( m_e \)) | \( 9.10938356 \times 10^{-31} \, \text{kg} \) |
| De Broglie wavelength (\( \lambda \)) | \( 20 \times 10^{-9} \, \text{m} \) |
| Wave number (\( k_0 \)) | \( \frac{2\pi}{\lambda} \) |
| Simulation domain size (\( L_x, L_y \)) | \( 1 \times 10^{-6} \, \text{m} \) |
| Number of grid points (\( N_x, N_y \)) | 512 |
| Grid spacing (\( \Delta x, \Delta y \)) | \( \frac{L_x}{N_x}, \frac{L_y}{N_y} \) |
| Time step (\( \Delta t \)) | \( 3 \times 10^{-14} \, \text{s} \) |
| Number of time steps (\( N_t \)) | 501 |
| Gaussian packet width (\( \sigma_x, \sigma_y \)) | \( 60 \times 10^{-9} \, \text{m} \) |
| Initial packet position (\( x_0, y_0 \)) | \( -2 \times 10^{-7} \, \text{m}, 0 \, \text{m} \) |
| Initial wave vector (\( k_{x0}, k_{y0} \)) | \( k_0, 0 \) |
| Harmonic potential frequency (\( \omega_x, \omega_y \)) | \( 5 \times 10^{11} \, \text{rad/s} \) |
| Potential energy (\( V(x, y) \)) | \( 0.5 m_e \left(\omega_x^2 x^2 + \omega_y^2 y^2 \right) \) |
Table I: Values of the simulation parameters.
Simulation Code in Python
import numpy as np
import matplotlib.pyplot as plt
# 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
# Harmonic potential parameters
omega_x = 5e11 # Oscillation frequency in x-direction (rad/s)
omega_y = 5e11 # Oscillation frequency in y-direction (rad/s)
# Define harmonic potential
V = 0.5 * m_e * (omega_x ** 2 * X ** 2 + omega_y ** 2 * Y ** 2)
# 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))
# 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()
for t in range(Nt):
# (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))
# Visualization of wave packet evolution
if t % 50 == 0:
plt.figure(figsize=(6, 5))
Rho = np.abs(psi) ** 2
plt.imshow(Rho / np.max(Rho) + V / np.max(V), extent=(-Lx / 2, Lx / 2, -Ly / 2, Ly / 2), cmap='hot')
plt.colorbar()
plt.title(f'Time step: {t}')
plt.xlabel(r'$x$ (m)', fontsize=14)
plt.ylabel(r'$y$ (m)', fontsize=14)
plt.tick_params(which="major", axis="both", direction="in", top=True, right=True, length=5, width=1, labelsize=12)
plt.pause(1) # Pause for 1 second
plt.close()
plt.show()
Simulation Results
The simulation illustrates the evolution of a quantum wave packet in a two-dimensional harmonic potential, described by the TDSE. The initial wave packet, \( \psi(x, y, t=0) \), is a localized Gaussian with specified widths (\( \sigma_x, \sigma_y \)) and initial momentum (\( k_{x0}, k_{y0} \)). Under the harmonic potential, the wave packet undergoes periodic oscillations. Using the Split-Step Fourier Method, the simulation accurately captures kinetic and potential energy contributions, revealing quantum features such as spreading, interference, and oscillation. The evolving probability density \( |\psi(x, y, t)|^2 \) provides insights into confined quantum systems, relevant to quantum mechanics and optics.
Animation in Python
The following code generates an animation of the wave function evolution for a particle in a harmonic potential.
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
# Harmonic potential parameters
omega_x = 5e11 # Oscillation frequency in x-direction (rad/s)
omega_y = 5e11 # Oscillation frequency in y-direction (rad/s)
# Define harmonic potential
V = 0.5 * m_e * (omega_x ** 2 * X ** 2 + omega_y ** 2 * Y ** 2)
# 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))
# 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=(8, 6))
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=12)
ax.set_ylabel(r'$y$ (m)', fontsize=12)
ax.tick_params(which="major", axis="both", direction="in", top=True, right=True, length=5, width=1, labelsize=10)
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=12)
return [im]
ani = animation.FuncAnimation(fig, update, frames=Nt, repeat=False)
# Save the animation
ani.save('electron_harmonic_simulation.mp4', writer='ffmpeg', fps=20)
plt.show()
Summary
In this article, we simulated the motion of a quantum particle in a two-dimensional harmonic potential using the Time-Dependent Schrödinger Equation (TDSE) and the Split-Step Fourier Method (SSFM).
The simulation shows how a Gaussian wave packet behaves under harmonic confinement, with periodic oscillations and quantum effects like spreading and interference. The SSFM accurately calculates the effects of kinetic and potential energy in both real and momentum spaces. Key parameters include the particle’s de Broglie wavelength, harmonic potential frequency, and grid size.
These results are useful for understanding systems like trapped ions, optical lattices, quantum computing, and precision measurements.
References
- Max. (2024). Simulation of Matter-Wave Diffraction Through a Grating. MatterWaveX. Retrieved from https://matterwavex.com/simulation-of-matter-wave-diffraction-through-a-multi-slit-grating/
- Feit, M. D., Fleck, J. A., & Steiger, A. (1982). Solution of the Schrödinger equation by a spectral method. Journal of Computational Physics, 47(3), 412–433. https://doi.org/10.1016/0021-9991(82)90091-2
- Griffiths, D. J., & Schroeter, D. F. (2018). Introduction to Quantum Mechanics (3rd ed.). Cambridge University Press.
- Wineland, D. J., & Blatt, R. (2008). Trapped ions in quantum information processing and quantum computation. Nature, 453(7198), 1008–1015. https://doi.org/10.1038/nature07124
- Tannor, D. J. (2007). Introduction to Quantum Mechanics: A Time-Dependent Perspective. University Science Books.
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.
