Last Updated on August 19, 2024 by Max
Electron matter-wave diffraction demonstrates the wave-particle duality of electrons by forming interference patterns, much like light diffraction. This simulation models electron diffraction through a multi-slit grating using the time-dependent Schrödinger equation and the split-step Fourier method.
The simulation provides valuable insights into quantum mechanics by capturing the wave packet’s evolution. The study includes detailed visualizations and analyses of the diffraction patterns, showcasing the effectiveness of numerical methods in exploring quantum phenomena, with practical applications in electron microscopy and quantum computing.
This article details the mathematical framework, numerical methods, and simulation results, making it a valuable resource for researchers and students.
Numerical Methods
Time-Dependent Schrödinger Equation
The Time-Dependent Schrödinger Equation (TDSE) in two dimensions is expressed as:
\[i\hbar \frac{\partial \psi(x, y, t)}{\partial t} = -\frac{\hbar^2}{2m_e} \left(\frac{\partial^2 \psi(x, y, t)}{\partial x^2} + \frac{\partial^2 \psi(x, y, t)}{\partial y^2}\right) + V(x, y) \psi(x, y, t)\]
where \( \psi(x, y, t) \) is the wave function, \( i \) is the imaginary unit, \( \hbar \) is the reduced Planck constant, \( m_e \) is the mass of the particle (electron), and \( V(x, y) \) is the potential energy function that influences the dynamics of the particle.
This equation describes the quantum evolution of a particle in a 2D potential. Numerical methods such as the Split-Step Fourier Method (SSFM) enable the study of the system’s time-dependent behavior.
Initial Wave Packet
The initial wave packet is defined as a Gaussian function with a specified width and initial position in space. The wave packet in two dimensions is given by,
\[\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] \exp\left[i (k_{x0} x + k_{y0} y)\right],\]
where \( \sigma_x \) and \( \sigma_y \) are the widths of the Gaussian packet, \( (x_0, y_0) \) is the initial position, and \( (k_{x0}, k_{y0}) \) is the initial wave vector. The distribution of the initial wave packet (probability density \(|\psi(x,y,t)|^2)\) is shown in Figure 1.
Normalization
The wave packet is normalized so that the total probability density equals 1:
\[\psi(x, y,t) \rightarrow \frac{\psi(x, y,t)}{\sqrt{\sum |\psi(x, y,t)|^2}}\]
Potential Energy Function (Multi-Slit Grating)
The potential energy \( V(x, y) \) representing the multi-slit grating is constructed as a series of slits and barriers:
\[ V(x, y) = \begin{cases} V_0 & \text{if } x_1 \leq x \leq x_2 \text{ and } y \text{ is in a barrier region} \\ 0 & \text{otherwise} \end{cases} \]
where \( V_0 \) is the height of the potential. The distribution of \( V(x, y) \) is shown in Figure 1.
Split-Step Fourier Method
The time-dependent Schrödinger equation is solved using the split-step Fourier method, which divides the evolution into small time steps and alternates between real space (where potential energy is applied) and Fourier space (where kinetic energy is applied). The spatial domain is discretized into a grid with points \( x_i \) and \( y_j \), and time is discretized into steps \( t_n \). The following steps outline the split-step Fourier method used to solve the Schrödinger equation.
- Half-step evolution in real space (Potential Energy):
The wave function is evolved in real space according to the potential energy:
\[\psi(x, y, t + \frac{\Delta t}{2}) = \psi(x, y, t) \cdot \exp\left(-\frac{i V(x, y) \Delta t}{2 \hbar}\right)\] - Fourier transform (to Momentum Space):
The wave function is transformed to momentum space using the Fast Fourier Transform (FFT):
\[\tilde{\psi}(k_x, k_y, t) = \mathcal{F}{\psi(x, y, t)}\] - Full-step evolution in momentum space (Kinetic Energy):
The wave function is evolved in momentum space according to the kinetic energy:
\[\tilde{\psi}(k_x, k_y, t + \Delta t) = \tilde{\psi}(k_x, k_y, t) \cdot \exp\left(-\frac{i \hbar (k_x^2 + k_y^2) \Delta t}{2 m_e}\right)\] - Inverse Fourier transform (back to Real Space):
The wave function is transformed back to real space:
\[\psi(x, y, t + \Delta t) = \mathcal{F}^{-1}{\tilde{\psi}(k_x, k_y, t + \Delta t)}\] - Final half-step evolution in real space (Potential Energy):
The wave function undergoes a final half-step evolution in real space:
\[\psi(x, y, t + \Delta t) = \psi(x, y, t + \frac{\Delta t}{2}) \cdot \exp\left(-\frac{i V(x, y) \Delta t}{2 \hbar}\right)\]
All simulation parameters related to the initial wave packet, grating potential barrier, and spatial and temporal discretization are listed in Table I.
| Parameter | Value |
|---|---|
| Reduced Planck constant (ℏ) | 1.0545718e-34 J·s |
| Mass of electron (mₑ) | 9.10938356e-31 kg |
| De Broglie wavelength (λ) | 20e-9 m |
| Wave number (k₀) | 2π / λ |
| Simulation domain size (Lₓ, Lᵧ) | 1e-6 m |
| Number of grid points (Nₓ, Nᵧ) | 512 |
| Grid spacing (Δx, Δy) | Lₓ / Nₓ, Lᵧ / Nᵧ |
| Time step (Δt) | 3e-14 s |
| Number of time steps (Nₜ) | 501 |
| Gaussian packet width (σₓ, σᵧ) | 60e-9 m |
| Initial packet position (x₀, y₀) | -2e-7 m, 0 m |
| Initial wave vector (kₓ₀, kᵧ₀) | k₀, 0 |
| Slit width (a) | 3e-8 m |
| Barrier width (b) | 3e-8 m |
| Grating thickness (d) | 2e-8 m |
| Potential height (V₀) | 5 * ℏ² * k₀² / (2 * mₑ) |
| Number of slits (Nₛ) | 40 |
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
# Grating parameters
slit_width = 3e-8 # Width of each slit (meters)
barrier_width = 3e-8 # Width of the barrier between slits (meters)
grating_thickness = 2e-8 # Thickness of the grating (meters)
h0 = 5 * hbar**2 * k0**2 / (2 * m_e) # Height of the potential for the barrier
num_slits = 40 # Number of slits
# 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 (Multi-slit grating)
V = np.zeros_like(X)
x1, x2 = -grating_thickness / 2, grating_thickness / 2
y_coords = np.linspace(-0.45 * num_slits * barrier_width - 0.5 * num_slits * slit_width,
0.45 * num_slits * barrier_width + 0.5 * num_slits * slit_width, 2 * num_slits)
for i in range(Nx):
for j in range(Ny):
if x1 <= X[i, j] <= x2:
for k in range(0, len(y_coords), 2):
if y_coords[k] <= Y[i, j] <= y_coords[k + 1]:
V[i, j] = h0
# 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()
Visualization
At selected time steps, the probability density \( |\psi(x, y, t)|^2 \) is visualized along with the potential energy \( V(x, y) \) to observe the diffraction pattern.
These steps iteratively evolve the wave function, allowing the simulation of the diffraction of electron matter-waves by a multi-slit grating.
Discussion of Simulation Results
- Diffraction Pattern: The simulation successfully reproduces the electron diffraction pattern, as shown in Figure 2, with clear interference fringes analogous to those observed in optical diffraction.
- Wave-Packet Evolution: The wave packet evolves as expected, with distinct diffraction occurring at the multi-slit grating.
- Parameter Sensitivity: The results are sensitive to the slit width, barrier width, and number of slits, with variations in these parameters affecting the fringe spacing and intensity. These effects can be explored by modifying the corresponding parameters in the provided Python code.
- Practical Applications: The interference patterns relate to electron microscopy, where they enhance imaging resolution, and quantum computing, where they underpin quantum coherence and superposition necessary for qubit operations.
Animation of the Diffraction Pattern
The video below presents an animation of the time evolution of the diffraction pattern in real space \((x, y)\).
Animation Code in Python
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
# Grating parameters
slit_width = 3e-8 # Width of each slit (meters)
barrier_width = 3e-8 # Width of the barrier between slits (meters)
grating_thickness = 2e-8 # Thickness of the grating (meters)
h0 = 5 * hbar**2 * k0**2 / (2 * m_e) # Height of the potential for the barrier
num_slits = 40 # Number of slits
# 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 (Multi-slit grating)
V = np.zeros_like(X)
x1, x2 = -grating_thickness / 2, grating_thickness / 2
y_coords = np.linspace(-0.45 * num_slits * barrier_width - 0.5 * num_slits * slit_width,
0.45 * num_slits * barrier_width + 0.5 * num_slits * slit_width, 2 * num_slits)
for i in range(Nx):
for j in range(Ny):
if x1 <= X[i, j] <= x2:
for k in range(0, len(y_coords), 2):
if y_coords[k] <= Y[i, j] <= y_coords[k + 1]:
V[i, j] = h0
# 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('electron_diffraction_simulation.mp4', writer='ffmpeg', fps=20)
plt.show()
Conclusion
The simulation of electron matter-wave diffraction through a multi-slit grating clearly illustrates wave-particle duality, generating distinct interference patterns. Using the split-step Fourier method, this study emphasizes the value of numerical techniques in exploring quantum phenomena. The findings have practical applications in electron microscopy and quantum computing, making this work a valuable resource for understanding quantum mechanics.
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.
