Angular spectrum method
The angular spectrum method is a fundamental technique in scalar diffraction optics for propagating a monochromatic scalar wave field from an initial transverse plane (typically at z=0) to a parallel observation plane (at z>0) by decomposing the field into a continuous spectrum of plane waves via a two-dimensional Fourier transform.[1] This representation, also known as the angular spectrum representation, expresses the field as an integral over spatial frequencies k_x and k_y, where each component propagates with a z-directed wave vector k_z = \sqrt{k^2 - k_x^2 - k_y^2} (with k = 2\pi / \lambda the wavenumber in the medium), enabling exact computation of diffraction effects under the Helmholtz equation.[1] The method distinguishes between propagating waves (where k_z is real) and evanescent waves (where k_z is imaginary, leading to exponential decay), providing a unified framework for both near- and far-field propagation.[2] Introduced in the context of Fourier optics by Joseph W. Goodman in his 1968 textbook Introduction to Fourier Optics, the angular spectrum method builds on earlier diffraction theories like those of Kirchhoff and Rayleigh-Sommerfeld, offering a frequency-domain alternative that aligns with linear systems analysis.[3] The propagation kernel derives from the Helmholtz equation for the transverse field component, yielding the transfer function H(k_x, k_y, z) = \exp(i k_z z) for propagating components, which is applied after Fourier transformation of the input field U(x, y, 0).[2] In numerical implementations, fast Fourier transforms (FFTs) facilitate efficient computation on discrete grids, though care must be taken to handle evanescent waves and avoid aliasing via techniques like zero-padding or low-pass filtering.[3] The method's versatility extends beyond optics to acoustics and other wave phenomena, with applications in simulating diffraction patterns, holographic imaging, optical device design (such as photonic circuits and fiber delay lines), and modeling near-field effects like total internal reflection.[4] Recent advances, including band-extended and scaled variants, enhance accuracy for long-distance propagation and off-axis calculations, addressing limitations in traditional Fresnel or Fraunhofer approximations.[4] Its computational efficiency and physical insight into wave decomposition have made it a cornerstone for both theoretical analysis and practical simulations in modern optical engineering.[1]Fundamentals
Definition and overview
The angular spectrum method is a technique for modeling the propagation of scalar wave fields, such as light or acoustic waves, by decomposing the field into a continuous spectrum of plane waves with different propagation angles. This approach allows for the precise description of wave evolution in free space or homogeneous media without approximations beyond the scalar wave equation.[5] Physically, the method represents the wave field at an initial plane as a superposition of obliquely propagating plane waves, where each component is characterized by spatial frequencies that correspond to specific angles relative to the normal direction.[5] These plane waves propagate independently according to their wave vectors, and their recombination at a subsequent plane yields the propagated field. The underlying mathematical tool for this decomposition is the Fourier transform, which maps the spatial distribution to an angular spectrum.[5] The method originated in the mid-20th century amid developments in Fourier optics and diffraction theory, with foundational work by J. A. Ratcliffe in 1956[6] applying it to ionospheric propagation problems. Key contributions to its use in optical diffraction were made by Joseph W. Goodman in 1968, who integrated it into the framework of Fourier optics for analyzing imaging and coherence.[5] It offers general advantages as an exact solution for both paraxial and non-paraxial propagation in homogeneous media, while enabling computationally efficient numerical simulations through fast Fourier transform algorithms.[5]Relation to Fourier optics
The angular spectrum method is fundamentally grounded in the principles of Fourier optics, which conceptualizes diffraction as a linear filtering process within the spatial frequency domain. In this framework, the optical field distribution across a plane is represented by its two-dimensional Fourier transform, where each spatial frequency component corresponds to a plane wave propagating at an angle determined by the wave vector components k_x and k_y. This decomposition allows propagation to be treated as a multiplicative operation on these frequency components, aligning directly with Fourier optics' emphasis on frequency-domain analysis for understanding wave behavior.[7] A key connection arises through the convolution theorem, which underpins much of Fourier optics. Free-space propagation of the field can be expressed as a spatial convolution between the initial field and the impulse response function, or propagator, that describes how a point source spreads over distance z. In the angular spectrum domain, this convolution transforms into a straightforward multiplication by the phase factor \exp(i k_z z), where k_z = \sqrt{k^2 - k_x^2 - k_y^2} and k = 2\pi / \lambda is the wave number. This duality enables efficient computation and conceptual clarity, as the method leverages the theorem to model diffraction without direct spatial-domain integration.[7] Spatial frequency, defined as cycles per unit length (with dimensions of inverse length), serves as a prerequisite concept linking the method to propagation angles via the transverse wave vector components: k_x = k \sin\theta_x and k_y = k \sin\theta_y, where \theta_x and \theta_y are the angles relative to the optical axis. This mapping highlights how higher spatial frequencies correspond to steeper propagation angles.[8] Unlike conventional Fourier optics, which typically assumes the paraxial approximation for small angles (where \sin\theta \approx \theta and evanescent contributions are neglected), the angular spectrum method provides a more general formulation by including all spatial frequencies. For components where k_x^2 + k_y^2 > k^2, k_z becomes imaginary, resulting in evanescent waves that decay exponentially along the propagation direction rather than oscillating. This inclusion extends the applicability to non-paraxial scenarios, such as near-field effects or high-numerical-aperture systems, while retaining the core frequency-domain insights of Fourier optics.[8]Theoretical formulation
Angular spectrum representation
The angular spectrum representation provides a fundamental decomposition of the scalar wave field into a superposition of plane waves, each characterized by its direction of propagation. In this formulation, the complex scalar field u(x, y, 0) at a reference plane z = 0 is expressed as the two-dimensional inverse Fourier transform of the angular spectrum A(f_x, f_y), where f_x and f_y denote spatial frequencies in the x- and y-directions, respectively: u(x, y, 0) = \iint_{-\infty}^{\infty} A(f_x, f_y) \exp\left[i 2\pi (f_x x + f_y y)\right] \, df_x \, df_y. This representation originates from the plane wave expansion of the field and is a cornerstone of Fourier optics.[9] The angular spectrum A(f_x, f_y) itself is defined as the two-dimensional Fourier transform of the field at z = 0: A(f_x, f_y) = \iint_{-\infty}^{\infty} u(x, y, 0) \exp\left[-i 2\pi (f_x x + f_y y)\right] \, dx \, dy. The spatial frequencies f_x and f_y correspond to the directional components of the wave vectors, related to the propagation angles \theta_x and \theta_y by \sin \theta_x \approx \lambda f_x and \sin \theta_y \approx \lambda f_y under the paraxial approximation for small angles, where \lambda is the wavelength.[7][9] The components of the angular spectrum can be interpreted in terms of the longitudinal wave number k_z = \sqrt{(2\pi / \lambda)^2 - (2\pi f_x)^2 - (2\pi f_y)^2}. For propagating waves, where the transverse spatial frequency magnitude \sqrt{f_x^2 + f_y^2} < 1 / \lambda, k_z is real, corresponding to plane waves that carry energy away from the reference plane without decay. In contrast, evanescent waves arise when \sqrt{f_x^2 + f_y^2} > 1 / \lambda, making k_z imaginary; these components decay exponentially with distance from the plane and are associated with near-field effects, such as those in total internal reflection or subwavelength imaging scenarios.[7][9] This representation assumes the field is defined over an infinite transverse plane at z = 0, ensuring the integrals converge and the decomposition is complete. For practical scenarios with finite apertures, such as limited illumination or observation areas, edge effects introduce artifacts like Gibbs ringing in the spectrum, requiring careful handling in applications.[7] The angular spectrum at subsequent planes can then be obtained by multiplying A(f_x, f_y) by a propagation phase factor, enabling the modeling of free-space evolution.Free-space propagation
The angular spectrum method provides an exact analytical framework for describing the propagation of scalar electromagnetic fields in homogeneous, free-space media, governed by the Helmholtz equation. The field distribution u(x, y, [z](/page/Z)) at a longitudinal distance [z](/page/Z) from the initial plane (z = 0) evolves through the phase modulation of its angular spectrum components, representing a superposition of plane waves with varying transverse wavevectors. This approach treats propagation as a linear filtering operation in the spatial frequency domain, where each spectral component advances with its corresponding longitudinal phase factor derived from the dispersion relation of waves in free space.[8][7] The propagation transfer function H(f_x, f_y, z) multiplies the initial angular spectrum A(f_x, f_y, 0) to yield the spectrum at distance z: A(f_x, f_y, z) = A(f_x, f_y, 0) \exp\left[i k z \sqrt{1 - \lambda^2 (f_x^2 + f_y^2)}\right], where k = 2\pi / \lambda is the wavenumber, \lambda is the wavelength, and f_x, f_y are the spatial frequencies in the x and y directions, respectively. This expression arises from the plane-wave decomposition, where the longitudinal component of the wavevector k_z = k \sqrt{1 - \lambda^2 (f_x^2 + f_y^2)} ensures satisfaction of the wave equation (\nabla^2 + k^2) u = 0. For propagating waves, where \lambda^2 (f_x^2 + f_y^2) < 1, k_z is real, leading to oscillatory phase advancement.[8][7] The field at the propagated plane is reconstructed via the inverse two-dimensional Fourier transform of the evolved spectrum: u(x, y, z) = \iint_{-\infty}^{\infty} A(f_x, f_y, z) \exp\left[i 2\pi (f_x x + f_y y)\right] \, df_x \, df_y. This integral, often termed the angular spectrum propagator, sums the contributions of all plane waves, each tilted at angles \theta_x \approx \lambda f_x and \theta_y \approx \lambda f_y for small angles, to form the diffracted field. Equivalently, the propagation can be expressed directly in the spatial domain as a convolution with the propagator kernel, though the frequency-domain formulation is computationally preferable for its separability.[8][7] Evanescent waves arise when \lambda^2 (f_x^2 + f_y^2) > 1, making the argument of the square root negative; in this regime, k_z = -i \kappa with \kappa = k \sqrt{\lambda^2 (f_x^2 + f_y^2) - 1} > 0, transforming the transfer function to \exp[-\kappa z] for forward propagation (z > 0). These components, corresponding to high spatial frequencies or supercritical angles, decay exponentially away from the initial plane and do not contribute to far-field propagation, effectively filtering near-field details. This behavior is crucial for understanding resolution limits in imaging systems.[8][7] The method is rigorously valid for any propagation distance z in non-absorbing, homogeneous media, providing an exact solution to the scalar Helmholtz equation without approximations on the field or distance. In the paraxial limit of small angles (\lambda^2 (f_x^2 + f_y^2) \ll 1), the square root approximates to \sqrt{1 - \xi} \approx 1 - \xi/2 with \xi = \lambda^2 (f_x^2 + f_y^2), reducing the transfer function to the quadratic phase factor of the Fresnel diffraction integral and recovering the paraxial approximation. This exactness distinguishes the angular spectrum approach from approximate methods like Fresnel or Fraunhofer diffraction, enabling accurate modeling across near-, intermediate-, and far-field regimes.[8][7]Numerical implementation
Discrete Fourier transform approach
The discrete Fourier transform (DFT) approach provides a straightforward numerical method for implementing the angular spectrum propagation in free space, leveraging fast Fourier transform (FFT) algorithms for efficiency. This technique discretizes the continuous angular spectrum representation, where the input complex field u(x, y, 0) is sampled on a uniform 2D grid of size N \times N with spatial sampling intervals \Delta x and \Delta y (often equal for isotropic cases). The propagation to a distance z follows a three-step algorithm. First, the 2D DFT of the discretized input field yields the discrete angular spectrum A(m, n), where m, n = 0, 1, \dots, N-1.[10] In the second step, the discrete angular spectrum is multiplied element-wise by the transfer function H(m, n) = \exp\left[ i k z \sqrt{1 - \lambda^2 (f_x(m)^2 + f_y(n)^2)} \right], with wave number k = 2\pi / \lambda, spatial frequencies f_x(m) = m / (N \Delta x) and f_y(n) = n / (N \Delta y), and \lambda the wavelength; this discrete form approximates the continuous free-space propagator. The resulting spectrum is then transformed back via a 2D inverse DFT to obtain the output field u(x, y, z) on the same grid.[10][11] Appropriate sampling is critical to prevent aliasing and ensure accurate representation of diffracted fields. The spatial sampling must satisfy the Nyquist criterion \Delta x \leq \lambda / 2 to capture the highest spatial frequencies in the input field. Additionally, to avoid aliasing in the propagated field due to the nonlinear phase in the transfer function, the sampling interval should fulfill \Delta x \leq \lambda z / (N \Delta x), where N \Delta x is the total grid extent; the corresponding frequency sampling is \Delta f_x = 1 / (N \Delta x).[12][11] Zero-padding techniques enhance the implementation by extending the grid size beyond the input aperture, typically by appending zeros to reach a larger N_{\text{pad}} > N. This increases the effective field of view, improves output resolution by effectively reducing \Delta x in the reconstructed plane, and mitigates finite-grid artifacts such as wrap-around errors from periodic DFT assumptions. Padding factors of 2–4 are common, balancing accuracy and computational cost.[10][12] In practice, this DFT-based method is routinely implemented in software environments like MATLAB, utilizing built-in functions such asfft2 for the forward transform and ifft2 for the inverse, or in Python via NumPy's np.fft.fft2 and np.fft.ifft2. Both forward (z > 0) and backward (z < 0) propagation are handled by adjusting the sign in the transfer function exponent; for backward cases, conjugating the phase term equivalently reverses the direction while preserving numerical stability.[11][10]