Frozen Gaussian Approximation (FGA)

• FGA is a semiclassical computational method that represents high-frequency wavefields as integrals over fixed-width Gaussian packets propagating along classical trajectories.
• The method bypasses severe mesh-size constraints and beam-spreading instabilities of traditional solvers, ensuring robust performance near caustics.
• Rigorous convergence results and diverse implementations, including Eulerian, Lagrangian, and Monte Carlo schemes, validate FGA across seismic, quantum, and elastic wave applications.

The Frozen Gaussian Approximation (FGA) is a semiclassical computational methodology for efficiently capturing high-frequency wave phenomena in linear partial differential equations with highly oscillatory solutions. FGA represents the wavefield as a sum or integral over phase-space of fixed-width (frozen) Gaussian wavepackets whose centers evolve along classical trajectories, with phase and amplitude propagated by ODEs derived via asymptotic expansions. FGA bypasses the severe mesh-size constraints of conventional grid-based solvers and avoids the beam-spreading instabilities and caustic singularities that limit geometric optics and Gaussian beam methods. Rigorous convergence theorems establish FGA as an $\mathcal{O}(\varepsilon)$-accurate approximation in the semiclassical parameter, and the method admits tractable Eulerian, Lagrangian, and stochastic (Monte Carlo) realizations across a range of hyperbolic and dispersive systems—including periodic media, elastic waves, seismic modeling, quantum molecular dynamics, and nonadiabatic dynamics.

1. Mathematical Formulation and Phase-Space Representation

FGA is built upon an integral ansatz expressing the solution to a high-frequency wave equation as a superposition over phase-space of fixed-width Gaussian wavepackets. For a semiclassical Schrödinger equation of the form

$$i\,\varepsilon\,\partial_t\psi^\varepsilon(t,x) = -\frac{\varepsilon^2}{2}\,\Delta_x\psi^{\varepsilon} + V(x)\,\psi^\varepsilon,$$

with $\varepsilon \ll 1$, the FGA takes

$$\psi^{\text{FGA}}(t,x) = (2\pi\varepsilon)^{-d} \int_{\mathbb{R}^d}\int_{\mathbb{R}^d} a(t,q,p) \exp\left[\frac{i}{\varepsilon} \Phi(t,x;q,p)\right] dq\,dp,$$

where the phase function is

$$\Phi(t,x;q,p) = S(t,q,p) + p \cdot (x - Q(t;q,p)) + \frac{i}{2}|x - Q(t;q,p)|^2,$$

and $(Q,P)$ solve classical Hamiltonian flow equations. In periodic media, the integrand is modulated by Bloch eigenfunctions, leading to an FGA-Bloch ansatz for multiscale problems (Delgadillo et al., 2015).

The key properties distinguishing FGA:

• Frozen width: Each Gaussian maintains a fixed (typically minimal-uncertainty) width, with no beam-spreading Riccati equation as in Gaussian beam methods.
• Phase-space propagation: Centers $(Q,P)$ evolve by Hamilton's equations for an appropriate band-dependent or system-dependent classical Hamiltonian, independent across phase-space points.
• Amplitude transport: Amplitudes $a(t,q,p)$ satisfy first-order transport ODEs (Liouville-type or their matrix generalizations), incorporating Jacobian contributions from the linearization of the Hamiltonian flow.

In scalar high-frequency wave and hyperbolic systems, the ansatz generalizes to vector-valued or band-decomposed forms, often with explicit projection onto spectral subspaces or Bloch bands (Lu et al., 2010, Lu et al., 2010, Delgadillo et al., 2015).

2. Asymptotic Derivation and Evolution Equations

The FGA ansatz is derived via asymptotic analysis, inserting into the governing PDE and matching coefficients at consecutive powers of the semiclassical parameter $\varepsilon$. Leading-order consistency yields Hamilton's equations for the centers, while next-to-leading order imposes Liouville-type transport on amplitudes. For a single Bloch band in periodic media (Delgadillo et al., 2015):

• Hamiltonian flow:

$$\dot{Q}_n = \nabla_p E_n(P_n), \quad \dot{P}_n = -\nabla_q U(Q_n).$$

• Phase (action):

$$\dot{S}_n = P_n\cdot\dot{Q}_n - (E_n(P_n) + U(Q_n)), \qquad S_n(0) = 0.$$

• Amplitude transport:

