Gaussian Beamlet Decomposition (GBD)

• Gaussian Beamlet Decomposition (GBD) is a computational method that represents arbitrary optical fields as coherent sums of paraxial Gaussian beamlets.
• It employs analytic propagation via ABCD matrices and tight phase-space tiling to capture complex diffraction and aberration effects in optical systems.
• Accuracy is controlled by beamlet density and grid parameters, with advanced implementations achieving significant speedups using GPU acceleration.

Gaussian Beamlet Decomposition (GBD) is a class of computational algorithms used to approximate the propagation of optical fields—particularly in high-frequency and highly oscillatory regimes—by representing arbitrary wavefronts as superpositions of spatially localized, paraxial Gaussian beams. These techniques combine analytic propagation of Gaussian beams with frame decompositions that tile phase space, yielding a versatile and efficient alternative to full diffraction integrals for modeling physical optics phenomena, especially in large and complex optical systems. GBD is rigorously founded in wave asymptotics and finds applications ranging from wave equation parametrices to physical optics modeling of astronomical telescopes, interferometers, and high-contrast imaging coronagraphy.

1. Foundations and Formulation

A core principle of GBD is the expansion of an arbitrary complex optical field as a coherent sum of spatially distributed fundamental Gaussian beams—termed "beamlets"—each characterized by a center (launch point), waist size, complex curvature, and amplitude-weighting. For a field on a plane (e.g., an entrance pupil):

$$E(x, y) \approx \sum_{i=1}^{G} b_i\, E_i(x - x_{i,0}, y - y_{i,0})$$

where each $E_i$ is a fundamental (possibly astigmatic) Gaussian function centered at $(x_{i,0}, y_{i,0})$, typically with common waist size $w_0$ but potentially varying in direction or curvature. The expansion coefficients $b_i$ are determined by projection or least-squares fit to the sampled input field on a chosen decomposition grid. The complex beam parameter for each direction takes the form:

$$q(z)^{-1} = \frac{1}{R(z)} + i\frac{\lambda}{\pi w^2(z)}$$

where $R(z)$ is the radius of curvature, $w(z)$ is the beam radius at the $1/e$ amplitude contour, and $\lambda$ the wavelength.

Higher-order (multi-scale) GBD frameworks extend the decomposition using tight frames in phase space, such as the wave-atom parametrix, by partitioning frequency and spatial support, enabling robust performance for oscillatory fields (Berra et al., 2017).

2. Propagation via ABCD Matrices and Ray Methods

Each Gaussian beamlet is propagated independently through the optical system using paraxial ray-transfer (ABCD) matrices. For an optical system described by the ABCD matrix:

$E_i$0

the complex beam parameter is updated as:

$E_i$1

For general astigmatic or misaligned systems, this update generalizes to $E_i$2 complex matrix equations for the local curvature tensor $E_i$3 (Ashcraft et al., 2021, Ashcraft et al., 2023). The ray-centric description allows GBD to capture system-level aberrations, misalignments, and vignetting directly, as the underlying rays follow the geometric propagation and each beamlet's field structure is updated accordingly.

A central computational advantage arises because, once beamlet parameters are propagated, their fields at the output plane can be evaluated analytically and summed coherently—bypassing the need for FFT-based convolution per surface and enabling efficient scaling for systems with many optical elements (Ashcraft et al., 2024).

3. Multi-Scale and Phase-Space Decomposition

Advanced GBD approaches construct tight frames based on localized Gaussian functions in space and frequency, e.g.,

$E_i$4

covering phase space by dyadic frequency partition $E_i$5 and spatial lattice $E_i$6 (Berra et al., 2017). This structure yields a decomposition supporting multi-resolution analysis, crucial for capturing both the oscillatory and smooth components of the field and enabling fine control over phase-space localization. Operators with slowly-varying symbols act nearly diagonally on this frame, and the expansion naturally adapts to the local frequency content. This underpins rigorous parametrices for the wave equation and efficient modeling of boundary value problems.

4. Accuracy, Error Metrics, and Benchmarking

GBD's accuracy is governed by grid parameters (number of beamlets, waist scaling, window/aperture coverage) and the nature of the target field (smooth, clipped, aberrated, etc.). Quantitative error metrics used in benchmarking include the discretized normalized mean-square error (DNMSE):

$E_i$7

