Paper deep dive
From Basis to Basis: Gaussian Particle Representation for Interpretable PDE Operators
Zhihao Li, Yu Feng, Zhilu Lai, Wei Wang
Intelligence
Status: succeeded | Model: Gemma-4-26B-A4B | Prompt: intel-v1 | Confidence: 91%
Last extracted: 7/20/2026, 12:37:33 PM
Summary
The paper introduces the Gaussian Particle Operator (GPO), a neural operator framework for solving Partial Differential Equations (PDEs) that uses a learnable Gaussian particle basis to represent fluid fields. This approach offers interpretable, mesh-agnostic state representation with near-linear computational complexity, addressing limitations of traditional neural operators and Transformers regarding interpretability and high-frequency structure modeling.
Entities (8)
Relation Signals (7)
Gaussian Particle Operator → uses → Gaussian Basis
confidence 95% · We propose representing fields with a Gaussian basis... Building on this representation, we introduce a Gaussian Particle Operator
Gaussian Particle Operator → achieves → Near-linear complexity
confidence 92% · achieves near-linear complexity in N for a fixed modal budget
Gaussian Particle Operator → implements → Petrov-Galerkin Attention
confidence 92% · PG Gaussian Attention enables global cross-scale coupling... acts in modal space
Gaussian Particle Operator → evaluateson → NS3D
confidence 90% · On standard PDE benchmarks... NS3D... our method attains state-of-the-art competitive accuracy
Gaussian Particle Operator → evaluateson → NS2D
confidence 90% · On standard PDE benchmarks... NS2D... our method attains state-of-the-art competitive accuracy
Gaussian Basis → provides → Interpretability
confidence 90% · providing intrinsic interpretability... directly visualizable state
Gaussian Particle Operator → evaluateson → ERA5
confidence 88% · On standard PDE benchmarks and real datasets (including ERA5...)
Cypher Suggestions (0)
No Cypher suggestions yet.
Abstract
Abstract:Learning PDE dynamics for fluids increasingly relies on neural operators and Transformer-based models, yet these approaches often lack interpretability and struggle with localized, high-frequency structures while incurring quadratic cost in spatial samples. We propose representing fields with a Gaussian basis, where learned atoms carry explicit geometry (centers, anisotropic scales, weights) and form a compact, mesh-agnostic, directly visualizable state. Building on this representation, we introduce a Gaussian Particle Operator that acts in modal space: learned Gaussian modal windows perform a Petrov-Galerkin measurement, and PG Gaussian Attention enables global cross-scale coupling. This basis-to-basis design is resolution-agnostic and achieves near-linear complexity in N for a fixed modal budget, supporting irregular geometries and seamless 2D-to-3D extension. On standard PDE benchmarks and real datasets, our method attains state-of-the-art competitive accuracy while providing intrinsic interpretability.
Tags
Links
- Source: https://arxiv.org/abs/2602.21551v1
- Canonical: https://arxiv.org/abs/2602.21551v1
Trouble viewing inline? Open PDF directly →
Full Text
64,058 characters extracted from source content.
Expand or collapse full text
From Basis to Basis: Gaussian Particle Representation for Interpretable PDE Operators Zhihao Li Yu Feng Zhilu Lai Wei Wang Abstract Learning PDE dynamics for fluids increasingly relies on neural operators and Transformer-based models, yet these approaches often lack interpretability and struggle with localized, high-frequency structures while incurring quadratic cost in spatial samples. We propose to represent fields with a Gaussian basis, where learned atoms carry explicit geometry (centers, anisotropic scales, weights) and form a compact, mesh-agnostic, directly visualizable state. Building on this representation, we introduce a Gaussian Particle Operator that acts in modal space: learned Gaussian modal windows perform a Petrov–Galerkin measurement, a PG Gaussian Attention effects global cross-scale coupling. This basis-to-basis design is resolution-agnostic and achieves near-linear complexity in N for fixed modal budget, supporting irregular geometries and seamless 2D→ 3D extension. On standard PDE benchmarks and real datasets, our method attains state-of-the-art–competitive accuracy while providing intrinsic interpretability. Machine Learning, ICML 1 Introduction Fluid-governed PDEs (Wazwaz, 2002; Gurtin, 1982) underpin critical real-world systems, from numerical weather prediction and climate reanalysis to ocean circulation and engineering aerodynamics (McKeown et al., 2020; Shlesinger et al., 1987). Classical solvers (finite element/volume and spectral methods) (Wazwaz, 2002; Johnson, 2012; Klaasen and Troy, 1984) deliver high fidelity but face persistent challenges: strongly multi-scale dynamics, mesh dependence and complex geometries, stiffness in time integration, and high computational cost for long rollouts. Neural operators (Li et al., 2021; Kovachki et al., 2023; Lu et al., 2021) emerged as data-driven maps between function spaces, enabling resolution-agnostic surrogates; more recently, Transformer-based operators (Cao, 2021; Hao et al., 2023) leverage attention to capture long-range interactions and achieve strong empirical performance on diverse PDE tasks. However, these models still suffer from two key limitations: (i) poor interpretability—latent features and attention weights are rarely tied to physically meaningful modes; and (i) localization/frequency bias—global self-attention tends to favor low-rank, low-frequency correlations, making sharp fronts, vortical filaments, and other high-frequency structures harder to model, while naïvely scaling attention over N spatial samples incurs (N2)O(N^2) cost (Li et al., 2025b). We advocate representing fluid fields with a Gaussian (particle) basis rather than a fixed grid, hand-picked spectra (Gupta et al., 2021; Li et al., 2025a), or a monolithic implicit network (Serrano et al., 2023, 2024). Gaussian atoms carry explicit geometry—centers and (anisotropic) scales—which align naturally with coherent flow structures (vortices, filaments, fronts), afford multiscale locality, and are directly visualizable and differentiable. This basis is meshagnostic and compact, supports irregular boundaries, and extends seamlessly from 2D to 3D (Buhmann, 2000; Park and Sandberg, 1991). While prior neural representations often rely on global Fourier features, wavelets, or blackbox INRs, learning a particleized Gaussian basis as the primary state of the field has been scarcely explored and offers a clearer bridge to physical intuition. Concretely, a field is approximated by weighted Gaussians with μi _i (centers), σi _i (scales, possibly anisotropic), and wiw_i (mixture weights) learned from data; evaluating these atoms at query locations yields a compact coefficient vector that serves as the field’s interpretable latent state (basis). Figure 1: Overview of the Gaussian Basis Field framework. Given the physics field a(x)a(x), an encoder EθE_θ produces G Gaussian components per spatial location (mean μ, scale σ, and mixture weight w). These define a Gaussian field that is evaluated at queries to form the basis z, which is decoded by a fixed decoder fϕdecf_φ^dec to reconstruct the output field. We present an interpretable, resolution-agnostic neural operator that learns a Gaussian particle basis for fields and couples it with a Petrov–Galerkin Gaussian Attention layer, enabling basis-to-basis modeling with near-linear complexity and strong accuracy on 2D/3D and irregular domains. Our contributions are summarized as follows: (1) Gaussian Particle Representation. An encoder learns per-site Gaussians (μ,σ,w)(μ,σ,w); evaluating at arbitrary queries yields an interpretable, visualizable basis Z that is mesh-agnostic and extends seamlessly to 3D. (2) PG Gaussian Attention. Learned Gaussian modal windows perform PG-style measurement (N→GN\!→\!G), a G×G×G attention implements the modal kernel (global coupling), and the result is scattered back (G→NG\!→\!N), yielding a principled and interpretable operator. (3) Efficiency & Scalability. With a small modal budget G≪NG\! \!N, spatial transfers scale (N)O(N) and modal attention is independent of N, delivering near-linear growth with resolution and supporting multi-step operator stacking. (4) Empirical validation. Across standard PDE benchmarks and real datasets (including ERA5 and 3D/irregular domains), our approach attains state-of-the-art–competitive accuracy while providing intrinsic interpretability (particle and modal diagnostics), yielding improved spectral fidelity and rollout stability—demonstrating a favorable accuracy–interpretability trade-off. 2 Methodology Figure 2: Architecture of the Gaussian Particle Operator (GPO). The pipeline encodes a()a(x) into a Gaussian field (μ,σ,w)(μ,σ,w), evaluates a basis Z, applies the modal operator ψO_ψ, and decodes to u^() u(x); black arrows denote forward computation and red arrows denote gradients. 2.1 Gaussian Basis Representation From physics field to Gaussian field and basis. We represent a spatial field a:Ω⊂ℝd→ℝdaa: \!⊂\!R^d\!→\!R^d_a by a set of Gaussian particles placed at each sample location jx_j. Each particle is parameterized by a center μj,i∈ℝd _j,i\!∈\!R^d, axis-aligned scale σj,i∈ℝd _j,i\!∈\!R^d, and mixture weight wj,i∈[0,1]w_j,i\!∈\![0,1] with ∑i=1Gwj,i=1 _i=1^Gw_j,i=1. The associated (unnormalized) kernel is G(~;μj,i,σj,i)=exp(−12‖(~−μj,i)/σj,i‖22).G( x; _j,i, _j,i)= \! (- 12 \|( x- _j,i)/ _j,i \|_2^2 ). (1) Evaluating these kernels at query ~j x_j yields the Gaussian basis coefficients zj,i=wj,iG(~j;μj,i,σj,i),z_j,i=w_j,i\,G( x_j; _j,i, _j,i), (2) where j=[zj,1,…,zj,G]⊤∈ℝGz_j=[z_j,1,…,z_j,G] ^G. Physically, the Gaussian field acts as a mollified, locally supported expansion of a(⋅)a(·); the coefficients in Eq.(2) can be viewed as localized averages of a under data-adaptive windows (μ,σ)(μ,σ), while w distributes mass among overlapping particles. This basis is resolution-agnostic and naturally extends to irregular geometries. Encoder. Given samples (xj,aj)j=1N\(x_j,a_j)\_j=1^N, the encoder EθE_θ injects geometry by embedding coordinates and fusing them with field features. Concretely, we form ηj=[aj,γ(xj)] _j=[\,a_j,\;γ(x_j)\,], where γ(⋅)γ(·) is a fixed Fourier-feature map (or a small MLP), and compute ϕj _j =ReLU(Winηj), =ReLU(W_in _j), (3) μj _j =WμReLU(Uμϕj), =W_μ\,ReLU(U_μ _j), (4) σj _j =Softplus(WσReLU(Uσϕj)), =Softplus\! (W_σ\,ReLU(U_σ _j) ), (5) wj w_j =Softmax(WwReLU(Uwϕj)), =Softmax\! (W_w\,ReLU(U_w _j) ), (6) reshaped as μj,i,σj,i∈ℝd _j,i, _j,i ^d and wj,i∈[0,1]w_j,i∈[0,1] with ∑iwj,i=1 _iw_j,i=1. We use γ(x)γ(x) as a fixed positional encoding (Fourier features): γ(x)=[sin(2πBx),cos(2πBx)]γ(x)=[ (2π Bx), (2π Bx)], where B∈ℝm×dB ^m× d is sampled once from (0,σB2)N(0, _B^2) and kept fixed; the resulting coordinate embedding has dimension 2m2m. Gaussian basis evaluation. With (μ,σ,w)(μ,σ,w) predicted by EθE_θ, the weighted Gaussian evaluation Eq.(2) produces the per-site latent vector jz_j. Physically, μ encodes particle locations, σ controls receptive-field sizes (anisotropy along axes), and w balances overlapping contributions. Computationally, the map (μ,σ,w,~)↦(μ,σ,w, x) is local and embarrassingly parallel. Decoder. A lightweight MLP head fϕdec:ℝG→ℝcoutf^dec_φ:R^G\!→\!R^c_out regresses from jz_j to the field value at the query: a^(~j)=fϕdec(j). a( x_j)=f^dec_φ(z_j). (7) In practice, we use a two-layer perceptron with ReLU. Gaussian particle regularization. Because the constraints act on the particle parameters produced by the encoder, we impose them at the Gaussian-field level (conceptually tied to EθE_θ but applied after parameter prediction): ℒμ _μ =1N∑j=1N‖∑i=1Gwj,iμj,i−j‖22, = 1N _j=1^N \| _i=1^Gw_j,i _j,i-x_j \|_2^2, (8) ℒσ _σ =1NGd∑j,i,ℓ[[σj,i,ℓ−σmax]++[σmin−σj,i,ℓ]+], = 1NGd _j,i, [\,[ _j,i, - _ ]_++[ _ - _j,i, ]_+ ], (9) which promote spatial interpretability (centers near coordinates), avoid degenerate particles, and discourage overly peaky mixtures. Overview pipeline. Eqs.( 2–7) define the Gaussian-field training pipeline, and the complete forward/backward diagram is in Figure 1. We minimize the reconstruction loss together with ℒμL_μ, ℒσL_σ. Approximation capacity of the Gaussian basis. We record a standard density result: Lemma 2.1 (Density of Gaussian mixtures). Let Ω⊂ℝd \!⊂\!R^d be compact. Finite linear combinations of (possibly anisotropic) Gaussian kernels are dense in C(Ω)C( ) and in Lr(Ω)L^r( ) for 1≤r<∞1≤ r<∞. Consequently, for any continuous scalar field v and ε>0 >0, there exist G and parameters (μi,σi,ci)i=1G\( _i, _i,c_i)\_i=1^G with ci∈ℝc_i such that ∥v(⋅)−∑i=1Gciexp(−12∥(⋅−μi)/σi∥22)∥∞<ε. \|v(·)- _i=1^Gc_i\, \! (- 12\|(·- _i)/ _i\|_2^2 ) \|_∞< . Vector-valued fields admit componentwise approximation. (Proof in Appx.B.1) 2.2 Petrov–Galerkin Gaussian Attention Sec 2.1 defines local Gaussian bases: at each location j, we obtain G weighted coefficients j∈ℝGz_j ^G from particles (μ,σ,w)(μ,σ,w). To learn a resolution-agnostic operator, we adopt a Petrov–Galerkin (PG) view—expanding in a trial space and testing with a (possibly different) test space—and move from the spatial grid to modal space (Franca et al., 2006; Brooks and Hughes, 1982). Learned Gaussian modal windows first measure the field by pooling across locations, a global mode coupling acts on the G components, and the result is scattered back to locations. This PG pipeline yields an efficient and interpretable basis-to-basis operator. 2.2.1 From Petrov–Galerkin projection to a Gaussian-basis operator In Petrov–Galerkin, we expand the field in a trial space and evaluate it with a test space. Here, Gaussian particles form the trials, while Gaussian modal windows serve as discrete tests that pool local information into global modes. Trial functions (Gaussian basis). Let the unnormalized Gaussian particle (anchored at location j, component i) be ϕj,i()=exp(−12‖(−μj,i)/σj,i‖22),Σj,i=diag(σj,i2). _j,i(x)= \! (- 12 \|(x- _j,i)/ _j,i \|_2^2 ), _j,i=diag( _j,i^2). (10) With weighted evaluations (Sec. 2.1), each site provides j=[zj,1,…,zj,G]⊤∈ℝGz_j=[z_j,1,…,z_j,G] ^G, where zj,i=wj,iϕj,i(~j)z_j,i=w_j,i\, _j,i( x_j). Test functions (Gaussian modal windows). We distinguish two normalizations: (i) row-normalized assignments ∑gpj,g=1 _gp_j,g=1 (used to map each location to modes), and (i) column-normalized weights p¯j,g=pj,g/∑j′pj′,g p_j,g=p_j,g/ _j p_j ,g so that ∑jp¯j,g=1 _j p_j,g=1 (used for PG measurement). ψg(x)≈∑j=1Np¯j,gδ(x−x~j),p¯j,g=pj,g∑j′pj′,g. _g(x)≈ _j=1^N p_j,g\,δ(x- x_j), p_j,g= p_j,g _j p_j ,g. (11) Using a linear projection of coefficients j=jWz∈ℝDs_j=z_jW_z ^D, the PG measurement (test of the trial field) yields modal tokens tg=∑j=1Np¯j,gsj,sj=zjWz.t_g= _j=1^N p_j,g\,s_j, s_j=z_jW_z. (12) Modal coupling and scatter. Let κ:1,…,G2→ℝD×Dκ:\1,…,G\^2 ^D× D be a (learned) coupling kernel over modes. PG updates the modal state and scatters it back: Ug=∑g′=1Gκ(g,g′)tg′∈ℝD, U_g= _g =1^Gκ(g,g )\,t_g ^D, (13) ~j=(∑g=1Gpj,gUg)Wout∈ℝG. z_j= ( _g=1^Gp_j,g\,U_g )W_out ^G. (14) Stacking sites gives Z~∈ℝN×G Z ^N× G. Algebraically, Z~≈A⊤(ZWz)Wout, \;\; Z\;≈\;A\,K\,A \,(Z\,W_z)\,W_out, (15) where A[j,g]=pj,gA[j,g]=p_j,g and [g,g′]K[g,g ] encodes modal coupling. Thus PG supplies the structure: test (measure) → couple → scatter. Table 1: Performance comparison with baselines on benchmarks. L2L_2 loss is recorded. Model NS2D NS3D ERA5-temp ERA5-wind u Carra Airfoil Turbulent PlanetSWE (Geo-)FNO 3.24E-02 5.07E-01 7.09E-03 1.02E-01 3.50E-01 5.12E-03 4.82E-01 9.23E-02 M2NO 2.76E-02 4.65E-01 3.35E-03 7.15E-02 3.62E-01 1.98E-03 4.53E-01 8.96E-02 CNO 9.27E-02 8.12E-01 1.24E-02 1.86E-01 6.89E-01 4.65E-03 7.91E-01 1.57E-01 LSM 3.11E-02 3.80E-01 5.86E-03 8.23E-02 4.05E-01 5.37E-03 4.16E-01 9.58E-02 AMG 3.12E-02 4.32E-01 6.11E-03 8.79E-02 3.21E-01 3.09E-03 4.47E-01 8.14E-02 GT 8.81E-02 5.39E-01 5.44E-03 1.55E-01 3.73E-01 8.26E-03 5.41E-01 1.23E-01 GNOT 7.19E-01 1.01E+00 1.55E-02 3.49E-01 7.57E-01 1.53E-02 9.86E-01 2.14E-01 Transolver 3.76E-02 5.29E-01 4.18E-03 1.06E-01 3.76E-01 9.46E-03 5.33E-01 1.05E-01 ONO 4.26E-02 8.83E-01 1.45E-02 3.49E-01 7.25E-01 5.26E-03 8.27E-01 1.69E-01 LNO 4.81E-02 3.68E-01 7.32E-03 1.31E-01 4.36E-01 4.39E-03 5.69E-01 9.12E-02 GPO 3.02E-02 3.44E-01 2.26E-03 6.68E-02 2.97E-01 1.12E-03 3.95E-01 8.01E-02 2.2.2 Attention as a parameterization of the PG operator We now instantiate Eq.( 15) with a multi-head attention layer that is global in modal space and local in the N↔GN\! \!G transfers. Let Z∈ℝN×GZ ^N× G and particle parameters (μ,σ,w)∈ℝN×G×d×ℝN×G×d×ℝN×G(μ,σ,w) ^N× G× d×R^N× G× d×R^N× G. Learned Gaussian modal windows. For head h, form a per-site descriptor ξj=[j,j,μj,σj]∈ℝG(2d+2) _j=[\,z_j,w_j, _j, _j\,] ^G(2d+2) and project to hj(h)∈ℝDh_j^(h)\!∈\!R^D. A softmax over modes produces windows pj,g(h)=softmaxg(Wp(h)hj(h)),p^(h)_j,g=softmax_g\! (W_p^(h)h^(h)_j ), (16) which instantiate the PG test functions in discrete form. Wp(h)∈ℝG×DW_p^(h)\!∈\!R^G× D is the (head-specific) linear projection that maps the local embedding hj(h)h^(h)_j at location j to mode logits over the G Gaussian modes. PG measurement N→GN\!→\!G. Project coefficients j(h)=jWz(h)s^(h)_j=z_jW^(h)_z and compute tokens tg(h)=∑jpj,g(h)j(h)∑jpj,g(h)∈ℝD,T(h)=[t1h,…,tGh]∈ℝG×D,t^(h)_g= _jp^(h)_j,g\,s^(h)_j _jp^(h)_j,g ^D,T^(h)=[t_1^h,…,t_G^h] ^G× D, (17) which matches the PG measurement in Eq.( 12). Global modal coupling (G×G×G attention). Scaled dot-product attention parameterizes the kernel K: Q(h) Q^(h) =T(h)WQ(h),K(h)=T(h)WK(h),V(h)=T(h)WV(h), =T^(h)W^(h)_Q,K^(h)=T^(h)W^(h)_K,V^(h)=T^(h)W^(h)_V, (18) α(h) α^(h) =softmax(Q(h)K(h)⊤D), =softmax\! ( Q^(h)K^(h)^\! D ), (19) T~(h) T^(h) =α(h)V(h)∈ℝG×D. =α^(h)V^(h)\;∈\;R^G× D. (20) Here α(h)(g,g′)α^(h)(g,g ) plays the role of a data-driven modal coupling κ(g,g′)κ(g,g ). Scatter G→NG\!→\!N and readout. Using the same windows, scatter the coupled modes back and read out to G coefficients: yj(h) y^(h)_j =∑gpj,g(h)Ug(h)∈ℝD, = _gp^(h)_j,g\,U^(h)_g ^D, (21) ~j z_j =(∥h=1Hyj(h))Wout∈ℝG, = ( \|_h=1^Hy^(h)_j )W_out ^G, (22) Z~ Z =[~1⊤;…;~N⊤]. =[ z_1 ;…; z_N ]. (23) To stabilize training and preserve the per-site total mass (row-wise ℓ1 _1 sum), we first take a convex residual update with a mixing coefficient λ∈[0,1]λ\!∈\![0,1] and then renormalize each row: Z Z =(1−λ)Z+λZ~, =(1-λ)\,Z\;+\;λ\, Z, (24) Zj,:′ Z _j,: =∑g=1GZj,g∑g=1GZ^j,g+εZ^j,:,j=1,…,N, = _g=1^GZ_j,g _g=1^G Z_j,g+ \; Z_j,:\,, j=1,…,N, (25) where ε>0 \!>\!0 avoids division by zero. Eq.( 24) provides a conservative blend between the old and updated coefficients, while Eq.( 25) rescales each site’s coefficients so that ∑gZj,g′=∑gZj,g _gZ _j,g= _gZ_j,g. Complexity. Per head, forming the window logits requires projecting ξj∈ℝG(2d+2) _j ^G(2d+2) to ℝDR^D, costing (BHN⋅G(2d+2)⋅D)O(BHN· G(2d+2)· D). The two transfers N↔GN\! \!G cost (BHNGD)O(BHNGD), and modal attention costs (BHG2D)O(BHG^2D) (independent of N). Overall, for fixed G the dominant term is linear in N, with explicit dependence on spatial dimension d. Expressivity of the modal operator. We formalize that PG Gaussian Attention can approximate a broad class of continuous operators: Theorem 2.2 (Universal approximation in modal form). Let :Lp(Ω;ℝcin)→Lq(Ω;ℝcout)T:L^p( ;R^c_in)\!→\!L^q( ;R^c_out) be continuous on bounded sets and admit either a Mercer/Hilbert–Schmidt kernel K(,′)K(x,x ) or a low-rank factorization ≈Φ(⋅)Φ(⋅)⊤T\!≈\! (·)\,K\, (·)^\! with continuous Φ:Ω→ℝm : \!→\!R^m. Then, for any ε>0 >0, there exist a modal budget G and parameters Θ of our encoder, Gaussian modal windows, PG Gaussian Attention, and decoder such that ‖Θ−‖<ε\|G_ -T\|< (operator norm on bounded subsets). (Proof in Appx.B.2) 2.3 Gaussian Particle Operator: Overall Framework 2.3.1 Neural Operator Formulation Let Ω⊂ℝd \!⊂\!R^d be the domain, a:Ω→ℝcina: \!→\!R^c_in the input field, and u:Ω→ℝcoutu: \!→\!R^c_out the target field. We model the map a↦ua u by a neural operator Θ=fϕdec∘(ψ)∘n∘(⋅;Πθ(⋅))∘Eθ,G_ \;=\;f^dec_φ\ \ (O_ψ )^\! n\ \ Z (\,·\;;\, _θ(·) )\ \ E_θ, (26) where: • EθE_θ (encoder) predicts Gaussian particles Πθ(a)=(μθ,σθ,wθ) _θ(a)= ( _θ, _θ,w_θ ) on the context locations xjj=1N\x_j\_j=1^N (using both aja_j and xjx_j); • (⋅;Π)Z(·; ) evaluates the Gaussian basis and returns per-location, G-dimensional coefficients Z∈ℝN×GZ ^N× G with zj,i=wj,iexp(−12‖(j−μj,i)/σj,i‖22);z_j,i\;=\;w_j,i\, \! (- 12 \|(x_j- _j,i)/ _j,i \|_2^2 ); (27) • ψO_ψ is the Gaussian-basis operator (Sec. 2.2) acting on Z and parameterized by PG Gaussian Attention; it can be applied n times: Z(0)=Z, Z^(0)=Z, (28) Z(k+1)=ψ(Z(k);Πθ(a)), Z^(k+1)=O_ψ (Z^(k); _θ(a) ), (29) • fϕdecf^dec_φ (decoder) maps the updated basis to the output field values: u^(j)=fϕdec(Zj,:(n))\ u(x_j)=f^dec_φ\! (Z^(n)_j,: ). The construction is resolution-agnostic: for any query set x~mm=1M\ x_m\_m=1^M (2D/3D grids or irregular meshes), we keep the particle field Πθ(a) _θ(a) predicted on the context locations fixed and simply re-evaluate zj,i(x~m)=wj,iexp(−12‖(x~m−μj,i)/σj,i‖22)z_j,i( x_m)=w_j,i\, \! (- 12\|( x_m- _j,i)/ _j,i\|_2^2 ), followed by the same ψO_ψ and fϕdecf^dec_φ. 2.3.2 Pipeline overview As shown in Fig. 2, given an input field a()a(x), the encoder EθE_θ produces per-site Gaussian particles Πθ(a)=(μ,σ,w) _θ(a)=(μ,σ,w), i.e., a Gaussian field. At query locations jj=1N\x_j\_j=1^N, we then evaluate the Gaussian basis by Eq.(27) to obtain Z∈ℝN×GZ ^N× G. The Gaussian-basis operator ψO_ψ acts on Z in modal space and can be applied for n stages as in Eq.(28) to capture multi-step coupling, yielding Z(n)Z^(n). Finally, the decoder fϕdecf^dec_φ maps Zj,:(n)Z^(n)_j,: to u^(j) u(x_j). During training, the target u()u(x) may also be encoded by EθE_θ to provide an auxiliary Gaussian-field supervision signal. 3 Experiments Figure 3: Interpretable visualization on an in-distribution sample (above) and an out-of-distribution sample (below). Left to right: ground truth, reconstruction from the Gaussian basis, absolute error, and learned Gaussian particles overlaid (ellipses: center μ, axes ∝σ σ, color/size ∝w w). Benchmarks. We evaluate on two Navier–Stokes surrogates, three simulated PDE systems, and two real reanalyses to span 2D→ 3D, planar→ , and regular→ domains. NS2D (Kovachki et al., 2023) is an incompressible periodic box sampled on 64×6464×64; NS3D (Takamoto et al., 2022) extends to a periodic cube on 64364^3, stressing 3D scalability. Airfoil (Li et al., 2023) considers 2D airfoil flow on irregular geometries (e.g., NACA profiles) to test mesh-/shape-generalization. From The Well (Ohana et al., 2024), PlanetSWE is a forced rotating shallow-water system on the sphere (equiangular 256×512256×512), and turbulent_radiative_layer_2D is a Kelvin–Helmholtz-driven turbulent mixing layer with radiative cooling on a Cartesian 384×128384×128 grid. ERA5 (Hersbach et al., 2023) uses one month on the 0.25∘0.25 global grid (721×1440721×1440), with variables 2 m temperature (t) and 10 m zonal wind (u). CARRA (Schyberg et al., 2020) uses one Arctic month on its native regional grid (989×789989×789) with an irregular land/sea/ice mask, and we retain only 10 m meridional wind (v10v_10). We train one-step operators and assess multi-step rollouts on temporal datasets, using native grids/meshes, and apply latitude/area weighting on spherical domains (ERA5, PlanetSWE) and the native mask for CARRA. Baselines. We compare against a broad set of neural operators spanning spectral, convolutional, multiresolution, graph-based, and attention-based designs. We include physics-inspired / structure-aware operators: (Geo-)FNO (Li et al., 2021, 2023) (Fourier neural operator on regular grids and its geometry-aware variant), CNO (Raonic et al., 2023) (convolutional neural operator), LSM (Wu et al., 2023) (learned spectral mixing), M2NO (Li et al., 2025a) (multiresolution neural operator), and AMG (Li et al., 2025b) (multi-graph neural operator for arbitrary geometries). We also compare with Transformer-based operators that parameterize the operator via attention (global coupling): Galerkin Transformer (Cao, 2021), GNOT (Hao et al., 2023), ONO (Xiao et al., 2024), LNO (Wang and Wang, 2024), and Transolver (Wu et al., 2024). All baselines use identical data splits, losses, and rollout protocols, and we keep model capacity comparable; 2D/3D/irregular variants are used where applicable. Implementations. We evaluate all models using the relative L2L_2 error on held-out sets. Inputs/targets are normalized per variable using training statistics; models are trained on the normalized data, and all metrics are computed after inverse normalization. Training uses AdamW (Loshchilov and Hutter, 2019) with an initial learning rate of 10−310^-3 and a StepLR scheduler that reduces the learning rate at fixed intervals. All experiments are run on a single NVIDIA RTX 4090 GPU. Per-dataset and per-baseline hyperparameters are provided in Appx.C. Figure 4: Layer-wise evolution of the Gaussian particle field. Particle activations after successive PG Gaussian Attention layers. 3.1 Benchmark Performance Table 1 reports L2L_2 errors across synthetic Navier–Stokes surrogates, real reanalyses, and additional PDE systems spanning irregular geometries (Airfoil), Cartesian turbulence with radiative effects, and spherical dynamics (PlanetSWE). Overall, the results are consistent with the inductive biases of each model family while demonstrating that GPO is a strong and stable performer across all regimes. (i) Regular-grid fluid surrogates (NS2D/NS3D). Spectral or multiresolution priors are advantageous on clean periodic grids: Fourier- and multiscale-based baselines are consistently strong, reflecting their efficiency at representing globally coherent modes. GPO matches or surpasses these methods on both 2D and 3D Navier–Stokes, indicating that its Gaussian particle basis can recover the benefits of modal representations without committing to a fixed Fourier parameterization. This supports our motivation that learned localized atoms coupled through a Petrov–Galerkin measurement can serve as an effective alternative to hand-crafted spectral kernels, while remaining applicable beyond periodic grids. (i) Large-scale geophysical reanalyses (ERA5/CARRA). On global and masked regional grids, performance favors models that combine (a) non-Euclidean or mask-aware handling and (b) sufficient long-range coupling. Purely spectral designs degrade under spherical weighting and irregular masks, while attention-style operators improve by modeling global dependencies but may incur unfavorable scaling. GPO achieves the best overall accuracy on these reanalyses, consistent with our contributions: the Gaussian particle representation provides localized, interpretable support that adapts naturally to masked domains, and PG Gaussian Attention enables global interaction in a modal space, reducing sensitivity to grid irregularity and avoiding quadratic dependence on the number of spatial samples. (i) Geometry generalization and heterogeneous physics (Airfoil / Turbulent / PlanetSWE). The additional benchmarks highlight transfer across geometry and operator classes. On Airfoil, which stresses irregular boundaries and shape variation, GPO outperforms geometry-aware spectral and attention baselines, suggesting that particle-based latent states are well-suited to representing boundary-localized structures. On Turbulent, where sharp gradients and multi-scale mixing are dominant, GPO is consistently better than both convolutional and transformer baselines, aligning with our claim that localized Gaussian atoms capture high-frequency structures more effectively than globally parameterized kernels. On PlanetSWE, which requires global coupling on the sphere, GPO also leads, indicating that PG-attention provides sufficient long-range interaction while maintaining robustness across coordinate systems. Summary. Across all benchmarks, GPO is competitive on regular synthetic settings and shows clear advantages on irregular, masked, and non-Cartesian domains. This empirical pattern aligns with our central thesis: representing PDE states as compact, localized, and interpretable Gaussian particles and evolving them with modal Petrov–Galerkin attention yields a general-purpose operator that remains accurate under changing geometry, resolution, and physics. 3.2 Interpretable Diagnostics 3.2.1 Gaussian Reconstruction Setup. Before evaluating the operator, we validate that the learned Gaussian particle basis (Sec. 2.1) is faithful and interpretable. We train the encoder–decoder to reconstruct the field using weighted Gaussian evaluation, together with particle regularizers (center alignment and a scale-range barrier). For both in-distribution (ID) and out-of-distribution (OOD) samples, we visualize the ground truth, reconstruction, absolute error, and particle overlays (ellipses: center =μ=μ, axes ∝σ σ, color/size ∝w w). Observations. Fig. 3 shows that (i) particles consistently concentrate on coherent structures (fronts, filaments, vortices), with anisotropic footprints aligned to local flow directions; (i) reconstruction errors remain localized near sharp gradients and subgrid-scale filaments; and (i) under OOD shifts, particle geometry and weights adapt smoothly, preserving large- and meso-scale structures with only moderate degradation. This supports that Gaussian particles form a visualizable, mesh-agnostic state representation that remains stable beyond the training distribution, and provides an explicit trial space for subsequent operator updates. 3.2.2 Layer-wise Particle Dynamics Setup. During prediction, we apply the Gaussian Particle Operator (Sec. 2.2) for n stages. The encoder-fixed particle geometry (μ,σ)(μ,σ) defines trial atoms, while PG Gaussian Attention updates per-site coefficients zj(k)∈ℝGz^(k)_j ^G at stage k. We visualize layer-wise particle fields by overlaying particle footprints and plotting a scalar activation summary (e.g., A(k)(xj)=∑g=1Gzj,g(k)A^(k)(x_j)= _g=1^Gz^(k)_j,g), making the redistribution of modal energy directly observable. Interpretation. In Fig. 4, the evolution of activation maps admits a mechanistic reading: (i) local homogenization resembles diffusion among nearby modes; (i) directional transfers reflect advection-like transport aligned with anisotropic particle geometry; and (i) localized growth/decay indicates cross-scale energy exchange. Comparing the input a(x)a(x) and target u(x)u(x), the progressive adjustment from Layer 1→ n tracks the emergence and transport of the fine-scale structures that differentiate u from a. This supports our second contribution: PG-attention yields interpretable global coupling operating in a compact modal space. 2244668810101212141416161818202010−1.510^-1.510−1.410^-1.410−1.310^-1.310−1.210^-1.210−1.110^-1.110−110^-1Rollout step ttRelative L2L_2 errorLSMFNOTransolverGPO Figure 5: Rollout stability schematic: error versus rollout horizon. 3.2.3 Rollout Stability Setup. To probe long-horizon behavior, we perform autoregressive rollouts and report the relative L2L_2 error as a function of rollout step t (Figure 5). We compare representative spectral and attention baselines against GPO under identical rollout protocols. Observations. All methods exhibit error accumulation over time, but GPO shows consistently slower growth and reduced drift across the horizon. This behavior is consistent with our design: the particle geometry anchors the representation to localized structures, while PG measurement and modal coupling mitigate compounding aliasing and stabilize global interactions without directly attending over all spatial points. Practically, this indicates that the learned dynamics remain coherent under repeated composition, which is critical for real forecasting workloads. 3.2.4 Spectral Fidelity Setup. Beyond pointwise errors, we assess physical plausibility via the kinetic energy spectrum E(k)E(k) on NS2D (Fig. 6), computed from Fourier modes and shown on a log–log scale with a reference inertial-range slope. Observations. Compared to baselines, GPO tracks the ground-truth spectrum more closely across both low and intermediate wavenumbers, and exhibits a delayed high-k roll-off, indicating less over-smoothing of small-scale content. In particular, spectral baselines tend to lose energy early in the high-frequency range, while attention models partially recover mid-range structure but still attenuate high-k modes. The improved spectral alignment of GPO supports our central thesis: representing states with localized Gaussian particles and evolving them through PG-attention better preserves multi-scale structure, yielding dynamics that are not only accurate but also physically consistent. Additional visualizations are provided in Appx. E. 3.3 Model Analysis Ablations (Table 3) confirm that the Gaussian Field is essential, while even a single Gaussian per site improves over MLP baselines. The best results arise from combining the Gaussian Field with PG Gaussian Attention. Increasing the modal budget G generally helps but exhibits diminishing returns; we therefore adopt a moderate G for a balanced accuracy–efficiency trade-off. Complexity (Table 4) shows that GPO maintains low memory use and competitive runtime, scaling near-linearly with the number of query points (unlike spatial self-attention), and offering favorable trade-offs across 2D/3D and irregular domains. See Appx.D for full ablations and analyses. 10010^010110^110210^210−910^-910−810^-810−710^-710−610^-610−510^-510−410^-410−310^-310−210^-210−110^-1Wavenumber kkEnergy spectrum E(k)E(k)GTFNOTransolverGPORef. slope k−5/3k^-5/3 Figure 6: NS2D energy spectrum (log–log). 4 Conclusion We introduced the Gaussian Particle Operator (GPO), a resolution-agnostic neural operator that represents fields with an interpretable Gaussian particle basis and performs basis-to-basis updates via Petrov–Galerkin Gaussian Attention. The pipeline makes intermediate objects—particles (μ,σ,w)(μ,σ,w), modal windows, and inter-modal couplings—directly visualizable, while remaining robust across regular grids, irregular masks, and 2D→ 3D settings. Empirically, GPO is consistently strong on large-scale geophysical datasets, with improved spectral behavior and stable rollouts. Limitations and future work Scalability. Accuracy and cost depend on the modal budget G and head width. While the N↔GN\! \!G transfers are linear in N, storing and applying NGNG windows can become a bottleneck at extreme resolutions. Promising directions include adaptive particles, sparse routing/pruning, structured or low-rank coupling in modal space. Physics integration. Training is largely data-driven with lightweight regularization, and does not enforce invariants or hard constraints (e.g., divergence-free flow or boundary conditions). Future work will integrate PDE structure into the particle basis and PG updates, aiming for operators that are both accurate and more mechanistically interpretable. References A. N. Brooks and T. J.R. Hughes (1982) Streamline upwind/petrov-galerkin formulations for convection dominated flows with particular emphasis on the incompressible navier-stokes equations. Computer Methods in Applied Mechanics and Engineering 32 (1), p. 199–259. External Links: ISSN 0045-7825, Document Cited by: §2.2. M. D. Buhmann (2000) Radial basis functions. Acta numerica 9, p. 1–38. Cited by: §1. S. Cao (2021) Choose a transformer: fourier or galerkin. In NeurIPS, p. 24924–24940. Cited by: Appendix A, §1, §3. L. P. Franca, G. Hauke, and A. Masud (2006) Revisiting stabilized finite element methods for the advective–diffusive equation. Computer Methods in Applied Mechanics and Engineering 195 (13), p. 1560–1572. External Links: ISSN 0045-7825, Document, Link Cited by: §2.2. G. Gupta, X. Xiao, and P. Bogdan (2021) Multiwavelet-based operator learning for differential equations. In NeurIPS, p. 24048–24062. Cited by: Appendix A, §1. M. E. Gurtin (1982) An introduction to continuum mechanics, mathematics in science and engineering. Academic Press, Cambridge. Cited by: §1. Z. Hao, C. Ying, Z. Wang, H. Su, Y. Dong, S. Liu, Z. Cheng, J. Zhu, and J. Song (2023) GNOT: a general neural operator transformer for operator learning. ArXiv abs/2302.14376. Cited by: Appendix A, §1, §3. J. He, X. Liu, and J. Xu (2024) MgNO: efficient parameterization of linear operators via multigrid. ArXiv abs/2310.19809. Cited by: Appendix A. H. Hersbach, B. Bell, P. Berrisford, G. Biavati, A. Horányi, J. Muñoz Sabater, J. Nicolas, C. Peubey, R. Radu, I. Rozum, D. Schepers, A. Simmons, C. Soci, D. Dee, and J.-N. Thépaut (2023) ERA5 hourly data on single levels from 1940 to present. Copernicus Climate Change Service (C3S) Climate Data Store (CDS). External Links: Document, Link Cited by: §3. C. Johnson (2012) Numerical solution of partial differential equations by the finite element method. Courier Corporation, North Chelmsford. Cited by: §1. G. A. Klaasen and W. C. Troy (1984) Stationary wave solutions of a system of reaction-diffusion equations derived from the fitzhugh-nagumo equations. SIAM Journal on Applied Mathematics 44 (1), p. 96–110. External Links: ISSN 00361399, Link Cited by: §1. N. B. Kovachki, Z. Li, B. Liu, K. Azizzadenesheli, K. Bhattacharya, A. M. Stuart, and A. Anandkumar (2023) Neural operator: learning maps between function spaces with applications to pdes. J. Mach. Learn. Res. 24, p. 89:1–89:97. Cited by: Appendix A, §1, §3. Z. Li, Z. Lai, X. Zhang, and W. Wang (2025a) M2NO: an efficient multi-resolution operator framework for dynamic multi-scale pde solvers. External Links: 2406.04822, Link Cited by: Appendix A, §1, §3. Z. Li, H. Song, D. Xiao, Z. Lai, and W. Wang (2025b) Harnessing scale and physics: A multi-graph neural operator framework for pdes on arbitrary geometries. In KDD (1), p. 729–740. Cited by: Appendix A, §1, §3. Z. Li, N. B. Kovachki, K. Azizzadenesheli, B. Liu, K. Bhattacharya, A. M. Stuart, and A. Anandkumar (2020a) Neural operator: graph kernel network for partial differential equations. ArXiv abs/2003.03485. Cited by: Appendix A. Z. Li, D. Z. Huang, B. Liu, and A. Anandkumar (2023) Fourier neural operator with learned deformations for pdes on general geometries. J. Mach. Learn. Res. 24, p. 388:1–388:26. Cited by: §3, §3. Z. Li, N. B. Kovachki, K. Azizzadenesheli, B. Liu, A. M. Stuart, K. Bhattacharya, and A. Anandkumar (2020b) Multipole graph neural operator for parametric partial differential equations. In NeurIPS, Cited by: Appendix A. Z. Li, N. B. Kovachki, K. Azizzadenesheli, B. Liu, K. Bhattacharya, A. M. Stuart, and A. Anandkumar (2021) Fourier neural operator for parametric partial differential equations. In ICLR, Cited by: Appendix A, §1, §3. I. Loshchilov and F. Hutter (2019) Decoupled weight decay regularization. In ICLR (Poster), Cited by: §3. L. Lu, P. Jin, G. Pang, Z. Zhang, and G. E. Karniadakis (2021) Learning nonlinear operators via deeponet based on the universal approximation theorem of operators. Nat. Mach. Intell. 3 (3), p. 218–229. Cited by: Appendix A, §1. R. McKeown, R. Ostilla-Mónico, A. Pumir, M. P. Brenner, and S. M. Rubinstein (2020) Turbulence generation through an iterative cascade of the elliptical instability. Science Advances 6 (9), p. eaaz2717. External Links: Document, Link Cited by: §1. R. Ohana, M. McCabe, L. Meyer, R. Morel, F. J. Agocs, M. Beneitez, M. Berger, B. Burkhart, S. B. Dalziel, D. B. Fielding, D. Fortunato, J. A. Goldberg, K. Hirashima, Y. Jiang, R. R. Kerswell, S. Maddu, J. Miller, P. Mukhopadhyay, S. S. Nixon, J. Shen, R. Watteaux, B. R. Blancard, F. Rozet, L. H. Parker, M. D. Cranmer, and S. Ho (2024) The well: a large-scale collection of diverse physics simulations for machine learning. In NeurIPS, Cited by: §3. J. Park and I. W. Sandberg (1991) Universal approximation using radial-basis-function networks. Neural Computation 3 (2), p. 246–257. External Links: ISSN 0899-7667, Document, Link Cited by: §1. B. Raonic, R. Molinaro, T. D. Ryck, T. Rohner, F. Bartolucci, R. Alaifari, S. Mishra, and E. de Bézenac (2023) Convolutional neural operators for robust and accurate learning of pdes. In NeurIPS, Cited by: §3. H. Schyberg, X. Yang, M. A. Ø. Køltzow, B. Amstrup, Å. Bakketun, E. Bazile, J. Bojarova, J. E. Box, P. Dahlgren, S. Hagelin, M. Homleid, A. Horányi, J. Høyer, Å. Johansson, M. A. Killie, H. Körnich, P. Le Moigne, M. Lindskog, T. Manninen, P. Nielsen Englyst, K. P. Nielsen, E. Olsson, B. Palmason, C. Peralta Aros, R. Randriamampianina, P. Samuelsson, R. Stappers, E. Støylen, S. Thorsteinsson, T. Valkonen, and Z. Q. Wang (2020) Arctic regional reanalysis on single levels from 1991 to present. Copernicus Climate Change Service (C3S) Climate Data Store (CDS). External Links: Document, Link Cited by: §3. L. Serrano, L. L. Boudec, A. K. Koupaï, T. X. Wang, Y. Yin, J. Vittaut, and P. Gallinari (2023) Operator learning with neural fields: tackling pdes on general geometries. In NeurIPS, Cited by: §1. L. Serrano, T. X. Wang, E. L. Naour, J. Vittaut, and P. Gallinari (2024) AROMA: preserving spatial structure for latent PDE modeling with local neural fields. In NeurIPS, Cited by: §1. M. F. Shlesinger, B. J. West, and J. Klafter (1987) Lévy dynamics of enhanced diffusion: application to turbulence. Phys. Rev. Lett. 58, p. 1100–1103. External Links: Document, Link Cited by: §1. M. Takamoto, T. Praditia, R. Leiteritz, D. MacKinlay, F. Alesiani, D. Pflüger, and M. Niepert (2022) PDEBench: an extensive benchmark for scientific machine learning. In NeurIPS, Cited by: §3. T. Wang and C. Wang (2024) Latent neural operator for solving forward and inverse PDE problems. In NeurIPS, Cited by: §3. A. Wazwaz (2002) Partial differential equations: methods and applications. Balkema Publishers, Leiden. Cited by: §1. H. Wu, T. Hu, H. Luo, J. Wang, and M. Long (2023) Solving high-dimensional pdes with latent spectral models. In ICML, Proceedings of Machine Learning Research, Vol. 202, p. 37417–37438. Cited by: §3. H. Wu, H. Luo, H. Wang, J. Wang, and M. Long (2024) Transolver: a fast transformer solver for pdes on general geometries. ArXiv abs/2402.02366. Cited by: Appendix A, §C.1, §3. Z. Xiao, Z. Hao, B. Lin, Z. Deng, and H. Su (2024) Improved Operator Learning by Orthogonal Attention. ArXiv abs/2310.12487. Cited by: Appendix A, §3. Appendix A Related Work Neural operators. Classical neural operator methods aim to learn mappings between function spaces directly from data, typically by parameterizing a resolution-agnostic kernel or lifting to a latent space and learning integral transforms. Representative approaches include the Fourier Neural Operator (FNO), which performs global convolution via spectral multipliers to approximate operator kernels in Fourier space (Li et al., 2021; Kovachki et al., 2023), and DeepONet, which decomposes an operator into branch/trunk networks to separately encode input functions and query coordinates (Lu et al., 2021). Variants extend these ideas with multiresolution bases (Li et al., 2020b; Gupta et al., 2021; He et al., 2024; Li et al., 2025a), graph or kernelized message passing (Li et al., 2025b), and learned Green’s functions (Li et al., 2020a). These models are purely data driven: while many designs are physics-inspired, their internal representations are typically opaque. In particular, the learned latent bases and mixing weights are not tied to interpretable physical primitives, which limits diagnostic insight and the ability to attribute predictions to physically meaningful components. Transformer-based methods. A recent line of work adopts Transformers to parameterize neural operators, replacing hand-crafted kernel parameterizations with data-driven attention. Examples include Galerkin Transformers, which align attention with variational forms (Cao, 2021), GNOT (Hao et al., 2023) that leverage attention for long-range coupling, Transolver (Wu et al., 2024), which introduces slice-based attention for efficient global mixing, and operator networks that stack attention with physics-informed objectives (Xiao et al., 2024). Empirically, with sufficiently large training corpora and careful scaling, attention-based operators often match or surpass traditional neural operators in expressive power and generalization to out-of-distribution forcings and grids. However, these gains come with two well-known limitations. (i) Lack of interpretability: standard attention weights are not anchored to physically interpretable trial/test functions, making it difficult to ascribe predictions to identifiable modes or localized mechanisms. (i) Frequency bias: global self-attention tends to emphasize low-rank, global correlations (low-frequency structure), while recovering sharp, localized, or high-frequency phenomena often requires architectural add-ons or extensive data augmentation. As a result, pure Transformer operators provide limited physical attribution and may under-represent fine-scale features without additional inductive biases. Appendix B Expressivity of the Gaussian Field and GPO B.1 Expressivity of the Gaussian Field Lemma B.1 (Density of Gaussian mixtures). On compact Ω , finite mixtures of anisotropic Gaussians are dense in C(Ω)C( ) (and dense in Lr(Ω)L^r( ) for 1≤r<∞1\!≤\!r\!<\!∞). Hence for any continuous scalar field v and ε>0 >0, there exist G,μi,σi,wii=1G,\ _i, _i,w_i\_i=1^G such that ∥v(⋅)−∑i=1Gwiexp(−12∥(⋅−μi)/σi∥2)∥∞<ε\|v(·)- _i=1^Gw_i (- 12\|(·- _i)/ _i\|^2)\|_∞< . Sketch. Standard universal approximation results for radial basis functions/Gaussian kernels. Proof. We give a constructive proof based on Gaussian mollification and Riemann sums. Step 1: Approximate identity via Gaussian mollifiers. Let Ω⊂ℝd ^d be compact and let v∈C(Ω)v∈ C( ). By Tietze’s extension theorem there exists v~∈Cc(ℝd) v∈ C_c(R^d) such that v~|Ω=v v|_ =v. For Σ∈ℝd×d ^d× d symmetric positive definite, set the (unnormalized) Gaussian ϕΣ(x)=exp(−12x⊤Σ−1x). _ (x)\;=\; \! (- 12\,x ^-1x ). Let Σϵ↓0\ _ε\_ε 0 be any family with ‖Σϵ‖→0\| _ε\|→ 0. Since Gaussians form an approximate identity, the normalized mollification v~∗ϕΣϵ/∫ℝdϕΣϵ v* _ _ε/ _R^d _ _ε converges to v~ v uniformly on compact sets as ϵ↓0ε 0 (uniform continuity of v~ v and standard approximate-identity properties). Because the normalization constant is a positive scalar depending only on Σϵ _ε, we can absorb it into the mixture weights later. Hence, for any η>0η>0 there exists ϵ0 _0 such that for all 0<ϵ≤ϵ00<ε≤ _0, supx∈Ω|(v~∗ϕΣϵ)(x)−v~(x)|<η2. _x∈ |\,( v* _ _ε)(x)- v(x)\, |< η2. (30) Step 2: Riemann-sum approximation of the convolution (finite mixture). Fix such an ϵε, write Σ=Σϵ = _ε, and denote the convolution (v~∗ϕΣ)(x)=∫ℝdv~(y)ϕΣ(x−y)y.( v* _ )(x)\;=\; _R^d v(y)\, _ (x-y)\,dy. Since v~ v is compactly supported and continuous while ϕΣ _ is continuous and rapidly decaying, the integrand is continuous with compact support in y uniformly in x∈Ωx∈ . Hence Riemann sums approximate the integral uniformly in x: there exists a finite set of nodes μii=1G⊂ℝd\ _i\_i=1^G ^d with associated positive quadrature weights Δii=1G\ _i\_i=1^G such that supx∈Ω|(v~∗ϕΣ)(x)−∑i=1Gv~(μi)ϕΣ(x−μi)Δi|<η2. _x∈ |( v* _ )(x)\;-\; _i=1^G v( _i)\, _ (x- _i)\, _i |< η2. (31) Define mixture weights wi:=v~(μi)Δiw_i:= v( _i)\, _i (real-valued; the lemma does not restrict their sign), and note that each term is exactly a (shared-covariance) Gaussian atom exp(−12‖(x−μi)‖Σ−12) \! (- 12\|(x- _i)\|_ ^-1^2 ), i.e., ∑i=1GwiϕΣ(x−μi)=∑i=1Gwiexp(−12(x−μi)⊤Σ−1(x−μi)). _i=1^Gw_i\, _ (x- _i)\;=\; _i=1^Gw_i \! (- 12(x- _i) ^-1(x- _i) ). Step 3: Uniform approximation on Ω . Combining (30) and (31), supx∈Ω|v~(x)−∑i=1Gwiexp(−12(x−μi)⊤Σ−1(x−μi))| _x∈ | v(x)- _i=1^Gw_i \! (- 12(x- _i) ^-1(x- _i) ) | ≤ ≤ supx∈Ω|v~(x)−(v~∗ϕΣ)(x)|+supx∈Ω|(v~∗ϕΣ)(x)−∑iwiϕΣ(x−μi)| _x∈ | v(x)-( v* _ )(x) |+ _x∈ |( v* _ )(x)- _iw_i _ (x- _i) | < < η. η. Restricting back to Ω (where v~=v v=v) yields ∥v(⋅)−∑i=1Gwiexp(−12∥(⋅−μi)∥Σ−12)∥∞<η. \|v(·)- _i=1^Gw_i \! (- 12\|(·- _i)\|_ ^-1^2 ) \|_∞<η. Since η>0η>0 was arbitrary, finite mixtures of (possibly anisotropic) Gaussians are dense in C(Ω)C( ). Anisotropy and vector-valued extension. We used a common covariance Σ for clarity; allowing mode-dependent Σi _i only increases expressivity, so the same result holds with per-atom anisotropy. For vector-valued v, apply the scalar result componentwise. LrL^r density. Because C(Ω)C( ) is dense in Lr(Ω)L^r( ) for 1≤r<∞1≤ r<∞ on compact Ω , the uniform approximation implies LrL^r approximation, completing the proof. ∎ B.2 Expressivity of GPO Theorem B.2 (Universal approximation in modal form). Let :→T:X be a continuous operator on compacta that admits a Hilbert–Schmidt (Mercer-type) kernel K(,′)K(x,x ) or, more generally, a low-rank factorization ≈Φ(⋅)Φ(⋅)⊤T≈ (·)\,K\, (·) with continuous features Φ:Ω→ℝm : ^m. Then, for any ε>0 >0, there exist G and network parameters Θ such that ‖Θ−‖→<ε\|G_ -T\|_X < . Sketch. By Lemma B.1 and universal approximation of MLPs, windows p(,g)p(x,g) and latent features S(Z())S(Z(x)) approximate Φ() (x); attention realizes a trainable K on the G modes. The scatter-and-decoder emulate the output feature map. Increasing G and widths yields density in the space of continuous operators. Proof. We prove the claim for operators on compact domains by reducing to a finite–rank Mercer approximation and showing that each stage of our pipeline can approximate the corresponding finite–dimensional objects arbitrarily well. Throughout, ∥⋅∥→\|·\|_X denotes the operator norm on bounded subsets. Step 0: Mercer (or low–rank) truncation. Assume T is continuous on bounded sets and admits either a Hilbert–Schmidt kernel K(,′)K(x,x ) or, more generally, a low–rank factorization ≈Φ(⋅)Φ(⋅)⊤T≈ (·)\,K\, (·) with continuous Φ:Ω→ℝm : ^m. In the Mercer case, by spectral theory, K(,′)=∑r=1∞λrφr()φr(′),λr≥0,φr⊂C(Ω),K(x,x )= _r=1^∞ _r\, _r(x)\, _r(x ), _r≥ 0,\;\ _r\⊂ C( ), and the partial sums define finite–rank operators mf()=∑r=1mλrφr()∫φr(′)f(′)′T_mf(x)= _r=1^m _r _r(x) _r(x )f(x )\,dx with ‖−m‖→0\|T-T_m\|→ 0 as m→∞m→∞ (uniform on compacta). In the given low–rank form, select m and continuous Φm:Ω→ℝm _m: ^m, m∈ℝm×mK_m ^m× m such that ‖−m‖<ε/3,mf()=Φm()m∫ΩΦm(′)⊤f(′)′. \|T-T_m \|< /3, _mf(x)= _m(x)\,K_m _ _m(x ) f(x )\,dx . (32) Step 1: Approximating the feature maps by Gaussian basis + MLPs. By Lemma B.1 (density of Gaussian mixtures) and universal approximation of MLPs, for any δ>0δ>0 there exist: (i) a pointwise encoder/evaluator producing Z()∈ℝGZ(x) ^G from Gaussian particles (μ,σ,w)(μ,σ,w) and a small MLP S such that the trial features Ψ()∈ℝD (x) ^D, defined by Ψ()=S(Z()) (x)=S(Z(x)), satisfy sup∈Ω‖Ψ()−Φm()‖2<δ; _x∈ \| (x)- _m(x) \|_2<δ; (33) (i) head–wise Gaussian modal windows p(,g)≥0p(x,g)≥ 0 with ∑g=1Gp(,g)=1 _g=1^Gp(x,g)=1, implemented by linear maps on [Z(),(μ,σ,w)()][Z(x),(μ,σ,w)(x)] and a softmax, such that the test functionals ℳg(f)=∫Ωp(,g)f()M_g(f)= _ p(x,g)\,f(x)\,dx approximate the m target coordinates ∫Φm()⊤f() _m(x) f(x)\,dx after a fixed linear readout. Concretely, there exists W∈ℝm×GW ^m× G with ‖W[ℳg(⋅)]g=1G−∫Φm(⋅)⊤(⋅)‖<C1δ\|\,W\,[M_g(·)]_g=1^G- _m(·) (·)\,\|<C_1δ. (One can view Wp(,⋅)Wp(x,·) as a learned quadrature/test family for the m coordinates.) Step 2: Discrete PG measurement and quadrature error. Given a discretization jj=1N\x_j\_j=1^N with empirical measure converging to the sampling measure on Ω , the N→GN\!→\!G aggregation used in Sec. 2.2 forms tokens tg=∑j=1Np(j,g)Ψ(j)∑j=1Np(j,g)∈ℝD.t_g\;=\; _j=1^Np(x_j,g)\, (x_j) _j=1^Np(x_j,g) ^D. By uniform continuity of Ψ and p(⋅,g)p(·,g) on compact Ω , Riemann (or Monte Carlo) sums converge to the integrals. Hence there exists N0N_0 so that for all N≥N0N≥ N_0, ‖[tg]g=1G−[∫p(,g)Ψ()∫p(,g)]g=1G‖<C2δ. \| [\;t_g\; ]_g=1^G\;-\; [\; p(x,g)\, (x)\,dx p(x,g)\,dx\; ]_g=1^G \|<C_2δ. (34) Post-multiplying by W and using (33) shows that the vector of m measured coordinates is within C3δC_3δ of ∫Φm()⊤f() _m(x) f(x)\,dx for any f in a bounded set. Step 3: Implementing the modal coupling by attention + linear maps. We next show that the G×G× G modal attention stage can realize the finite linear map mK_m (up to basis changes) to arbitrary precision. Using the head projections WzW_z and WoutW_out, the attention block computes T~=α(TWV),Y=(T~)Wout, T\;=\;α\, (TW_V ), Y\;=\; ( T )W_out, where T∈ℝG×DT ^G× D stacks the tokens tgt_g, α is the softmax attention matrix, and WV,WoutW_V,W_out are learned linear maps. Since softmax can approximate a Kronecker–delta (by sending on–diagonal logits to +∞+∞ and off–diagonal to −∞-∞), we can set α≈IGα≈ I_G arbitrarily closely. Then Y≈T(WVWout)Y≈ T(W_VW_out). Because WV,WoutW_V,W_out are unconstrained, their product can approximate any target matrix M∈ℝD×mM ^D× m to arbitrary precision. Choosing M to implement the composition WmW\,K_m (after the measurement map from Step 2), we obtain a block that emulates v↦mv _mv in the m–dimensional modal coordinates. (If desired, one may keep α nontrivial and absorb its effect into the surrounding linear maps; the argument is unchanged.) Step 4: Scatter and pointwise decoding. The G→NG\!→\!N scatter re-distributes the mixed modal features back to locations via the same windows p(,g)p(x,g), followed by a pointwise decoder MLP fϕdec:ℝG→ℝcoutf^dec_φ:R^G ^c_out. Since MLPs are universal approximators on compacta, the composition can approximate the desired output feature map ↦Φm()x _m(x) (or its linear image) uniformly, matching the form in (32). Step 5: Error aggregation. Let εm=‖−m‖<ε/3 _m=\|T-T_m\|< /3 be the truncation error. Pick δ>0δ>0 sufficiently small and N sufficiently large so that: (i) the feature/window approximations introduce at most CδCδ error in the measured coordinates (Steps 1–2), (i) the attention+linear block approximates the modal coupling mK_m within CδCδ uniformly on bounded sets (Step 3), and (i) the scatter+decoder approximates the output features within CδCδ uniformly (Step 4). By stability (continuity) of all stages, ‖Θ−‖≤‖Θ−m‖⏟≤Cδ+‖m−‖⏟εm<Cδ+ε/3. \|G_ -T \|\;≤\; \|G_ -T_m \|_≤ Cδ\;+\; \|T_m-T \|_ _m\;<\;Cδ+ /3. Choosing δ so that Cδ<2ε/3Cδ<2 /3 yields ‖Θ−‖<ε\|G_ -T\|< . Combining the steps completes the proof. ∎ Appendix C Implementation Details C.1 Baseline implementations Baseline models (Geo-FNO, LSM, Galerkin Transformer, GNOT, ONO, Transolver) are adapted from the Neural-Solver-Library (Wu et al., 2024) reference implementation at https://github.com/thuml/Neural-Solver-Library. Other models are adapted from their official repositories. Unless otherwise noted, we keep an identical training schedule across baselines: AdamW optimizer, initial learning rate 10−310^-3 with a StepLR scheduler (step size and decay factor as in the library’s default per dataset), up to 500 epochs with validation early stopping, the same data normalization/inverse-normalization protocol, and matched rollout/evaluation settings. C.2 GPO configurations The dataset-specific configurations of GPO are summarized in Table 2. We provide the source code of GPO in the Supplementary Material. Table 2: Model configurations of GPO. Benchmarks Model Configurations hidden_dim num_layers num_heads num_gaussians NS2D 128 8 8 32 NS3D 64 8 4 16 ERA5-temp 64 4 4 16 ERA5-wind u 64 4 4 16 Carra 64 4 4 16 Airfoil 64 4 4 16 Turbulent 128 8 8 32 PlanetSWE 64 8 4 16 Appendix D Model Analysis D.1 Ablation Study Table 3 reports the L2L_2 error under controlled variants (parameter counts are adjusted to be comparable). (i) Necessity of the Gaussian Field. Replacing the Gaussian Field with plain MLP encoder/decoder (w/o Gaussian Field) degrades accuracy markedly (7.44×10−27.44×10^-2), and removing the PG operator while keeping the Gaussian Field (w/o PG Operator) is even worse (8.57×10−28.57×10^-2). Notably, even a single Gaussian per site (num_gaussian=1) already improves to 6.28×10−26.28×10^-2, indicating that particleized Gaussian evaluation is a beneficial inductive bias beyond a black-box MLP. (i) Synergy of PG Operator and Gaussian Field. Combining the Gaussian basis with the PG Gaussian Attention yields the full GPO (baseline: 3.90×10−23.90×10^-2), demonstrating that the PG measurement→ coupling→ complements the local particle representation; each component alone is insufficient. (i) Effect of the number of Gaussians. Increasing num_gaussian consistently reduces error (from 4.21×10−24.21×10^-2 at G=16G=16 to 3.84×10−23.84×10^-2 at G=64G=64), but with diminishing returns; considering cost (Sec. D.2), we adopt G=16/32G=16/32 as a practical trade-off between efficiency and accuracy. Table 3: Ablation results comparing the L2L_2 error of different configurations. Model Configuration L2L_2 Error w/o PG Operator 8.57E-02 w/o Gaussian Field 7.44E-02 num_gaussian = 1 6.28E-02 num_gaussian = 16 4.21E-02 num_gaussian = 64 3.84E-02 GPO (Baseline) 3.90E-02 D.2 Computational Complexity Empirical measurements (Table 4, 64×64×364×64×3, batch 16) corroborate the analysis: GPO attains low memory footprint (2,313 MiB) and competitive time (44.66 s/epoch train; 1.67 s/epoch inference) with a modest parameter count (6.10 MB), outperforming attention baselines in training speed (Galerkin/Transolver/ONO/GNOT) and GPU memory, while remaining close to spectral baselines at inference. Although FNO is fastest on this small grid, GPO’s cost grows near–linearly with N and remains stable when moving to higher resolutions or 3D, where spatial attention becomes prohibitive and FFT memory/IO costs rise. By aggregating locally (N↔GN\! \!G) and coupling globally only in modal space (G×G×G), GPO delivers resolution–agnostic efficiency: linear scaling in N, controllable quadratic dependence on G, and favorable memory/time trade–offs across 2D/3D and irregular domains. Table 4: Computational efficiency comparison across models (measured with input size 64×64×3, batch size 16). Model Param Count Param (MB) GPU Mem (MiB) Train (s/epoch) Inference (s/epoch) FNO 640,305 4.84 949 28.27 0.5 LSM 19,187,457 73.23 2,875 48.42 1.73 Galerkin Transformer 1,096,321 4.18 4,301 65.29 2.96 GNOT 2,485,901 9.48 8,643 139.42 6.09 Transolver 3,069,889 11.71 4,917 97.03 4.10 ONO 1,596,673 6.09 6,163 94.80 4.27 GPO (Ours) 1,598,257 6.10 2,313 44.66 1.67 Appendix E Additional Visualizations Figure 7: Interpretable visualization of ERA5 on an in-distribution sample (above) and an out-of-distribution sample (below). Left to right: ground truth, reconstruction from the Gaussian basis and learned Gaussian particles overlaid (ellipses: center μ, axes ∝σ σ, color/size ∝w w).