Path-Integral Simulation of Matter Wave Diffraction: Step-by-Step Guide

By /

Last Updated on December 22, 2025 by Dr. Sushanta Barman

Matter wave diffraction is a direct consequence of the wave nature of particles such as electrons, atoms, and ions. When the de Broglie wavelength (\(\lambda_{dB} = h/p,\) with \(h\) being the Planck constant and \(p\) being the momentum) of a particle becomes comparable to the dimensions of a slit or grating, interference effects emerge that cannot be explained using classical particle trajectories.

Diffraction is commonly analyzed using the Schrödinger equation or wave-optics analogies. While these approaches are rigorous, they often hide the physical role of phase accumulation during propagation. The path integral formulation provides an alternative and intuitive description by explicitly expressing diffraction as a coherent sum over all allowed propagation paths.

In this tutorial, we present a simple and practical implementation of matter wave diffraction using the path integral (Huygens–Fresnel) approach. The focus is on physical understanding and numerical implementation, rather than formal derivations. The method is demonstrated using a one-dimensional simulation of electron diffraction from a multi-slit grating. It is suitable for undergraduate students learning quantum wave propagation and interference.

Physical System and Geometry

We consider one-dimensional transverse diffraction of an electron matter wave with a de Broglie wavelength \(\lambda_{dB} = 50 \mathrm{nm}\), as shown in Fig. 1. The electron propagates along the \(x\)-direction, while diffraction occurs in the transverse \(y\)-direction.

A thin multi-slit transmission grating is placed in the plane \(x=0\), and the diffracted wave is observed on a detection screen located at a distance \(L_{\mathrm{screen}}\) downstream, as illustrated in Fig. 1.

In this model, only the transverse component of the electron wave function, \(\psi(y)\), is simulated, while longitudinal motion is treated implicitly through free-space propagation.

Figure 1. Schematic of the physical system and geometry used in the simulation. An incident electron matter wave propagates along the longitudinal \(x\)-direction and encounters a thin multi-slit grating placed at \(x=0\). Each open slit acts as a secondary source of probability amplitude, in accordance with the Huygens–Fresnel principle. The diffracted matter wave propagates to a detection screen located at a distance \(L_{\mathrm{screen}}\), where constructive and destructive interference lead to intensity maxima and minima in the transverse \(y\)-direction.

Incident Wave Function

At the grating plane (\(x=0\)), the incident electron beam is modeled by a Gaussian transverse profile,
\[\psi_{\text{in}}(y)= \exp\left[-\frac{(y-y_0)^2}{2\sigma_y^2}\right].\]

This represents a collimated beam centered at \(y_0\) with width \(\sigma_y\). The wave function is normalized numerically so that
\[\int |\psi_{\text{in}}(y)|^2\,dy = 1.\]

The transverse intensity profile \(|\psi_{in}(y)|^2\) of the incident electron wave packet at the grating plane is shown in Fig. 2. The wave packet is centered at \(y_0 = 0\) and has a width \(\sigma_y = 1\,\mu\text{m}\).

Figure 2. Incident transverse intensity profile of the electron beam at the grating plane. The Gaussian-shaped curve represents the normalized transverse probability density \( |\psi_{\mathrm{in}}(y)|^2 \), centered at \(y=0\), corresponding to a well-collimated incoming matter wave before interaction with the grating.

Grating Transmission Function (T(y))

The grating is modeled as a thin binary mask using a transmission function \(T(y)\).

\[T(y)=\begin{cases}1, & \text{inside a slit opening},\\[4pt] 0, & \text{inside an opaque bar}.\end{cases}\]

For a grating with slit width \(w\) (opening part), barrier width \(b\) (opaque part), and period \(d=w+b\), the slits are placed periodically along \(y\), as shown in Fig. 1. The wave function immediately after the grating is,

\[\psi_{\text{grating}}(y) = \psi_{\text{in}}(y)\,T(y).\]

This step enforces the physical constraint that only paths passing through the open slits contribute. These open slits act as secondary sources, as shown in Fig. 1. All paths blocked by opaque regions are removed at this stage.

The transmitted transverse intensity profile \(|\psi_{\mathrm{grating}}(y)|^2\) immediately after the grating is shown in Fig. 3.

Figure 3. Transmitted transverse intensity immediately after the grating. The profile \(|\psi_{\mathrm{grating}}(y)|^2\) shows the incident Gaussian beam after passing through the multi-slit grating. Intensity appears only at the slit positions. The opaque regions block the beam completely.

