Frozen Gaussian Sampling (FGS) Algorithm

• Frozen Gaussian Sampling (FGS) is a mesh-free Monte Carlo algorithm that approximates quantum dynamics by evolving localized Gaussian packets along classical trajectories.
• It employs importance sampling and ODE-driven packet propagation to reconstruct physical observables with sampling error independent of the semiclassical parameter.
• FGS effectively overcomes the limitations of mesh-based methods in high-dimensional simulations, and it has been applied to open quantum systems, wave equations, and nonadiabatic dynamics.

The Frozen Gaussian Sampling (FGS) algorithm is a mesh-free, Monte Carlo numerical method originally conceived for the simulation of high-frequency and semiclassical quantum dynamics. The approach leverages the frozen Gaussian approximation, representing quantum or wavefield solutions as superpositions of phase-space localizations that evolve along classical or stochastic trajectories. By performing Monte Carlo sampling over the space of Gaussian packet initial data, FGS attains observables with sampling error independent of the semiclassical parameter, thus fundamentally overcoming the "curse of dimensionality" that afflicts traditional grid- or mesh-based solvers in the small-parameter regime. It has been adapted for Markovian open quantum systems, scalar wave equations, semiclassical Schrödinger equations, nonadiabatic dynamics at metal surfaces, and related open system models (Xu et al., 16 Dec 2025, Chai et al., 2022, Huang et al., 2022, Xie et al., 2021, Lu et al., 2016, Wang et al., 2 Aug 2024).

1. Theoretical Foundations

FGS is based on the frozen Gaussian approximation (FGA), wherein the solution (wavefunction, Wigner distribution, or density operator) is expressed as an integral over phase-space of localized Gaussian packets. Each packet evolves according to an ensemble of ODEs or piecewise-deterministic processes derived from semiclassical Hamiltonian or dissipative dynamics.

In Markovian open quantum systems, as formalized for example via the Wigner–Fokker–Planck (WFP) equation, the density operator $\hat{\rho}(t)$ dynamics are recast in phase-space as a PDE coupled to Lindblad or Fokker–Planck terms. The frozen Gaussian ansatz represents the Wigner function $W(x, \xi, t)$ as a sum over evolved Gaussian packets:

$$W(t, z) \approx M^{-1}\sum_{j=1}^M W_{\text{FGS}}(t, z; q_j, p_j),$$

with each $W_{\text{FGS}}$ derived by evolving a Gaussian packet initiated at $(q_j, p_j)$ under the ODEs deduced from the WFP structure. This approach directly ties the FGS method to the theoretical framework of Gaussian superposition and classical or stochastic classical flows, as pioneered in the semiclassical limit (Xu et al., 16 Dec 2025, Xie et al., 2021).

2. Algorithmic Structure and Implementation

FGS unfolds in three principal stages:

1. Initialization: A set of $M$ i.i.d. tuples $(q_j, p_j)$ are sampled from an importance distribution $\pi$ on phase space, typically Gaussian and parameterized by the initial state's mean and covariance. The variance-reducing choice $\pi \propto |A(0, q, p)|$ is optimal for Gaussian or WKB initial data (Xie et al., 2021, Chai et al., 2022).
2. Packet Propagation: For each packet, solve the appropriate ODE system for packet center $(q_j(t), p_j(t))$, covariance $G_j(t)$, and amplitude $A_j(t)$. These ODEs encode drift, diffusion, dissipation, and potential quantum corrections. In nonadiabatic or open-system contexts, surface-hopping and Lindblad dissipative terms are included via additional SDEs, jump processes, or matrix-valued blocks (Xu et al., 16 Dec 2025, Huang et al., 2022, Lu et al., 2016).
3. Statistical Reconstruction: Physical observables are estimated by averaging the contributions of all packets. The expectation of any quadratic or polynomial observable $a(x, \xi)$ is given by

$$\langle a \rangle_{\mathrm{FGS}} \approx (1/M) \sum_{j=1}^M \int a(z) W_{\mathrm{FGS}}(t,z;q_j,p_j) \, dz,$$

where analytic integration is possible for various $a$ due to the explicit quadratic nature of the packet (Xu et al., 16 Dec 2025, Xie et al., 2021).

Variance reduction techniques include importance sampling (weights $w_j=\pi_0/\pi_j$ if $\pi \neq \pi_0$) and control variates (using exact solutions for harmonic portions of the potential) (Xu et al., 16 Dec 2025).

3. Error Analysis and Computational Complexity

A hallmark of FGS methodology is that the sampling error in observables is independent of the semiclassical parameter (e.g., $\varepsilon$ or $1/k$), in contrast to mesh/grid methods whose cost and error scale catastrophically with shrinking parameter (Xu et al., 16 Dec 2025, Chai et al., 2022, Xie et al., 2021).

For any Lipschitz or polynomial observable,

$\mathrm{RMS~error}~\epsilon_M = O(M^{-1/2}),$

where the prefactor is independent of $\varepsilon$ or $k$. For Gaussian initial data and observables, the variance bound is sharp and unaffected by the oscillatory regime (Chai et al., 2022, Xie et al., 2021).

