Angular Spectrum Method: Fixing Backward Fourier Propagation with a Safe Phase Conjugate for a Gaussian Beam

Learn why your Fourier-based angular spectrum code shows a shrinking Gaussian beam and how using a conjugated exp(i*kz*dz) corrects direction cleanly. See how.

Fourier-based beam propagation is a compact way to advance complex fields in space, but sign conventions and phase factors can easily send you in the wrong direction. A common symptom is a Gaussian profile that appears to shrink when you expect it to expand. Below is a minimal walk-through of the issue and a practical fix that stays within the original numerical setup.

Problem setup

The goal is to propagate a Gaussian beam a distance dz using the angular spectrum method. The computation runs, but the beam appears to propagate backwards: the profile shrinks around dz = z0, as if the field were moving toward its waist. Switching fft/ifft did not change the outcome. Flipping the sign in the complex phase exp(±1j · kz · dz) produced overflow warnings and invalid FFT values.

import numpy as np
from numpy.fft import fft2, ifft2, fftfreq
from matplotlib import pyplot as plt
EXP = np.exp
PI = np.pi
SQRT = np.sqrt
ATAN2 = np.arctan2
PHASE = np.angle
plt.rcParams["figure.autolayout"] = True
class RayField:
    def __init__(self, z_r: float, lam: float):
        self.z_r = z_r
        self.lam = lam
        self.k0 = 2 * PI / self.lam
        self.w0 = SQRT(self.lam * self.z_r / PI)
        self.theta0 = self.w0 / self.z_r
        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.z_r) ** 2)
    def R(self, z: np.array):
        zsafe = np.where(z == 0.0, 1e-20, z)
        curv = zsafe * (1 + (self.z_r / zsafe) ** 2)
        curv = np.where(z == 0.0, np.inf, curv)
        return curv
    def gouy(self, z: np.array):
        return ATAN2(z, self.z_r)
class GaussianMode(RayField):
    def __init__(self, z_r: float, lam: float):
        super().__init__(z_r, lam)
    def field(self, z=0.0, rho=0.0, phi=0.0):
        z = np.asarray(z, dtype=float)
        k = self.k0
        w0 = self.w0
        a0 = self.a0
        w = self.w(z)
        R = self.R(z)
        g = self.gouy(z)
        U = a0 * (w0 / w) * EXP(-rho ** 2 / w ** 2 + (-k * z - k * (rho ** 2) / (2 * R) + g) * 1.0j)
        return U
lam = 800 * (10 ** -9)
z_r = PI * (10 ** -3) / 2
step = 0.006
def propagate_spectrum(U, dz, npts, spacing):
    fx = fftfreq(npts, d=spacing)
    fy = fftfreq(npts, d=spacing)
    FX, FY = np.meshgrid(fx, fy)
    kx = 2 * PI * FX
    ky = 2 * PI * FY
    kz = SQRT(k ** 2 - kx ** 2 - ky ** 2 + 0j)
    spec = fft2(U)
    spec_out = spec * EXP(1j * kz * dz)  # observed to send the beam "backwards"
    U_out = ifft2(spec_out)
    return U_out
model = GaussianMode(z_r, lam)
N = 1000
k = model.k0
W0 = model.w0
xs = np.linspace(-4 * W0, 4 * W0, N)
ys = np.linspace(-4 * W0, 4 * W0, N)
dx = 8 * W0 / (N - 1)
X, Y = np.meshgrid(xs, ys)
RHO = SQRT(X ** 2 + Y ** 2)
PHI = ATAN2(Y, X)
U0 = model.field(z_r, RHO, PHI)
U1 = propagate_spectrum(U0, step, N, dx)
def show(u):
    mag = abs(u)
    ang = PHASE(u)
    plt.figure(figsize=(12, 5))
    plt.subplot(1, 2, 1)
    plt.pcolor(X, Y, mag, shading='auto', cmap='viridis')
    plt.title('Magnitude')
    plt.xlabel('x')
    plt.ylabel('y')
    plt.colorbar()
    plt.subplot(1, 2, 2)
    plt.pcolor(X, Y, ang, shading='auto', cmap='twilight')
    plt.title('Phase')
    plt.xlabel('x')
    plt.ylabel('y')
    plt.colorbar()
    plt.tight_layout()
    plt.show()
show(U1)

What’s going on

The direction of propagation in an angular-spectrum implementation is governed by the complex phase factor multiplying the spectrum. In this setup, simply replacing 1j by -1j in the exponential led to overflow and invalid FFT values. That behavior is consistent with a phase change that turns some terms into exponentially growing components, which numerically explodes. The practical observation is that switching fft and ifft did not alter the behavior, while a direct sign flip in the exponential produced runtime warnings. The result is a field that looks like it is advancing toward the waist instead of away from it.

Working fix

Conjugating the propagation factor flips the direction without triggering the overflow seen with a raw sign change. In other words, keep the kernel as exp(1j · kz · dz) and apply a complex conjugate to it before multiplying the spectrum. With that, the calculation runs cleanly and the beam propagates away from its waist.

def propagate_spectrum(U, dz, npts, spacing):
    fx = fftfreq(npts, d=spacing)
    fy = fftfreq(npts, d=spacing)
    FX, FY = np.meshgrid(fx, fy)
    kx = 2 * PI * FX
    ky = 2 * PI * FY
    kz = SQRT(k ** 2 - kx ** 2 - ky ** 2 + 0j)
    spec = fft2(U)
    phase = EXP(1j * kz * dz).conj()  # flips direction safely
    spec_out = spec * phase
    U_out = ifft2(spec_out)
    return U_out

This change is localized and does not require altering the rest of the pipeline. Using the conjugate avoids the numeric blow-up encountered when changing the sign inside the exponential.

Why this matters

With Fourier propagation, tiny mismatches in phase conventions accumulate into big visual differences. A wrong direction may look like a crisping of the profile at the wrong distance, leading you to incorrect conclusions about focusing behavior. Attempts to flip signs ad hoc can also destabilize the computation. Using the complex conjugate of the propagation factor offers a minimal, robust correction within the same formulation.

Takeaways

If a Fourier propagation seems to run backward, focus on the phase factor that multiplies the angular spectrum. In this setup, multiplying by exp(1j · kz · dz).conj() corrects the direction and avoids the overflow observed when changing 1j to -1j. This keeps the numerical method stable and the physical picture consistent with a beam moving away from its waist.