Path-Integral Propagation to the Screen

Here, the path-integral formulation is implemented in its scalar propagator form, which is equivalent to the Huygens–Fresnel principle in the paraxial, non-relativistic limit.

Each point within the slit openings acts as a secondary source of probability amplitude. The wave function at the screen is obtained by summing contributions from all allowed paths:

\[\psi_{\text{screen}}(y_s)= \int_{-\infty}^{\infty} dy \; K(y_s,y)\,\psi_{\text{grating}}(y),……(1)\]

where \(y_s\) denotes the transverse coordinate on the detection screen (see Fig. 1), and \(K(y_s,y)\) is the free-space propagation kernel that accounts for phase accumulation and geometric spreading during propagation from the grating plane to the screen.

For free-space propagation, the kernel is chosen as
\[K(y_s,y)= \frac{e^{ik_0 R}}{R},\qquad R=\sqrt{(y_s-y)^2 + L_{\text{screen}}^2}……(2)\]

  • The phase factor \(e^{ik_0R}\) represents phase accumulation along each path.
  • The factor \(1/R\) accounts for geometric spreading.

Overall multiplicative constants in the propagator are omitted, since the wave function is normalized after propagation and only relative intensities are of interest.

This expression is the quantum-mechanical analogue of the Huygens–Fresnel principle.

Numerical Implementation

The integral in Equation (1) is evaluated numerically on a discrete grid as,
\[\psi_{\text{screen}}(y_s) \approx \sum_j \frac{e^{ik_0 R(y_s,y_j)}}{R(y_s,y_j)} \psi_{\text{grating}}(y_j)\,dy. \]

After propagation, the wave function is normalized as,
\[ \int |\psi_{\text{screen}}(y)|^2\,dy = 1.\]

The measurable diffraction pattern is given by the intensity,
\[ I(y) = |\psi_{\text{screen}}(y)|^2.\]

Simulation Parameters

All simulation parameters used in the path-integral simulation are listed in Table I.

Parameter Value
Particle Electron
de Broglie wavelength (\(\lambda_{dB}\)) 50 nm
Wave number (\(k_0 = 2\pi/\lambda_{dB}\)) \(1.26\times10^{8}\,\mathrm{m^{-1}}\)
Transverse domain size (\(L_y\)) 100 μm
Number of transverse grid points (\(N_y\)) 20 000
Transverse grid spacing (\(\Delta y\)) 5 nm
Incident beam width (\(\sigma_y\)) 1 μm
Beam center position (\(y_0\)) 0
Slit width (\(w\)) 500 nm
Barrier width (\(b\)) 500 nm
Grating period (\(d = w + b\)) 1 μm
Number of slits (\(N\)) 100
Screen distance (\(L_{\mathrm{screen}}\)) 0.1 mm
Propagation method Path-integral (Huygens–Fresnel)
Propagation kernel \(K = e^{ik_0 R}/R\)

Table I: Values of the simulation parameters used in the path-integral simulation of electron matter wave diffraction. The chosen grid resolution adequately resolves the slit geometry and ensures stable numerical evaluation of the diffraction integral.

Python Code

Wave Packet Simulation
📄 Show Code 

import numpy as np
import matplotlib.pyplot as plt

# ------------------------------------------------------------
# 1. Physical constants and beam parameters
# ------------------------------------------------------------
hbar = 1.0545718e-34       # Reduced Planck constant (J·s)
m_e  = 9.10938356e-31      # Electron mass (kg)

# de Broglie wavelength and corresponding wave number k0 = 2π/λ
wavelength = 50e-9         # meters
k0 = 2.0 * np.pi / wavelength

# ------------------------------------------------------------
# 2. Transverse grid (y-direction)
# ------------------------------------------------------------
Ly = 1e-4                  # Total transverse simulation window (m)
Ny = 20000                 # Number of grid points along y
dy = Ly / Ny               # Grid spacing (m)
y  = np.linspace(-Ly/2, Ly/2, Ny)

# ------------------------------------------------------------
# 3. Incoming transverse wave at the grating plane (x=0)
# ------------------------------------------------------------
sigma_y = 1e-6           # 1σ width of the incident Gaussian envelope (m)
y0 = 0.0                   # Beam center position in y (m)

psi_in = np.exp(-(y - y0)**2 / (2.0 * sigma_y**2)).astype(np.complex128)

# Normalize the incident transverse wavefunction: ∫|ψ|^2 dy = 1
norm_in = np.sqrt(np.sum(np.abs(psi_in)**2) * dy)
psi_in /= norm_in

