In this article, we will explore the simulation of double-slit diffraction for matter waves (electrons) using the Split-Step Fourier Method (SSFM).

This method is widely used for solving the time-dependent SchrÃ¶dinger equation, especially in cases involving complex potential landscapes.

**Time-Dependent SchrÃ¶dinger Equation**

The time evolution of a quantum wave packet is governed by the time-dependent SchrÃ¶dinger equation (TDSE):

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

where \( \psi(\mathbf{r}, t) \) is the wave function, \( \hbar \) is the reduced Planck constant, \( m \) is the mass of the particle, and \( V(\mathbf{r}) \) is the potential energy.

**Split-Step Fourier Method (SSFM)**

The Split-Step Fourier Method allows us to separate the effects of the kinetic and potential energy operators by alternating between real and momentum space. The TDSE can be approximated using the Trotter-Suzuki decomposition [1]:

$$\psi(\mathbf{r}, t + \Delta t) \approx e^{-i\hat{H} \Delta t / \hbar} \psi(\mathbf{r}, t) \approx e^{-i \hat{V} \Delta t / 2\hbar} e^{-i \hat{T} \Delta t / \hbar} e^{-i \hat{V} \Delta t / 2\hbar} \psi(\mathbf{r}, t),$$

where \( \hat{T} \) is the kinetic energy operator and \( \hat{V} \) is the potential energy operator. More details on the split-step Fourier method can be found here [2]

**System Under Investigation**

In this simulation, we investigate the diffraction of an electron wave packet in two dimensions \((x,y)\) through a double-slit, as shown in Figure 1. This classic demonstration showcases the fundamental principles of quantum mechanics. This system highlights the wave-particle duality of electrons, showing how they exhibit interference patterns similar to light waves when passing through slits.

**Electron Wave Packet**

The electron wave packet is initialized as a Gaussian function in both the \( x \) and \( y \) directions. The initial parameters for the wave packet are as follows:

**Width of the Gaussian packet \((\sigma_x, \sigma_y)\)**: \(50 \, \text{nm}\)**Initial position \((x_0, y_0)\)**: \((-200 \, \text{nm}, 0 \, \text{nm})\)**Initial wave vector \((k_{x0}, k_{y0})\)**: Corresponding to a de Broglie wavelength of \(20 \, \text{nm}\)

The wave function at \(t = 0\) is given by:

$$\psi_0(x, y, t=0) = \exp\left[-\left(\frac{(x – x_0)^2}{2\sigma_x^2} + \frac{(y – y_0)^2}{2\sigma_y^2}\right)\right] \exp(i (k_{x0} x + k_{y0} y)),$$

where the Gaussian terms represent the spatial localization of the wave packet, and the exponential term with the imaginary unit represents the initial momentum of the packet.

**Double Slit Potential**

The potential energy landscape in the simulation domain is designed to represent a double slit. The parameters of the double slit are:

**Slit width**: \(75 \, \text{nm}\)**Gap between slits**: \(50 \, \text{nm}\)**Thickness of the double-slit grating**: \(30 \, \text{nm}\)**Potential barrier height**: \(1 \times 10^{10} \, \text{J}\) (to effectively act as an impenetrable barrier for the electron wave packet)**Geometric potential barrier**: The double slit is centered in the simulation domain. The geometric potential energy \( V(x, y) \) is defined such that it is zero everywhere except in the regions defining the slits, where it takes a very high value (acting as barriers):

$$V(x, y) = \begin{cases}1 \times 10^{10} \, \text{J} & \text{if } (x, y) \text{ is within the obstruction region} \\0 & \text{otherwise}\end{cases}$$

**Discretization of Equations**

To implement this method numerically, we discretize the simulation domain and time evolution:

1.** Grid Parameters**:

- Simulation domain size: \( L_x = L_y = 1 \, \mu m \)
- Number of grid points: \( N_x = N_y = 512 \)
- Grid spacing: \( \Delta x = \frac{L_x}{N_x} \), \( \Delta y = \frac{L_y}{N_y} \)

2. **Fourier Space Components**:

- Wave numbers: $$ k_x = \frac{2\pi}{L_x} \text{FFTfreq}(N_x, \Delta x), \quad k_y = \frac{2\pi}{L_y} \text{FFTfreq}(N_y, \Delta y), $$ where
**FFTfreq**stands for Fast Fourier Transform frequencies. - Fourier space grid:

$$K_x, K_y = \text{meshgrid}(k_x, k_y), \quad K^2 = K_x^2 + K_y^2$$

3. **Time Evolution**:

- Time step: \( \Delta t = 3 \times 10^{-14} \, s \)
- Number of time steps: \( N_t = 501 \)
- Evolution using SSFM:
- Half step in real space:

