Angular Spectrum Propagation Drift: Too-Small Step Size, FFT/IFFT Error, and a Better Match to the Gaussian Beam

Learn why FFT-based angular spectrum propagation drifts at long range: tiny step sizes accumulate FFT/IFFT error. Larger dz restores Gaussian beam accuracy.

FFT-based angular spectrum propagation is a solid way to evolve a transverse optical field over distance, but it can bite when grids, edges, and step sizes collide. A common pattern is to add absorbing boundaries and split the propagation into small steps so the field never hits the box edges. If the propagated profile then drifts from the analytical solution only at long ranges, suspicion usually falls on the filter. In the case below, though, the filter is fine. The instability comes from an overly small propagation step that forces too many FFT/ifft cycles and accumulates numerical error.

Problem setup

The goal is to start from a Gaussian field sampled at z = z0, propagate it a distance d, and compare the numerical result at z = z0 + d with the analytical Gaussian beam at the same plane. A stepwise propagation with an absorbing boundary mask is used to keep energy away from the borders. For relatively short distances everything matches nicely. For longer distances, once the beam approaches the edges, the numerical result deviates from the analytical beam, even with a smooth mask in place.

Code example that reproduces the mismatch

The snippet below shows the core logic: a Gaussian beam definition, a single-step propagator in the spatial-frequency domain, a multi-step driver with an absorbing edge, and a test that compares the propagated field to the analytical one at the target plane. Variable names are in English and the logic matches the described setup.

import numpy as np
from numpy.fft import fft2, ifft2, fftfreq
from numpy import exp, pi, sqrt
arctg = np.arctan2
class BeamCore:
  def __init__(self, zR: float, lambda0: float):
    self.zR = zR
    self.lambda0 = lambda0
    self.k = 2 * pi / self.lambda0
    self.w0 = sqrt(self.lambda0 * self.zR / pi)
    self.theta0 = self.w0 / self.zR
    self.a0 = sqrt(2 / (pi * self.w0))
    self.i0 = self.a0 * self.a0.conjugate()
  def w(self, z: np.array):
    return self.w0 * sqrt(1 + (z / self.zR)**2)
  def R(self, z: np.array):
    z_safe = np.where(z == 0.0, 1e-20, z)
    rval = z_safe * (1 + (self.zR / z_safe)**2)
    rval = np.where(z == 0.0, np.inf, rval)
    return rval
  def gouy(self, z: np.array):
    return arctg(z, self.zR)
class GaussianBeam(BeamCore):
  def __init__(self, zR: float, lambda0: float):
    super().__init__(zR, lambda0)
  def field(self, z: np.array = 0.0, r: np.array = 0.0, phi: np.array = 0.0) -> np.array:
    z = np.asarray(z, dtype=float)
    w = self.w(z)
    R = self.R(z)
    zeta = self.gouy(z)
    U = self.a0 * (self.w0 / w) * exp(-r**2 / w**2 + (-self.k * z - self.k * (r**2) / (2 * R) + zeta) * 1.0j)
    return U
def step_propagate(U, dz, kz):
  UF = fft2(U)
  if dz > 0.0:
    UFprop = UF * exp(1j * kz * dz).conj()
  elif dz < 0.0:
    UFprop = UF * exp(1j * kz * (-dz))
  Uout = ifft2(UFprop)
  return Uout
def propagate_sliced(U, zdist, x1d, y1d, k0, edge_mask):
  nx, ny = U.shape
  dx = (x1d[-1] - x1d[0]) / (nx - 1)
  dy = (y1d[-1] - y1d[0]) / (ny - 1)
  fx = fftfreq(nx, d=dx)
  fy = fftfreq(ny, d=dy)
  FX, FY = np.meshgrid(fx, fy)
  kx = 2 * pi * FX
  ky = 2 * pi * FY
  kz = sqrt(k0**2 - kx**2 - ky**2 + 0j)
  Ustep = U
  while zdist != 0.0:
    U_faded = Ustep * edge_mask
    dz = 0.00005  # very small step
    dz = min(dz, zdist)
    Ustep = step_propagate(U_faded, dz, kz)
    zdist -= dz
  return Ustep