$$\frac{d}{dt}a_n + \nabla_p H_n \cdot \nabla_q a_n - \nabla_q H_n \cdot \nabla_p a_n + \frac{1}{2} a_n\,\mathrm{Tr}\left[\partial^2_{q p} H_n\right] = 0.$$

The "frozen" aspect refers to the identity-matrix imaginary part in the Gaussian quadratic exponent; no additional evolution of a beam-modulus matrix is required (in contrast to thawed Gaussian or Gaussian beam propagation) (Yang et al., 2013, Lu et al., 2010).

These evolution equations extend to matrix-valued amplitudes for hyperbolic systems (using projections onto system eigenvectors (Lu et al., 2010)) and include nonadiabatic transitions for quantum matrix Schrödinger systems through surface-hopping extensions (Lu et al., 2016, Lu et al., 2016).

3. Numerical Algorithms and Parallel Implementation

The FGA algorithm proceeds through:

• Initial decomposition: Projecting the initial condition onto the frozen Gaussian wavepackets, typically via the Fourier-Bros-Iagolnitzer (FBI) or windowed-Gabor transform in the continuous setting, or via quadrature/sampling in high dimensions (Yang et al., 2013, Delgadillo et al., 2015).
• Phase-space propagation: Each sampled phase-space point evolves independently by ODE integration along the prescribed Hamiltonian flow, with amplitude ODEs solved concurrently.
• Final reconstruction: The wavefield is reconstituted as a sum/integral of Gaussian contributions at the desired spatial points and time.

This decomposition enables embarrassingly parallel ODE computation—the evolution of each beam is decoupled—facilitating scalable implementations on distributed-memory or GPU architectures (Yang et al., 2013, Hateley et al., 2018).

For mesh-free, high-dimensional problems and observables, Monte Carlo ("frozen Gaussian sampling" or FGS) versions of FGA select beams/samples via importance sampling from the modulus of the initial Gaussian transform (Xie et al., 2021, Chai et al., 2022, Huang et al., 2022). Such FGS estimators offer independence of required sample size from $\varepsilon$ and improved scaling with spatial dimension, validated by rigorous sampling error bounds and practical tests (Xie et al., 2021, Chai et al., 2022).

4. Convergence, Error Control, and Theoretical Guarantees

Rigorous convergence analysis for FGA in scalar and strictly hyperbolic systems establishes $L^2$ accuracy $$\|u(t,\cdot) - u^{\text{FGA}}(t,\cdot)\| = \mathcal{O}(\varepsilon)$$, assuming asymptotically high-frequency initial data and smooth system coefficients or isolated spectral bands (Lu et al., 2010, Lu et al., 2010, Delgadillo et al., 2015):

• For strictly hyperbolic systems and semiclassical Schrödinger or elastic wave equations, primary error sources are in the initial decomposition and neglect of higher-order terms in the amplitude expansion.
• In periodic media, projection onto an isolated Bloch band is required, and the convergence theorem holds uniformly if the inter-band gap remains non-vanishing (Delgadillo et al., 2015).
• For matrix Schrödinger systems with surface hopping (nonadiabatic quantum dynamics), the FGA with stochastic path-integral representation achieves $\mathcal{O}(\varepsilon)$ accuracy in $L^2$ for bounded time intervals (Lu et al., 2016, Lu et al., 2016).

Higher-order corrections to the FGA expand the amplitude in $\varepsilon$, yielding $\mathcal{O}(\varepsilon^K)$ approximations (Lu et al., 2010).

In Monte Carlo (FGS) realizations, the statistical (sampling) error is provably independent of $\varepsilon$ for Gaussian initial data and grows only gently (sublinearly) in spatial dimension for observables, resulting in mesh-free algorithms suitable for high-$d$ applications (Xie et al., 2021, Chai et al., 2022).

5. Comparison to Related Methods and Extensions

FGA generalizes the Herman-Kluk propagator (quantum mechanics), retains key advantages over geometric optics and Gaussian beam methods, and connects directly to Wigner/Husimi phase-space representations. Key comparison points:

• Geometric optics (WKB): Fails at caustics, where amplitudes diverge. FGA remains valid and accurate through caustics by superposing complex Gaussians, smoothing out singularities (Lu et al., 2010).
• Gaussian beam: Requires evolving both centers and width matrices; susceptible to beam-spreading, breakdown at caustics, and Riccati numerical instabilities. FGA's frozen width circumvents instability and maintains localization (Yang et al., 2013, Hateley et al., 2018).
• Eulerian FGA: Reformulates ODE transport in (Q,P) as evolution equations for phase, action, and amplitude on a fixed phase-space grid, sidestepping issues of Lagrangian ray divergence or under-resolution (Lu et al., 2010).