# ------------------------------------------------------------
# 4. Multi-slit grating transmission function T(y)
#    Ideal binary mask: slit = 1, opaque bar = 0
# ------------------------------------------------------------
slit_width       = 500e-9   # Slit opening width (m)
barrier_width    = 500e-9   # Opaque bar width between slits (m)
num_slits        = 100       # Total number of slits

period = slit_width + barrier_width  # Grating period (m)

T = np.zeros_like(y, dtype=float)

# Place 'num_slits' slits centered around y = 0
start_index = -num_slits // 2
end_index   =  num_slits // 2

for n in range(start_index, end_index):
    y_center = n * period
    y_min = y_center - slit_width / 2.0
    y_max = y_center + slit_width / 2.0
    T[(y >= y_min) & (y <= y_max)] = 1.0

# Wavefunction immediately after the grating (aperture applied)
psi_grating = psi_in * T

# ------------------------------------------------------------
# 5. Path-integral (Huygens–Fresnel) propagation to screen
# ------------------------------------------------------------
# Propagation distance from grating plane to observation (screen) plane
L_screen = 1.0e-4          # meters (0.1 mm)

psi_screen = np.zeros_like(y, dtype=np.complex128)

for i, ys in enumerate(y):
    # Point-to-point distance from each grating coordinate y to screen coordinate ys
    R = np.sqrt((ys - y)**2 + L_screen**2)
    phase = k0 * R
    # Free-space Green's function (scalar spherical wave, up to an overall constant)
    G = np.exp(1j * phase) / R
    psi_screen[i] = np.sum(G * psi_grating) * dy

# Normalize the screen wavefunction: ∫|ψ|^2 dy = 1 (within the simulated window)
norm_scr = np.sqrt(np.sum(np.abs(psi_screen)**2) * dy)
psi_screen /= norm_scr

# Intensity on the screen, normalized to its peak value
I_screen = np.abs(psi_screen)**2
I_screen /= np.max(I_screen)

# ------------------------------------------------------------
# 6. Plots: incident transverse profile + transmitted profile + diffraction pattern
# ------------------------------------------------------------

# (A) Incident transverse intensity profile |psi_in|^2
I_in = np.abs(psi_in)**2
I_in /= np.max(I_in)

plt.figure(figsize=(7, 4))
plt.plot(y * 1e6, I_in, linewidth=1.5)
plt.xlabel(r'Transverse coordinate $y$ ($\mu$m)', fontsize=14)
plt.ylabel(r'Normalized intensity $|\psi_{\rm in}(y)|^2$', fontsize=14)
plt.title('Incident transverse intensity profile', fontsize=12)
plt.grid(True, alpha=0.3)
plt.tight_layout()
plt.show()

# (B) Transmitted intensity profile right after the grating |psi_grating|^2
I_gr = np.abs(psi_grating)**2
if np.max(I_gr) > 0:
    I_gr /= np.max(I_gr)

plt.figure(figsize=(7, 4))
plt.plot(y * 1e6, I_gr, linewidth=1.5)
plt.xlabel(r'Transverse coordinate $y$ ($\mu$m)', fontsize=14)
plt.ylabel(r'Normalized intensity $|\psi_{\rm grating}(y)|^2$', fontsize=14)
plt.title('Transmitted transverse intensity (just after grating)', fontsize=12)
plt.grid(True, alpha=0.3)
plt.xlim(-20,20)
plt.tight_layout()
plt.show()

# (C) Diffraction pattern at the screen
plt.figure(figsize=(7, 4))
plt.plot(y * 1e6, I_screen, linewidth=1.5)
plt.xlabel(r'Transverse coordinate $y$ ($\mu$m)', fontsize=14)
plt.ylabel(r'Normalized intensity $I(y)$', fontsize=14)
plt.title('Electron diffraction pattern at the screen\n'
          rf'$\lambda = {wavelength*1e9:.0f}\,\mathrm{{nm}},\ N_{{\rm slits}}={num_slits}$',
          fontsize=12)
plt.grid(True, alpha=0.3)
plt.tight_layout()
plt.show()

Explanation of the Simulation Code

Discretization of the Transverse Coordinate

The transverse coordinate \(y\) is discretized on a uniform grid of length \(L_y\) with \(N_y\) points. The grid spacing \(\Delta y = L_y/N_y\) is chosen sufficiently small to resolve the slit width and ensure accurate representation of the aperture edges. All integrals over \(y\) are evaluated numerically using Riemann summation.