as well as point-wise and summed relative errors (Zhao et al., 2022). For free-space propagation of unclipped Gaussian beams, DNMSE on the order of $E_i$8 to $E_i$9 is typical for grids of $(x_{i,0}, y_{i,0})$0 beamlets per dimension; for clipped beams and at hard edges, GBD often outperforms modal decompositions.

In the high-frequency regime, error bounds of $(x_{i,0}, y_{i,0})$1 for the parametrix, and $(x_{i,0}, y_{i,0})$2 for $(x_{i,0}, y_{i,0})$3-th order beams, are proven under non-trapping and smoothness assumptions, independent of the presence of caustics (Berra et al., 2017, Liu et al., 2013). For practical optical models, comparison to analytic Airy point spread functions demonstrates sum-normalized RMS errors of order $(x_{i,0}, y_{i,0})$4 for typical beamlet counts (Ashcraft et al., 2023).

5. Computational Considerations and Acceleration

Classical GBD implementations require $(x_{i,0}, y_{i,0})$5 complexity, as each beamlet's analytic contribution must be summed at each field point. This scales favorably compared to $(x_{i,0}, y_{i,0})$6 per-surface FFT-based Fresnel propagation when modeling multi-surface paths, since intermediate optical elements are traversed via inexpensive ABCD updates (Ashcraft et al., 2021).

Algorithmic innovations include the plane-evaluation method, wherein analytic expressions derived from the misaligned Collins integral allow all field points in the output plane to be computed with a single per-beamlet ray trace and a bulk analytic formula (Ashcraft et al., 2024). This yields up to $(x_{i,0}, y_{i,0})$7 acceleration in total runtime when leveraging GPU architectures, as demonstrated in contemporary open-source codes. Parallelization, advanced sampling schemes (e.g., Fibonacci spiral for circular apertures), and vectorized array operations further enhance throughput.

The overlap factor (OF), defined as $(x_{i,0}, y_{i,0})$8 for $(x_{i,0}, y_{i,0})$9 beamlets across width $w_0$0, is a critical parameter; optimal values depend on the desired spatial frequency fidelity (typically $w_0$1).

6. Applications, Extensions, and Limitations

GBD is widely utilized in physical optics simulation of complex astronomical instruments, including segmented-mirror telescopes and coronagraph systems. It enables rapid, high-fidelity modeling of point spread functions (PSFs), the study of wavefront errors, misalignment tolerances, and diffraction phenomena beyond the reach of traditional geometric or direct numerical wave propagation (Ashcraft et al., 2023, Ashcraft et al., 2021). In interferometry, GBD accurately treats clipped or aberrated beams and strongly curved macrosurfaces, surpassing mode-expansion methods in these regimes (Zhao et al., 2022).

Notable limitations are inherent to the paraxial approximation; extremely non-paraxial propagation or evanescent field contributions require alternative treatments. Capturing sharp-edge diffraction (e.g., spiders, segmentation gaps) with sufficient fidelity may necessitate high beamlet densities or enhancements such as truncated, flattened-Gaussian, or higher-order transverse mode bases. Polarization effects and full vectorial beamlet extensions are subjects of ongoing development (Ashcraft et al., 2023). Algorithmic flexibility and open-source implementations (e.g., "Poke" for PSF modeling) have enabled widespread adoption and continued methodological advancement (Ashcraft et al., 2024).

7. Relation to Alternative Approaches and Theoretical Context

GBD forms one of several families of high-frequency asymptotic and hybrid optics algorithms. Related approaches include:

• Modal expansions: Decomposition into Hermite–Gaussian or Laguerre–Gaussian modes (Xiao et al., 2018). These bases are complete for paraxial fields, but single-axis expansions can be less effective in complex, clipped, or aberrated scenarios, where the spatial grid-based GBD naturally resolves local field structure.
• Gaussian beam parametrices: In wave propagation, GBD connects to the construction of parametrices for wave and Helmholtz equations via boundary data or source-manifold integrations (Berra et al., 2017, Liu et al., 2013).
• Wave-atom/wave-packet methods: Multi-scale GBD elaborates on wave-atom based phase-space tilings, enabling controlled localization and efficient operator action in high-frequency problems (Berra et al., 2017).

The convergence properties, error control, and computational scalability make GBD a compelling methodology for simulating both linear wave equations in high-frequency limits and large-scale physical optics systems where traditional methods are computationally prohibitive. Researchers continue to refine beamlet representation, acceleration, and coupling to new physical effects for increasingly demanding applications.