Extensions of FGA include:

• Imaginary time propagation for quantum statistical problems, e.g., rare-gas clusters, giving access to thermodynamic and structural quantities at low temperatures (Cartarius, 2015).
• Fractional Schrödinger equations, which introduce non-quadratic kinetic energies; extension requires regularization to avoid singularities in amplitude transport (Chai et al., 2024).
• Nonadiabatic dynamics: Surface-hopping FGA (FGA-SH) for matrix PDEs governing quantum-classical systems, with stochastic trajectory ensembles and explicit branching algorithms (Lu et al., 2016, Lu et al., 2016, Huang et al., 2022).
• High-efficiency reconstruction: On-grid correction algorithms (e.g., frozen Gaussian grid-point correction, FGGC) address the bottleneck of reconstructing solutions from non-grid-aligned wavepackets, preserving $\mathcal{O}(\varepsilon)$ asymptotics (Chai et al., 30 Apr 2025).

6. Applications and Numerical Demonstrations

FGA has been systematically applied and validated in:

• High-frequency seismic modeling, forward and inverse, for both acoustic and elastic wave equations. FGA reproduces main arrival times and amplitudes in complex models such as Marmousi, with efficiency and tunability for fast large-scale wave simulations (Yang et al., 2013, Hateley et al., 2018, Hateley et al., 2018). Elastic FGA handles P/S wave amplitudes and interface conversion (Hateley et al., 2019, Hateley et al., 2018).
• Seismic tomography and full-waveform inversion (FWI): Facilitates efficient forward and adjoint operator computation, P/S conversion, and is highly parallelizable (Hateley et al., 2018).
• Deep learning data generation for seismic imaging: FGA enables rapid production of high-frequency synthetic displacement fields for neural network training, with correct traveltimes and qualitative amplitude behavior (Hateley et al., 2018).
• Periodic media (photonic crystals, solid-state, band-structure propagation): FGA captures multiscale oscillations efficiently via Bloch-FGA, with rigorous convergence for isolated bands (Delgadillo et al., 2015, Delgadillo et al., 2015).
• Quantum molecular dynamics: FGA- and FGA-SH-type approaches address nonadiabatic electron-nuclear coupling, validated using Tully models and extended to metals and open quantum systems via inchworm Dyson resummation (Huang et al., 2022, Lu et al., 2016, Wang et al., 2024).
• Classical and quantum statistical mechanics: Imaginary time FGA computes partition functions and structural quantities with small error relative to fully coupled path integrals (Cartarius, 2015).

7. Limitations, Open Problems, and Algorithmic Advances

While FGA is highly efficient and robust, several limitations and ongoing advances are notable:

• Amplitude accuracy: The frozen-width ansatz introduces $\mathcal{O}(\varepsilon)$ amplitude errors that can persist even under mesh refinement; these are typically negligible for phase-dominated observables but may require higher-order corrections for sensitive amplitude features.
• Handling near-zero momentum: Asymptotic justification requires high-frequency (large $|p|$) initial data; cutoffs are imposed to exclude low-frequency bands and near-equilibrium points (Lu et al., 2010, Chai et al., 2024).
• Beam-number scaling and reconstruction cost: In high dimensions or at small $\varepsilon$, the number of beams and associated summation complexity can be significant; FGGC and FGS address these via least-squares grid correction and mesh-free Monte Carlo (Chai et al., 30 Apr 2025, Xie et al., 2021).
• Nonlinear and strongly coupled systems: FGA remains primarily a linear-PDE and weakly nonlinear method; application to strongly nonlinear or fully coupled systems often requires further approximation or hybridization.
• Adaptive and surrogate approaches: Active research includes learning-based acceleration, adaptive phase-space refinement, and surrogate models for ODE+amplitude mapping to enable FGA in very high-dimensional settings (Chai et al., 30 Apr 2025).

Ongoing work focuses on higher-order FGA expansions, generalized FGA for complex dispersive and nonlinear systems, rigorous analysis of stochastic surface-hopping variants, and broadening the algorithmic toolbox for efficient semiclassical computation in scientific and engineering practice.