# Parameters
lambda0 = 800 * (10**-9)
zR = pi * (10**(-3)) / 2
beam = GaussianBeam(zR, lambda0)
w0 = beam.w0
k0 = beam.k
# Absorbing edge mask built on a larger grid
pad = 500
grid = 1000 + 2 * pad
X = np.linspace(-8 * w0, 8 * w0, grid)
Y = np.linspace(-8 * w0, 8 * w0, grid)
XX, YY = np.meshgrid(X, Y)
r_mask = sqrt(XX**2 + YY**2)
def gauss_edge(r, pad_size):
  sigma = 1
  r_shift = (r - 4 * w0) / w0
  return exp(-(r_shift**2) / (2 * sigma**2))
edge_mask = np.where(r_mask > 4 * w0, gauss_edge(r_mask, pad), 1)
# Working grid for the field and propagation
x = np.linspace(-8 * w0, 8 * w0, 2000)
y = np.linspace(-8 * w0, 8 * w0, 2000)
Xg, Yg = np.meshgrid(x, y)
r = sqrt(Xg**2 + Yg**2)
phi = arctg(Yg, Xg)
# Initial field at z = zR
U0 = beam.field(zR, r, phi)
# Propagate by 5 * zR to reach z = 6 * zR
U_num = propagate_sliced(U0, 5 * zR, x, y, k0, edge_mask)
# Analytical field at z = 6 * zR
U_ref = beam.field(6 * zR, r, phi)
# Simple pixel-distance measure
def pix_distance(a: np.array, b: np.array):
  v = abs((a - b).flatten())
  return np.dot(v, v) / len(v)
err = pix_distance(U_num, U_ref)
print(err)

With the tiny step above, the deviation at long range is large. The field starts clean, the edges are faded, yet the final profile does not align with the analytical Gaussian as it should.

What actually goes wrong

The absorbing boundary layer is not the culprit. The issue is the propagation step being too small. Each step applies an FFT and an inverse FFT to the field. When the number of steps is pushed unnecessarily high, the repeated transforms accumulate numerical error. That error builds up across the whole propagation sequence and dominates once the beam has expanded.

Fix: increase the propagation step size

By increasing the step size, you reduce the number of FFT/ifft pairs and the accumulated error drops. In this setup, using dz = 0.001 yields a markedly better match to the analytical field at the target plane.

def propagate_sliced(U, zdist, x1d, y1d, k0, edge_mask):
  nx, ny = U.shape
  dx = (x1d[-1] - x1d[0]) / (nx - 1)
  dy = (y1d[-1] - y1d[0]) / (ny - 1)
  fx = fftfreq(nx, d=dx)
  fy = fftfreq(ny, d=dy)
  FX, FY = np.meshgrid(fx, fy)
  kx = 2 * pi * FX
  ky = 2 * pi * FY
  kz = sqrt(k0**2 - kx**2 - ky**2 + 0j)
  Ustep = U
  while zdist != 0.0:
    U_faded = Ustep * edge_mask
    dz = 0.001  # larger step size improves results here
    dz = min(dz, zdist)
    Ustep = step_propagate(U_faded, dz, kz)
    zdist -= dz
  return Ustep

You can quantify the improvement with the same pixel-distance function. With dz = 0.001 the distance is about 2.35, whereas with dz = 0.00005 it jumps to roughly 64.49. You can adjust dz and observe how the distance behaves for your target range.

Why this matters

In angular spectrum propagation, step size is not just a convenience parameter. It controls the number of spectral transforms you perform. If the step is too small, the method spends accuracy on redundant updates. Even when you use a smooth absorbing mask and a grid large enough to contain the beam, an overly fine step can dominate the error budget simply by forcing too many FFT/ifft cycles.

Conclusion

When a stepwise FFT propagator deviates from an analytical Gaussian at long ranges, do not blame the absorbing boundary first. Check the propagation step. In the setup shown here, increasing dz from 0.00005 to 0.001 significantly reduces the mismatch, as confirmed by a simple pixel-distance metric. Start with a larger step, verify against the analytical profile, and only then refine as needed for your propagation distance and grid.