$$\psi(\mathbf{r}, t + \Delta t/2) = \psi(\mathbf{r}, t) \exp\left(-\frac{i V(\mathbf{r}) \Delta t}{2 \hbar}\right)$$ - Fourier transform to momentum space:

$$\tilde{\psi}(\mathbf{k}, t + \Delta t/2) = \mathcal{F}[\psi(\mathbf{r}, t + \Delta t/2)]$$ - Full step in momentum space:

$$\tilde{\psi}(\mathbf{k}, t + \Delta t) = \tilde{\psi}(\mathbf{k}, t + \Delta t/2) \exp\left(-\frac{i \hbar K^2 \Delta t}{2m}\right)$$ - Inverse Fourier transform to return to real space:

$$\psi(\mathbf{r}, t + \Delta t) = \mathcal{F}^{-1}[\tilde{\psi}(\mathbf{k}, t + \Delta t)]$$ - Another half-step in real space:

$$\psi(\mathbf{r}, t + \Delta t) = \psi(\mathbf{r}, t + \Delta t/2) \exp\left(-\frac{i V(\mathbf{r}) \Delta t}{2 \hbar}\right)$$

- Half step in real space:

**Visualization**

To visualize the wave packet evolution, we normalize and plot the wave function at regular intervals. The combined plot of the normalized wave function and potential energy helps in observing the diffraction pattern.

By following the above steps, we can simulate and visualize the diffraction of a wave packet through a double slit, illustrating the wave nature of particles and the fundamental principles of quantum mechanics. A simulation code in Python is provided below.

**Numerical 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 = 50e-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
# Initial Gaussian wave packet (ensure it's complex)
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))
# Double slit potential
slit_width = 75e-9 # Slit width (meters)
slit_gap = 50e-9 # Gap between slits (meters)
slit_thickness = 30e-9 # Slit thickness (meters)
V = np.zeros((Ny, Nx))
# Representing the two slits
for i in range(Nx):
for j in range(Ny):
if (Nx // 2 - slit_thickness / (2 * dx) <= i <= Nx // 2 + slit_thickness / (2 * dx)) and \
(j < Ny // 2 - (slit_gap + slit_width) / (2 * dy) or j > Ny // 2 + (slit_gap + slit_width) / (2 * dy) or \
(Ny // 2 - slit_gap / (2 * dy) <= j <= Ny // 2 + slit_gap / (2 * dy))):
V[j, i] = 1e10
# 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()
# Normalize the wave function
psi /= np.sqrt(np.sum(np.abs(psi) ** 2))
for t in range(Nt):
# Half step in real space with potential energy operator
psi *= np.exp(-1j * V * dt / (2 * hbar))
# Full step in momentum space with kinetic energy operator
psi_k = np.fft.fft2(psi)
psi_k *= np.exp(-1j * hbar * K2 * dt / (2 * m_e))
psi = np.fft.ifft2(psi_k)
# Another half step in real space with potential energy operator
psi *= np.exp(-1j * V * dt / (2 * hbar))
# Normalize the wave function
#psi /= np.sqrt(np.sum(np.abs(psi) ** 2) * dx * dy)
print (int(t))
# Visualization of wave packet evolution
if t % 50 == 0:
plt.figure(figsize=(6,5))
normalized_psi = np.abs(psi) ** 2 / np.max(np.abs(psi) ** 2)
normalized_V = V / np.max(V)
combined = normalized_psi + normalized_V
plt.imshow(combined, extent=(-Lx / 2, Lx / 2, -Ly / 2, Ly / 2), cmap='magma')
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 seconds
plt.close() # Close the plot
```

**Simulation Results**

The simulation of wave packet diffraction through a double-slit effectively demonstrates the wave nature of particles and quantum interference. A video animation of the time evolution of the electron wave packet is depicted in Figure 2. Initially, the wave packet is a localized Gaussian distribution centered at \((-200 \, \text{nm}, 0)\) with momentum directed towards the double slit. As the wave packet evolves, it interacts with the slits, producing a characteristic interference pattern.

**Conclusion**

This simulation provides a clear demonstration of the wave-particle duality and the impact of potential barriers on wave packets. The Split-Step Fourier Method offers a powerful technique for solving the TDSE and can be extended to more complex systems and potentials.

**References**

[1] McLachlan, Robert I., and G. Reinout W. Quispel. “Splitting methods.” *Acta Numerica *11 (2002): 341-434.

[2] Sushanta Barman, “Wave Packet Simulation Techniques in Matter Wave Optics“, *MatterWaveTech.Com*, February 8, (2024).