Computational complexity per sample per step is $O(d^3)$ for a $d$-dimensional system (mainly due to $4d^2$ ODEs for the packet covariance evolution in the open-system case). Overall runtime to maintain fixed accuracy in observables is $O(M d^3)$, with $M = O(\delta^{-2})$ for statistical error $\delta$, again circumventing the $\varepsilon^{-d}$ scaling endemic to mesh-based methods (Xu et al., 16 Dec 2025, Chai et al., 2022).

A summary of FGS scaling properties:

Source	Sampling error dependence	Mesh error dependence
FGS	$O(M^{-1/2})$, $\varepsilon$-independent	none
Grid/methods	$O(N^{-1/2})$, $N=O(\varepsilon^{-d})$	$O(\varepsilon^{-d})$

4. Applications Across Quantum and Wave Regimes

FGS has been deployed across a range of quantum and classical wave propagation problems:

• Markovian open quantum systems: FGS in the Wigner-Fokker-Planck framework allows long-time simulation of dissipative dynamics, capturing relaxation to steady states even in non-harmonic potentials where analytic results are lacking (Xu et al., 16 Dec 2025).
• Scalar high-frequency wave equations: FGS provides an efficient solver where the grid size required by conventional approaches would be $O(k^{-d})$, with FGS breaking this scaling for Gaussian data (Chai et al., 2022).
• Semiclassical Schrödinger equations: FGS provides a mesh-free, $\varepsilon$-insensitive method for high-dimensional quantum wave function evolution and observable computation, validated in up to seven dimensions (Xie et al., 2021).
• Electron transfer and nonadiabatic surface hopping at metal surfaces: FGS incorporates both semiclassical packet propagation and surface-hopping jumps, achieving accuracy independent of $\varepsilon$ and the metal band number $N$ (Huang et al., 2022, Lu et al., 2016).
• Reduced dynamics in open quantum models (e.g., Caldeira–Leggett model): FGS, integrated with diagram-resumming methods such as the inchworm algorithm, enables direct computation of reduced density operators in extended environments (Wang et al., 2 Aug 2024).

In all these regimes, mesh-free nature fully eliminates boundary-induced instabilities typical of grid methods, maintaining stability over long-time propagation scenarios.

5. Representative Numerical Results

Key benchmarks of FGS's performance include:

• Harmonic oscillator and double-well potentials: The error in mean position and momentum ($\langle x \rangle, \langle \xi \rangle$) decays as $M^{-1/2}$ and is unaffected as $\varepsilon$ is reduced by two orders of magnitude ($1/16 \to 1/128$) (Xu et al., 16 Dec 2025).
• Steady-state convergence in strongly non-harmonic traps: FGS demonstrates relaxation to a unique stationary Wigner function even for sextic triple-well and other potentials outside the analytically tractable class (Xu et al., 16 Dec 2025).
• Breakdown comparison to time-splitting spectral-grid methods: FGS remains stable for $T=8$ with $M=3200$ where grid methods fail due to boundary artifacts (Xu et al., 16 Dec 2025).
• Multidimensional observables: Accurate high-dimensional means achieved with sampling complexity $O(4^m)$ for $m$ up to $7$, outperforming grid methods exponentially (Xie et al., 2021).
• Nonadiabatic transitions and electron transfer: FGS captures quantum-thermal and interference effects missed by classical surface hopping, with convergence unchanged in the simultaneous limit $\varepsilon\to0, N\to\infty$ (Huang et al., 2022).

6. Advantages, Limitations, and Extensions

Advantages:

• Sampling error and computational cost for observables are independent of the semiclassical parameter, terminating the prohibitive scaling in $\varepsilon\rightarrow 0$.
• Fully mesh-free and boundaryless, eliminating reflection artifacts and permitting unbounded dispersion.
• Linearly parallelizable in $M$ for both simulation and observable estimation.
• Rigorous ODE theory ensures stability and positive-definite covariances through all propagation.

Limitations:

• The core FGS decomposition is exact only for (near-)Gaussian initial states or those admitting a Gaussian mixture representation; general initial data require further analysis (Xu et al., 16 Dec 2025).
• For accurate statistics $M=10^3$–$10^4$ may be needed, and for high $d$, $M$ increases exponentially, albeit with a base significantly lower than $\varepsilon^{-1}$ (Xie et al., 2021).
• In nonadiabatic or open-system variants, surface-hopping extensions remain weak-coupling and closed-system; Lindblad-type generalizations are ongoing work (Huang et al., 2022).

Extensions such as the integration of FGS with diagrammatic series acceleration (e.g., via the inchworm method) offer promising directions for handling strongly coupled environments and open quantum systems (Wang et al., 2 Aug 2024).

7. Relation to Broader Methodologies

FGS is tightly connected to the Herman–Kluk propagator, the path integral formalisms for semiclassical wave propagation, and the broader class of semiclassical and quantum-classical hybrid algorithms. By providing a rigorous Monte Carlo bridge across these domains, it permits practical simulation of multi-dimensional, highly oscillatory quantum and wave fields with provably controlled error in regimes long regarded as computationally intractable. FGS thus constitutes a central tool in the numerical analysis of high-frequency asymptotics, nonadiabatic chemistry, and dissipative quantum theory.

Principal references: (Xu et al., 16 Dec 2025, Chai et al., 2022, Huang et al., 2022, Xie et al., 2021, Lu et al., 2016, Wang et al., 2 Aug 2024).