Construction of the Incident Wave Function

The incident electron at the grating plane is modeled as a Gaussian transverse wave packet. The wave function is explicitly normalized so that the total probability \(\int |\psi(y)|^2 dy = 1\). This normalization ensures that the intensity remains physically meaningful throughout propagation.

Implementation of the Grating via a Transmission Function

The multi-slit grating is implemented as a thin binary transmission mask \(T(y)\). The mask is constructed by setting \(T(y)=1\) inside slit openings and \(T(y)=0\) in opaque regions. The effect of the grating is applied by pointwise multiplication of the incident wave function with \(T(y)\), which removes all contributions from blocked paths.

Path-Integral Propagation to the Screen

Propagation from the grating plane to the detection screen is performed using a free-space scalar Green’s function (Eq. (2)). For each screen coordinate \(y_s\), the contribution from every aperture point \(y\) is summed with a phase factor \(\exp(i k_0 R)\), where \(R\) is the geometric distance between the two points, as shown in Fig. 1. The \(1/R\) factor accounts for geometric spreading. This discrete summation is the numerical implementation of the path-integral over all allowed propagation paths.

Normalization and Intensity Calculation

After propagation, the screen wave function is renormalized to maintain unit total probability within the simulated window. The observable diffraction pattern is obtained from the squared modulus of the wave function, and the intensity is scaled by its maximum value for clear visualization.

Visualization

Three intensity profiles are plotted: the incident beam profile (Fig. 2), the transmitted profile immediately after the grating (Fig. 3), and the final diffraction pattern at the screen (Fig. 4). This sequence demonstrates the evolution of the wave function from beam preparation to diffraction and interference.

Figure 4. Electron diffraction pattern at the detection screen. The normalized intensity \(I(y)\) shows clear interference fringes produced by the multi-slit grating. The central peak is the brightest. The side peaks arise from constructive and destructive interference.

Results and Discussion

The simulation clearly demonstrates matter wave diffraction from a multi-slit grating using the path-integral approach. The incident beam has a smooth Gaussian transverse profile, as shown in Fig. 2. After the grating, the wave function is spatially filtered, with intensity confined to the slit openings (see Fig. 3). This confirms the correct implementation of the transmission function \(T(y)\).

At the detection screen, a well-defined diffraction pattern is observed, as shown in Fig. 4. A strong central maximum appears, followed by weaker side peaks. These fringes arise from constructive and destructive interference between probability amplitudes originating from different slits.

The finite slit width produces an overall diffraction envelope, while the periodic grating leads to sharp interference maxima. The results agree with standard diffraction theory and validate the path-integral method as an intuitive and effective tool for simulating matter wave interference.

Summary

This article presented a step-by-step simulation of matter wave diffraction using the path-integral approach. Electron diffraction from a multi-slit grating was modeled in one transverse dimension.

The incident beam was described by a Gaussian wave function. The grating was implemented using a binary transmission function that allows only open paths to contribute. Free-space propagation to the screen was performed using a Huygens–Fresnel kernel.

The resulting diffraction pattern showed clear interference fringes with a strong central maximum and weaker side peaks. These features arise from phase differences accumulated along different paths. The simulation agrees with standard diffraction theory.

The approach provides a simple and intuitive way to understand matter wave diffraction. It is well-suited for undergraduate teaching and numerical exploration of matter wave phenomena.

References

  1. Sushanta Barman, “de Broglie Hypothesis: The Wave-Particle Duality of Matter,” MatterWaveX.Com, August 25 (2024).
  2. Sushanta Barman, “Schrödinger Equation for Matter Wave Dynamics,” MatterWaveX.Com, September 7 (2024).
  3. Konstantin K. Likharev, “The Huygens Principle,” Phys.LibreTexts.Org, Accessed on December 21 (2025).
  4. Sushanta Barman, “Why Physicists Love Gaussian Wave Packets,” MatterWaveX.Com, May 25 (2025).
  5. Eric R. Jones, Roger A. Bach, and Herman Batelaan, “Path integrals, matter waves, and the double slit,” European Journal of Physics 36, 065048 (2015).
  6. Dennis V. Perepelitsa, “Path Integrals in Quantum Mechanics,” MIT Department of Physics, Cambridge, Massachusetts, Accessed on December 21 (2025).

Leave a Comment

Your email address will not be published. Required fields are marked *

Related Posts

Scroll to Top