Physics-Integrated Operator Learning via Gaussian Splatting Representations

# Physics-Integrated Operator Learning via Gaussian Splatting Representations
Source: [https://arxiv.org/html/2608.24049](https://arxiv.org/html/2608.24049)
Journal:ElsevierJunyi GuoJian\-Xun WangEmail:[jw2837@cornell\.edu](mailto:[email protected])Address:Sibley School of Mechanical and Aerospace Engineering, Cornell University, Ithaca, NY, USACorresponding author:Corresponding author\. Tel: \+1 540 3156512

###### Abstract

Neural operators provide efficient surrogates for spatiotemporal PDE systems, but purely data\-driven formulations often accumulate substantial errors during long\-horizon autoregressive prediction and may fail to exploit available governing\-equation structure\. Existing approaches incorporate physics primarily through residual\-based training objectives or PDE\-specific architectural constraints, which can introduce optimization difficulties or limit architectural generality\. In this work, we introduce a representation\-level approach to physics integration in which a feed\-forward Gaussian splatting \(FFGS\) representation serves as a continuous interface between discretized solution fields and governing operators\. The FFGS representation reconstructs the state as a continuous Gaussian field with closed\-form spatial derivatives, allowing available physical PDE operators to be integrated directly within the learned evolution map without introducing a physics\-residual loss\. We evaluate the framework across two\- and three\-dimensional PDE systems, including advection, diffusion, nonlinear self\-advection, and reaction dynamics\. Over long\-horizon autoregressive rollouts, the proposed framework reduces relativeℓ2\\ell\_\{2\}error by1\.5×1\.5\\times–2\.2×2\.2\\timescompared with the strongest purely data\-driven baseline across the benchmark suite, while consistently improving spectral fidelity\. The framework also remains effective when the governing equations are partially known, demonstrating robustness to incomplete physics\. These results demonstrate that continuous field representations can provide a practical interface for incorporating known physical structure into generic neural\-operator surrogates\.

## 1Introduction

Spatiotemporal dynamics governed by partial differential equations \(PDEs\) are ubiquitous in science and engineering, from convective and diffusive transport to nonlinear wave steepening and moving reaction interfaces\. Traditional numerical methods discretize the governing operator and advance the state by time integration; they resolve these dynamics accurately, but at a cost that becomes prohibitive in many\-query settings such as design optimization, inverse problems, data assimilation, and uncertainty quantification, where the forward model must be evaluated repeatedly across different initial conditions and physical parameters\. This has motivated a rapidly growing effort to build data\-driven surrogates that learn a fast approximation of the solution operator from precomputed trajectories and, once trained, stand in for the numerical solver at a small fraction of its cost\.

Neural operators \(NOs\)\[[1](https://arxiv.org/html/2608.24049#bib.bib1),[2](https://arxiv.org/html/2608.24049#bib.bib2)\], which use neural networks to parameterize mappings between infinite\-dimensional function spaces, are currently the leading class of such surrogates\. They have demonstrated strong predictive performance for a broad range of spatiotemporal systems, including weather dynamics\[[3](https://arxiv.org/html/2608.24049#bib.bib3),[4](https://arxiv.org/html/2608.24049#bib.bib4)\]and turbulent flows\[[5](https://arxiv.org/html/2608.24049#bib.bib5),[6](https://arxiv.org/html/2608.24049#bib.bib6)\]\. However, standard NOs are primarily trained from trajectory data and infer the underlying dynamics implicitly from finite observations\. This limitation becomes particularly evident in autoregressive prediction, where local approximation errors accumulate through repeated operator evaluations and progressively degrade long\-horizon accuracy\[[7](https://arxiv.org/html/2608.24049#bib.bib7)\]\. The degradation is further amplified when the deployed dynamics deviate from the training distribution\[[8](https://arxiv.org/html/2608.24049#bib.bib8),[9](https://arxiv.org/html/2608.24049#bib.bib9)\]\. These failure modes motivate the integration of governing\-equation knowledge into operator learning\.

A direct strategy is to augment supervised operator learning with auxiliary loss terms derived from the governing equations\[[10](https://arxiv.org/html/2608.24049#bib.bib10),[11](https://arxiv.org/html/2608.24049#bib.bib11)\]\. Such physics\-informed formulations inherit several optimization difficulties of PINNs\[[12](https://arxiv.org/html/2608.24049#bib.bib12)\]: the data and physics objectives exhibit disparate scales and gradient dynamics, making their relative weighting delicate and problem\-dependent\[[13](https://arxiv.org/html/2608.24049#bib.bib13),[14](https://arxiv.org/html/2608.24049#bib.bib14)\]\. Moreover, enforcing the governing equations through a residual requires repeated evaluation of spatial and temporal derivatives during training\. For coordinate\-based representations, these derivatives are commonly obtained through automatic differentiation, whose computational and memory costs increase with the number of collocation points and derivative order, thereby limiting the scalability of residual\-based physics integration\.

An alternative direction incorporates physical knowledge directly into the construction of the operator\. For example, conservation\-preserving neural operators enforce prescribed invariants through differential constructions\[[15](https://arxiv.org/html/2608.24049#bib.bib15)\], while invariant and equivariant operators encode symmetry properties of the governing system\[[16](https://arxiv.org/html/2608.24049#bib.bib16),[17](https://arxiv.org/html/2608.24049#bib.bib17)\]\. Other approaches introduce local differential interactions or physics\-motivated feature transformations into neural\-operator layers\[[18](https://arxiv.org/html/2608.24049#bib.bib18),[19](https://arxiv.org/html/2608.24049#bib.bib19)\]\. These methods demonstrate that physical knowledge can improve learned dynamics when incorporated beyond the training objective\. However, the corresponding constructions are often tailored to specific PDE systems, and incorporating a different governing equation may require redesigning the corresponding physics\-aware components\.

This PDE\-specific design motivates introducing governing knowledge through the field representation rather than through the neural operator architecture\. Such a representation must both reconstruct the physical field from its discretized state and provide direct access to the spatial differential quantities required by the governing operator\. Explicit continuous basis representations are well suited to this role because both the reconstructed field and its spatial derivatives can be evaluated analytically from the basis functions\. The available governing terms can then be evaluated directly on this continuous representation while retaining a generic neural operator architecture\.

The use of explicit continuous basis representations is well established in classical radial basis function \(RBF\) collocation methods, where the solution is represented by smooth basis functions and differential operators are applied analytically to the basis to enforce the governing equations at collocation points\[[20](https://arxiv.org/html/2608.24049#bib.bib20),[21](https://arxiv.org/html/2608.24049#bib.bib21)\]\. More recently, Gaussian primitives obtained by 3D Gaussian Splatting\[[22](https://arxiv.org/html/2608.24049#bib.bib22)\]have been explored for physical\-field representation and modeling\. Existing studies span Gaussian representations for compact field encoding\[[23](https://arxiv.org/html/2608.24049#bib.bib23)\], Gaussian\-basis operator learning\[[24](https://arxiv.org/html/2608.24049#bib.bib24)\], and PDE\-constrained Gaussian optimization\[[25](https://arxiv.org/html/2608.24049#bib.bib25)\]\.

The original Gaussian\-splatting paradigm, however, relies on per\-instance optimization of the primitive parameters\. Applying such a procedure to every evolving PDE state would introduce substantial overhead in autoregressive time stepping\. Feed\-forward Gaussian\-splatting \(FFGS\) methods developed in computer vision provide an efficient alternative by amortizing Gaussian construction into learned models that predict primitive parameters in a single forward pass\[[26](https://arxiv.org/html/2608.24049#bib.bib26),[27](https://arxiv.org/html/2608.24049#bib.bib27),[28](https://arxiv.org/html/2608.24049#bib.bib28)\]\. Related amortized Gaussian representations have also begun to appear in scientific applications, for example for reconstructing physical fields from sparse observations\[[29](https://arxiv.org/html/2608.24049#bib.bib29)\]\. However, the use of an amortized Gaussian representation as a continuous interface for evaluating governing operators within neural\-operator time evolution remains largely unexplored\.

To address this gap, we develop a physics\-integrated neural\-operator framework that employs FFGS as a continuous representation of physical fields\. The analytic Gaussian representation provides direct access to the spatial derivatives required to evaluate available governing terms, allowing their contribution to be incorporated into the learned evolution map while retaining a generic neural\-operator backbone\. We further extend the framework to a partial\-physics setting, where the governing form is known but its coefficients are identified from trajectory data\. We evaluate the framework across two\- and three\-dimensional systems involving advection, diffusion, nonlinear self\-advection, and reaction dynamics\. Across these settings, representation\-level physics integration consistently improves long\-horizon rollout accuracy and spectral fidelity over purely data\-driven baselines, with substantial gains retained when only partial governing information is available\. Component ablations further show that the embedded\-physics and neural\-operator components play complementary roles in the resulting dynamics\. Overall, these results establish FFGS as a practical continuous interface for incorporating fully or partially known governing structure into generic neural\-operator surrogates\.

The remainder of this paper is structured as follows\. In Sec\.[2](https://arxiv.org/html/2608.24049#S2), we introduce the physics\-integrated neural\-operator framework and formulate the FFGS representation, embedded\-physics component, and rollout\-training procedure\. In Sec\.[3](https://arxiv.org/html/2608.24049#S3), we describe the benchmark problems, training protocol, baseline models, and evaluation metrics\. In Sec\.[4](https://arxiv.org/html/2608.24049#S4), we report the numerical results and evaluate long\-horizon rollout accuracy, spectral fidelity, partial\-physics performance, and inference cost\. Sec\.[5](https://arxiv.org/html/2608.24049#S5)analyzes the contributions of the individual components and discusses the limitations and future directions of the framework, followed by conclusions in Sec\.[6](https://arxiv.org/html/2608.24049#S6)\.

## 2Methodology

![Refer to caption](https://arxiv.org/html/2608.24049v1/FFGS.png)\(a\)Stage 1:FFGS representation and local Gaussian rendering\.
![Refer to caption](https://arxiv.org/html/2608.24049v1/big_architecture.png)\(b\)Stage 2:Overall physics\-integrated surrogate\.

Figure 1:Method overview\. \(a\) The FFGS module encodes a normalized grid field into anisotropic Gaussian parameters and reconstructs a continuous field by local rendering\. \(b\) The composite one\-step surrogate combines an FFGS\-based embedded\-physics component with a neural\-operator component\. The embedded\-physics component supplies closed\-form spatial derivatives, applies the low\-pass filter, and evaluates the available governing terms inside the update map; gradients are stopped through this component during rollout training\.### 2\.1Overview

We consider a general class of nonlinear, coupled PDE systems for aCC\-component field in a parametric setting,

∂𝐮∂t\\displaystyle\\frac\{\\partial\\mathbf\{u\}\}\{\\partial t\}=ℱ⁡\[𝐮,∇𝐱𝐮,∇𝐱2𝐮,…;𝝀\]\\displaystyle=\\mathcal\{F\}\\\!\\left\[\\mathbf\{u\},\\nabla\_\{\\mathbf\{x\}\}\\mathbf\{u\},\\nabla\_\{\\mathbf\{x\}\}^\{2\}\\mathbf\{u\},\\ldots;\\boldsymbol\{\\lambda\}\\right\]\(1\)ℐ⁡\[𝐮,𝐮𝟎\]\\displaystyle\\mathcal\{I\}\[\\mathbf\{u\},\\mathbf\{u\_\{0\}\}\]=𝟎,t=0,𝐱∈Ω,\\displaystyle=\\mathbf\{0\},\\qquad t=0,\\ \\mathbf\{x\}\\in\\Omega,ℬ⁡\[𝐮,∇𝐱𝐮,…\]\\displaystyle\\mathcal\{B\}\[\\mathbf\{u\},\\nabla\_\{\\mathbf\{x\}\}\\mathbf\{u\},\\ldots\]=𝟎,𝐱∈∂Ω\.\\displaystyle=\\mathbf\{0\},\\qquad\\mathbf\{x\}\\in\\partial\\Omega\.where𝐮⁡\(𝐱,t\)∈ℝC\\mathbf\{u\}\(\\mathbf\{x\},t\)\\in\\mathbb\{R\}^\{C\}is the solution field over the spatial domainΩ⊂ℝd\\Omega\\subset\\mathbb\{R\}^\{d\}and temporal horizont∈\[0,T\]t\\in\[0,T\];ℱ⁡\[⋅;𝝀\]\\mathcal\{F\}\[\\cdot;\\boldsymbol\{\\lambda\}\]is a nonlinear operator acting on the field and its spatial derivatives, with𝝀\\boldsymbol\{\\lambda\}denoting the physical parameters; andℐ\\mathcal\{I\}andℬ\\mathcal\{B\}encode the initial and boundary conditions, respectively\. When the spatial argument is omitted,𝐮⁡\(t\)\\mathbf\{u\}\(t\)denotes the full field𝐮⁡\(⋅,t\)\\mathbf\{u\}\(\\cdot,t\)\.

The governing PDE induces a finite\-time solution operator𝒢Δ​t\\mathcal\{G\}\_\{\\Delta t\}through the temporal evolution generated byℱ\\mathcal\{F\}\. In an explicit time\-stepping manner, this gives

𝐮⁡\(t\+Δ​t\)\\displaystyle\\mathbf\{u\}\(t\+\\Delta t\)=𝒢Δ​t​\(𝐮⁡\(t\),𝝀\)\\displaystyle=\\mathcal\{G\}\_\{\\Delta t\}\\\!\\left\(\\mathbf\{u\}\(t\);\\boldsymbol\{\\lambda\}\\right\)\(2\)=𝐮⁡\(t\)\+∫tt\+Δ​tℱ⁡\[𝐮⁡\(τ\),∇𝐱𝐮​\(τ\),∇𝐱2𝐮​\(τ\),…;𝝀\]​dτ\.\\displaystyle=\\mathbf\{u\}\(t\)\+\\int\_\{t\}^\{t\+\\Delta t\}\\mathcal\{F\}\\\!\\left\[\\mathbf\{u\}\(\\tau\),\\nabla\_\{\\mathbf\{x\}\}\\mathbf\{u\}\(\\tau\),\\nabla\_\{\\mathbf\{x\}\}^\{2\}\\mathbf\{u\}\(\\tau\),\\ldots;\\boldsymbol\{\\lambda\}\\right\]\\,d\\tau\.Numerical evaluation of this evolution generally requires a time stepΔ​t\\Delta tconstrained by the stability and accuracy requirements of the underlying integration scheme\. Neural operators instead approximate the finite\-time solution operator𝒢Δ​t\\mathcal\{G\}\_\{\\Delta t\}with a parameterized operator𝒩θ\\mathcal\{N\}\_\{\\theta\}, allowing prediction over the larger temporal interval\.

To integrate the available physics through representation, we construct a physics\-based approximation of the finite\-time solution operator𝒢Δ​t\\mathcal\{G\}\_\{\\Delta t\}\. Specifically, an operatorFFGS\\operatorname\{FFGS\}reconstructs𝐮⁡\(t\)\\mathbf\{u\}\(t\)as a continuous field whose closed\-form spatial derivatives provide the quantities required to evaluate the governing operatorℱ\\mathcal\{F\}\. The resulting governing terms are then integrated overΔ​t\\Delta t, defining the embedded\-physics operatorΦΔ​tFFGS​\(𝐮⁡\(t\),𝝀\)\\Phi\_\{\\Delta t\}^\{\\mathrm\{FFGS\}\}\(\\mathbf\{u\}\(t\);\\boldsymbol\{\\lambda\}\)\. We combine this embedded\-physics operator with a neural operator to approximate𝒢Δ​t\\mathcal\{G\}\_\{\\Delta t\},

𝐮⁡\(t\+Δ​t\)=𝒢Δ​t​\(𝐮⁡\(t\),𝝀\)≈𝒩θ​\(𝐮​\(t\),𝝀\)⏟neural\-operator component\+ΦΔ​tFFGS​\(𝐮⁡\(t\),𝝀\)⏟embedded\-physics component\.\\mathbf\{u\}\(t\+\\Delta t\)=\\mathcal\{G\}\_\{\\Delta t\}\\\!\\left\(\\mathbf\{u\}\(t\);\\boldsymbol\{\\lambda\}\\right\)\\approx\\underbrace\{\\mathcal\{N\}\_\{\\theta\}\\\!\\left\(\\mathbf\{u\}\(t\);\\boldsymbol\{\\lambda\}\\right\)\}\_\{\\text\{neural\-operator component\}\}\+\\underbrace\{\\Phi\_\{\\Delta t\}^\{\\mathrm\{FFGS\}\}\\\!\\left\(\\mathbf\{u\}\(t\);\\boldsymbol\{\\lambda\}\\right\)\}\_\{\\text\{embedded\-physics component\}\}\.\(3\)The following sections constructΦΔ​tFFGS\\Phi\_\{\\Delta t\}^\{\\mathrm\{FFGS\}\}from the Gaussian representation and its closed\-form spatial derivatives, and then describe the rollout training of the composite operator in Eq\. \([3](https://arxiv.org/html/2608.24049#S2.E3)\)\.

### 2\.2FFGS differentiable field representation

The embedded\-physics component in Eq\. \([3](https://arxiv.org/html/2608.24049#S2.E3)\) requires evaluating the spatial differential quantities appearing in the governing operatorℱ\\mathcal\{F\}\. FFGS provides these quantities by representing the current state𝐮⁡\(t\)\\mathbf\{u\}\(t\)with a continuous Gaussian expansion whose spatial derivatives are available in closed form\. State\-dependent algebraic quantities, by contrast, are evaluated directly from the current state rather than from its FFGS reconstruction\. To construct the continuous representation, FFGS predicts a set of Gaussian primitives from𝐮⁡\(t\)\\mathbf\{u\}\(t\)and renders the corresponding basis expansion\.

Given the current state𝐮⁡\(t\)\\mathbf\{u\}\(t\), a convolutional encoderℰφ\\mathcal\{E\}\_\{\\varphi\}predictsGGanisotropic Gaussian primitives,

\{\(𝝁i,Ri,𝝈i,𝜶i\)\}i=1G=ℰφ​\(𝐮⁡\(t\)\)\.\\left\\\{\\left\(\\boldsymbol\{\\mu\}\_\{i\},R\_\{i\},\\boldsymbol\{\\sigma\}\_\{i\},\\boldsymbol\{\\alpha\}\_\{i\}\\right\)\\right\\\}\_\{i=1\}^\{G\}=\\mathcal\{E\}\_\{\\varphi\}\\\!\\left\(\\mathbf\{u\}\(t\)\\right\)\.\(4\)where𝝁i∈Ω\\boldsymbol\{\\mu\}\_\{i\}\\in\\Omegais the center of theii\-th Gaussian,RiR\_\{i\}specifies its orientation,𝝈i\\boldsymbol\{\\sigma\}\_\{i\}contains its principal\-axis scales, and𝜶i∈ℝC\\boldsymbol\{\\alpha\}\_\{i\}\\in\\mathbb\{R\}^\{C\}contains its channel\-wise amplitudes\. The orientation and scales define the symmetric positive\-definite covariance matrix

Σi=Ri​diag⁡\(𝝈i2\)​Ri⊤,\\Sigma\_\{i\}=R\_\{i\}\\operatorname\{diag\}\(\\boldsymbol\{\\sigma\}\_\{i\}^\{2\}\)R\_\{i\}^\{\\top\},\(5\)which controls the anisotropic spatial extent of each primitive and allows the representation to adapt to directional variations in the field\. Further parameterization details are provided in[A](https://arxiv.org/html/2608.24049#A1)\.

The feed\-forward encoder replaces the per\-state optimization otherwise required to fit the Gaussian primitives\. This amortization is important for autoregressive evolution, where a continuous representation must be reconstructed repeatedly as new states are encountered\. A convolutional architecture is used to predict primitive parameters on a spatial anchor lattice with weights shared across locations; architectural details are deferred to[A](https://arxiv.org/html/2608.24049#A1)\.

The predicted primitives define the Gaussian kernels and the corresponding continuous field reconstruction,

gi​\(𝐱\)=exp⁡\[−12​\(𝐱−𝝁i\)⊤​Σi−1​\(𝐱−𝝁i\)\],𝐮^​\(𝐱q,t\)=∑i∈𝒮⁡\(q\)𝜶i​gi​\(𝐱q\)\.g\_\{i\}\(\\mathbf\{x\}\)=\\exp\\\!\\left\[\-\\frac\{1\}\{2\}\(\\mathbf\{x\}\-\\boldsymbol\{\\mu\}\_\{i\}\)^\{\\top\}\\Sigma\_\{i\}^\{\-1\}\(\\mathbf\{x\}\-\\boldsymbol\{\\mu\}\_\{i\}\)\\right\],\\qquad\\hat\{\\mathbf\{u\}\}\(\\mathbf\{x\}\_\{q\},t\)=\\sum\_\{i\\in\\mathcal\{S\}\(q\)\}\\boldsymbol\{\\alpha\}\_\{i\}g\_\{i\}\(\\mathbf\{x\}\_\{q\}\)\.\(6\)where𝒮⁡\(q\)\\mathcal\{S\}\(q\)denotes the fixed local set of Gaussian primitives used to render the query point𝐱q\\mathbf\{x\}\_\{q\}\. Restricting the rendering to this neighborhood avoids evaluating allGGprimitives at every query point while exploiting the spatial decay of the Gaussian kernels\. The rendering window, periodic boundary treatment, and approximation error associated with the local truncation are detailed in[A](https://arxiv.org/html/2608.24049#A1)\.

The analytic structure of the Gaussian basis in Eq\. \([6](https://arxiv.org/html/2608.24049#S2.E6)\) allows its spatial derivatives to be evaluated in closed form\. For each primitive,

∇𝐱gi​\(𝐱\)\\displaystyle\\nabla\_\{\\mathbf\{x\}\}g\_\{i\}\(\\mathbf\{x\}\)=−gi​\(𝐱\)​Σi−1​\(𝐱−𝝁i\),\\displaystyle=\-g\_\{i\}\(\\mathbf\{x\}\)\\Sigma\_\{i\}^\{\-1\}\(\\mathbf\{x\}\-\\boldsymbol\{\\mu\}\_\{i\}\),\(7\)∇𝐱2​gi​\(𝐱\)\\displaystyle\\nabla\_\{\\mathbf\{x\}\}^\{2\}g\_\{i\}\(\\mathbf\{x\}\)=gi​\(𝐱\)​\[\(𝐱−𝝁i\)⊤​Σi−2​\(𝐱−𝝁i\)−tr⁡\(Σi−1\)\],\\displaystyle=g\_\{i\}\(\\mathbf\{x\}\)\\left\[\(\\mathbf\{x\}\-\\boldsymbol\{\\mu\}\_\{i\}\)^\{\\top\}\\Sigma\_\{i\}^\{\-2\}\(\\mathbf\{x\}\-\\boldsymbol\{\\mu\}\_\{i\}\)\-\\operatorname\{tr\}\(\\Sigma\_\{i\}^\{\-1\}\)\\right\],whereΣi−2=Σi−1​Σi−1\\Sigma\_\{i\}^\{\-2\}=\\Sigma\_\{i\}^\{\-1\}\\Sigma\_\{i\}^\{\-1\}\. By linearity, the corresponding derivatives of the reconstructed field are

∇𝐱u^c​\(𝐱q,t\)\\displaystyle\\nabla\_\{\\mathbf\{x\}\}\\hat\{u\}\_\{c\}\(\\mathbf\{x\}\_\{q\},t\)=∑i∈𝒮⁡\(q\)αi,c​∇𝐱gi​\(𝐱q\),\\displaystyle=\\sum\_\{i\\in\\mathcal\{S\}\(q\)\}\\alpha\_\{i,c\}\\nabla\_\{\\mathbf\{x\}\}g\_\{i\}\(\\mathbf\{x\}\_\{q\}\),\(8\)∇𝐱2​u^c​\(𝐱q,t\)\\displaystyle\\nabla\_\{\\mathbf\{x\}\}^\{2\}\\hat\{u\}\_\{c\}\(\\mathbf\{x\}\_\{q\},t\)=∑i∈𝒮⁡\(q\)αi,c∇𝐱2gi\(𝐱q\),c=1,…,C\.\\displaystyle=\\sum\_\{i\\in\\mathcal\{S\}\(q\)\}\\alpha\_\{i,c\}\\nabla\_\{\\mathbf\{x\}\}^\{2\}g\_\{i\}\(\\mathbf\{x\}\_\{q\}\),\\qquad c=1,\\ldots,C\.Thus, the closed\-form spatial derivatives required byℱ\\mathcal\{F\}are evaluated directly from the reconstructed continuous field, without finite\-difference stencils or automatic differentiation with respect to the spatial coordinates\.

Spatial differentiation can amplify high\-wavenumber reconstruction errors\. We therefore apply a low\-pass spectral filterΠκ\\Pi\_\{\\kappa\}to the differential quantities obtained from the FFGS representation before they enter the governing operator,

∇𝐱𝐮~​\(t\)=Πκ​∇𝐱𝐮^​\(t\),∇𝐱2​𝐮~​\(t\)=Πκ​∇𝐱2𝐮^​\(t\)\.\\widetilde\{\\nabla\_\{\\mathbf\{x\}\}\\mathbf\{u\}\}\(t\)=\\Pi\_\{\\kappa\}\\nabla\_\{\\mathbf\{x\}\}\\hat\{\\mathbf\{u\}\}\(t\),\\qquad\\widetilde\{\\nabla\_\{\\mathbf\{x\}\}^\{2\}\\mathbf\{u\}\}\(t\)=\\Pi\_\{\\kappa\}\\nabla\_\{\\mathbf\{x\}\}^\{2\}\\hat\{\\mathbf\{u\}\}\(t\)\.\(9\)The filter is implemented in Fourier space using either a sharp spectral cutoff or a smooth high\-wavenumber roll\-off, with benchmark\-specific settings detailed in[B](https://arxiv.org/html/2608.24049#A2)\. It is applied to the closed\-form spatial derivatives obtained from the FFGS reconstruction, suppressing high\-frequency reconstruction errors before the filtered derivatives are passed toℱ\\mathcal\{F\}\. The current grid\-resolved state𝐮⁡\(t\)\\mathbf\{u\}\(t\), in contrast, is retained directly for state\-dependent algebraic terms in the governing operator\. Thus, FFGS serves as a continuous differential interface, providing the required spatial derivatives without replacing the current state in purely algebraic terms\.

The encoder is trained independently by reconstructing instantaneous field snapshots on the observation grid\. The reconstruction objective is

ℒFFGS​\(φ\)=1\|ℳ\|​∑𝐮⁡\(t\)∈ℳ1Nd​C​∑j=1Nd‖𝐮^​\(𝐱j,t\)−𝐮⁡\(𝐱j,t\)‖22,\\mathcal\{L\}\_\{\\mathrm\{FFGS\}\}\(\\varphi\)=\\frac\{1\}\{\|\\mathcal\{M\}\|\}\\sum\_\{\\mathbf\{u\}\(t\)\\in\\mathcal\{M\}\}\\frac\{1\}\{N^\{d\}C\}\\sum\_\{j=1\}^\{N^\{d\}\}\\left\\\|\\hat\{\\mathbf\{u\}\}\(\\mathbf\{x\}\_\{j\},t\)\-\\mathbf\{u\}\(\\mathbf\{x\}\_\{j\},t\)\\right\\\|\_\{2\}^\{2\},\(10\)whereℳ\\mathcal\{M\}denotes a mini\-batch of field snapshots\. After this reconstruction stage, the parametersφ\\varphiare frozen\.

During the composite update, the frozen FFGS map constructs a continuous representation of the current state, from which the required closed\-form spatial derivatives are evaluated through Eqs\. \([7](https://arxiv.org/html/2608.24049#S2.E7)\)–\([8](https://arxiv.org/html/2608.24049#S2.E8)\)\. The resulting differential quantities are then spectrally filtered according to Eq\. \([9](https://arxiv.org/html/2608.24049#S2.E9)\)\. When assembling the available terms of the governing operatorℱ\\mathcal\{F\}, state\-dependent algebraic quantities are evaluated directly from the current state𝐮⁡\(t\)\\mathbf\{u\}\(t\), while spatial differential quantities are supplied by the filtered FFGS derivatives\. The resulting physics\-based right\-hand side is subsequently integrated over the intervalΔ​t\\Delta tto obtain the embedded\-physics operatorΦΔ​tFFGS​\(𝐮⁡\(t\),𝝀\)\\Phi\_\{\\Delta t\}^\{\\mathrm\{FFGS\}\}\\\!\\left\(\\mathbf\{u\}\(t\);\\boldsymbol\{\\lambda\}\\right\)appearing in Eq\. \([3](https://arxiv.org/html/2608.24049#S2.E3)\)\. The specific time\-integration scheme used in the implementation is detailed in[B](https://arxiv.org/html/2608.24049#A2)\.

### 2\.3Neural\-operator component and rollout training

The neural\-operator component𝒩θ\\mathcal\{N\}\_\{\\theta\}in Eq\. \([3](https://arxiv.org/html/2608.24049#S2.E3)\) is implemented using a standard operator\-learning backbone, with the specific architecture deferred to[C](https://arxiv.org/html/2608.24049#A3)\. When the physical parameters vary across trajectories,𝝀\\boldsymbol\{\\lambda\}is additionally provided to𝒩θ\\mathcal\{N\}\_\{\\theta\}as conditioning information; this dependence is omitted in fixed\-parameter settings\.

With the FFGS parametersφ\\varphifrozen, the physics\-integrated surrogate is trained through multi\-step autoregressive rollout\. Starting from a ground\-truth state𝐮⁡\(t\)\\mathbf\{u\}\(t\), we set𝐮​\(t\)pred=𝐮⁡\(t\)\\mathbf\{u\}\(t\)^\{\\mathrm\{pred\}\}=\\mathbf\{u\}\(t\)and recursively apply the one\-step map in Eq\. \([3](https://arxiv.org/html/2608.24049#S2.E3)\) forHHsteps,

𝐮pred​\(t\+k​Δ​t\)\\displaystyle\\mathbf\{u\}^\{\\mathrm\{pred\}\}\(t\+k\\Delta t\)=𝒢^Δ​t,θ​\(𝐮pred​\(t\+\(k−1\)​Δ​t\),𝝀\)\\displaystyle=\\widehat\{\\mathcal\{G\}\}\_\{\\Delta t,\\theta\}\\\!\\left\(\\mathbf\{u\}^\{\\mathrm\{pred\}\}\(t\+\(k\-1\)\\Delta t\);\\boldsymbol\{\\lambda\}\\right\)\(11\)=sg\[ΦΔ​tFFGS\(𝐮pred\(t\+\(k−1\)Δt\);𝝀\)\]\+𝒩θ\(𝐮pred\(t\+\(k−1\)Δt\);𝝀\),k=1,…,H\.\\displaystyle=\\operatorname\{sg\}\\\!\\left\[\\Phi\_\{\\Delta t\}^\{\\mathrm\{FFGS\}\}\\\!\\left\(\\mathbf\{u\}^\{\\mathrm\{pred\}\}\(t\+\(k\-1\)\\Delta t\);\\boldsymbol\{\\lambda\}\\right\)\\right\]\+\\mathcal\{N\}\_\{\\theta\}\\\!\\left\(\\mathbf\{u\}^\{\\mathrm\{pred\}\}\(t\+\(k\-1\)\\Delta t\);\\boldsymbol\{\\lambda\}\\right\),\\qquad k=1,\\ldots,H\.wheresg⁡\[⋅\]\\operatorname\{sg\}\[\\cdot\]denotes the stop\-gradient operator\. Training on recursively predicted states exposes the neural operator to the states generated by the composite dynamics during autoregressive evolution, rather than restricting optimization to one\-step transitions from ground\-truth inputs\.

The neural\-operator parametersθ\\thetaare optimized using the mean squared error accumulated over the rollout horizon,

ℒrollout​\(θ\)=1H​∑k=1H1Nd​C​∑j=1Nd‖𝐮pred​\(𝐱j,t\+k​Δ​t\)−𝐮⁡\(𝐱j,t\+k​Δ​t\)‖22\.\\mathcal\{L\}\_\{\\mathrm\{rollout\}\}\(\\theta\)=\\frac\{1\}\{H\}\\sum\_\{k=1\}^\{H\}\\frac\{1\}\{N^\{d\}C\}\\sum\_\{j=1\}^\{N^\{d\}\}\\left\\\|\\mathbf\{u\}^\{\\mathrm\{pred\}\}\(\\mathbf\{x\}\_\{j\},t\+k\\Delta t\)\-\\mathbf\{u\}\(\\mathbf\{x\}\_\{j\},t\+k\\Delta t\)\\right\\\|\_\{2\}^\{2\}\.\(12\)The governing dynamics therefore enter the training procedure through the embedded\-physics componentΦΔ​tFFGS\\Phi\_\{\\Delta t\}^\{\\mathrm\{FFGS\}\}in the rollout itself, rather than through an auxiliary physics\-residual term in the loss\.

The stop\-gradient in Eq\. \([11](https://arxiv.org/html/2608.24049#S2.E11)\) excludes the embedded\-physics component from the backward computational graph, avoiding backpropagation through the repeated FFGS evaluations and numerical integrations\. Gradients are propagated through the neural\-operator paths of the unrolled computation to optimizeθ\\theta, while the embedded\-physics component is reevaluated at every forward step\. At inference, the same composite map𝒢^Δ​t,θ\\widehat\{\\mathcal\{G\}\}\_\{\\Delta t,\\theta\}is applied autoregressively over the target horizon\.

### 2\.4Extension to partially known physics

The formulation above assumes that the physical parameters𝝀\\boldsymbol\{\\lambda\}entering the embedded\-physics component are known\. We further consider a partially known setting in which the governing form is prescribed but some of its coefficients are unavailable\. In this work, these unknown coefficients are identified from trajectory data before rollout training and subsequently used in the same physics\-integrated surrogate\.

We consider governing operators whose dependence on the unknown coefficients is linear\. Specifically, the coefficient\-dependent part ofℱ\\mathcal\{F\}can be written as

ℱ⁡\[𝐮,∇𝐱𝐮,∇𝐱2𝐮,…;𝝀\]=∑m=1Pλm​ℱm​\[𝐮,∇𝐱𝐮,∇𝐱2𝐮,…\],\\mathcal\{F\}\\\!\\left\[\\mathbf\{u\},\\nabla\_\{\\mathbf\{x\}\}\\mathbf\{u\},\\nabla\_\{\\mathbf\{x\}\}^\{2\}\\mathbf\{u\},\\ldots;\\boldsymbol\{\\lambda\}\\right\]=\\sum\_\{m=1\}^\{P\}\\lambda\_\{m\}\\mathcal\{F\}\_\{m\}\\\!\\left\[\\mathbf\{u\},\\nabla\_\{\\mathbf\{x\}\}\\mathbf\{u\},\\nabla\_\{\\mathbf\{x\}\}^\{2\}\\mathbf\{u\},\\ldots\\right\],\(13\)whereℱm\\mathcal\{F\}\_\{m\}denotes the prescribed form of themm\-th governing term andλm\\lambda\_\{m\}is its unknown coefficient\. The frozen FFGS representation provides the closed\-form spatial derivatives required to evaluate eachℱm\\mathcal\{F\}\_\{m\}, while the temporal derivative is estimated from neighboring trajectory frames\. Evaluating Eq\. \([13](https://arxiv.org/html/2608.24049#S2.E13)\) over the observed spatial grid therefore yields a linear regression system for𝝀\\boldsymbol\{\\lambda\}\.

Collecting the evaluated governing terms into a regression matrix𝐀\\mathbf\{A\}and the estimated temporal derivative into𝐛\\mathbf\{b\}, the coefficient estimate is obtained from

𝐀​𝝀≈𝐛,𝝀^=arg​min𝝀⁡‖𝐀​𝝀−𝐛‖22\.\\mathbf\{A\}\\boldsymbol\{\\lambda\}\\approx\\mathbf\{b\},\\qquad\\hat\{\\boldsymbol\{\\lambda\}\}=\\operatorname\*\{arg\\,min\}\_\{\\boldsymbol\{\\lambda\}\}\\left\\\|\\mathbf\{A\}\\boldsymbol\{\\lambda\}\-\\mathbf\{b\}\\right\\\|\_\{2\}^\{2\}\.\(14\)The construction of𝐀\\mathbf\{A\}follows directly from the prescribed governing terms and the FFGS\-evaluated spatial derivatives\. The temporal\-difference scheme, regression normalization, and trajectory\-wise calibration procedure used to construct𝐛\\mathbf\{b\}and solve Eq\. \([14](https://arxiv.org/html/2608.24049#S2.E14)\) are detailed in[D](https://arxiv.org/html/2608.24049#A4)\.

After identification,𝝀^\\hat\{\\boldsymbol\{\\lambda\}\}is fixed for the corresponding trajectory and used in place of𝝀\\boldsymbol\{\\lambda\}in the composite evolution\. In particular, the embedded\-physics component is evaluated asΦΔ​tFFGS​\(𝐮⁡\(t\),𝝀^\)\\Phi\_\{\\Delta t\}^\{\\mathrm\{FFGS\}\}\(\\mathbf\{u\}\(t\);\\hat\{\\boldsymbol\{\\lambda\}\}\), and the same coefficient estimate is supplied to𝒩θ\\mathcal\{N\}\_\{\\theta\}when coefficient conditioning is used\. No further modification of the FFGS representation or the composite update is required\. The partially known setting therefore differs from the prescribed\-coefficient setting only in how the physical parameters entering Eq\. \([3](https://arxiv.org/html/2608.24049#S2.E3)\) are obtained; subsequent rollout training follows Eq\. \([11](https://arxiv.org/html/2608.24049#S2.E11)\) with the identified coefficients\.

## 3Experimental setup

This section describes the benchmark suite and data generation, the model instantiation and training protocol, the baselines, and the evaluation metrics\. Architecture\-level implementation details are provided in the appendices\.

### 3\.1Benchmark suite and data generation

We evaluate the proposed surrogate in two regimes\. In the fixed\-coefficient full\-physics regime, both the governing form and coefficients are prescribed\. We use five benchmarks in this regime to cover pure transport, transport–diffusion, nonlinear self\-advection, reaction\-driven moving interfaces, and three\-dimensional transport–diffusion \(Table[1](https://arxiv.org/html/2608.24049#S3.T1)\)\. In the coefficient\-unknown regime, the governing form is provided but the coefficients are withheld and must be identified from trajectory data before rollout training\.

All quantities are nondimensional, and all fixed\-coefficient problems are posed on the unit periodic domain\[0,1\)d\[0,1\)^\{d\}\. Initial conditions are sampled as random truncated Fourier series with modes restricted to\|ki\|≤kmax=5\|k\_\{i\}\|\\leq k\_\{\\max\}=5along each spatial axis\. Each benchmark contains 256 training trajectories and 30 held\-out test trajectories\. Models are trained on the first 50 steps and evaluated by 200\-step autoregressive rollout\. The velocity is prescribed as𝐯=v​𝟏\\mathbf\{v\}=v\\mathbf\{1\}, where𝟏\\mathbf\{1\}is the vector of ones\. The symbolscc,DD,ν\\nu, andrrdenote the convection scaling, diffusivity, viscosity, and reaction rate, respectively\.

Table 1:Fixed\-coefficient full\-physics benchmarks and data\-generation parameters\.All fixed\-coefficient trajectories are generated with the pseudo\-spectral reference solver used in APEBench\[[30](https://arxiv.org/html/2608.24049#bib.bib30)\]\. The coefficients, grids, and frame intervals follow the corresponding configurations of the APEBench scenario and are listed in Table[1](https://arxiv.org/html/2608.24049#S3.T1)\. The train and test splits are fixed for each benchmark and shared by all compared methods\. Models are trained on the first 50 steps of each training trajectory and evaluated by 200\-step autoregressive rollout under held\-out initial conditions\. This protocol tests both generalization to unseen initial conditions and extrapolation beyond the temporal window used for training\. The advection–Allen–Cahn benchmark uses a shorter frame interval,Δ​t=0\.1\\Delta t=0\.1, corresponding to roughly ten frames per reaction time𝒪⁡\(1/r\)\\mathcal\{O\}\(1/r\)\.

For the coefficient\-unknown regime, we use a parameterized advection–diffusion setting generated with the same pseudo\-spectral reference solver, with governing form

∂tu=−𝐯⋅∇u\+D∇2u\.\\partial\_\{t\}u=\-\\mathbf\{v\}\\cdot\\nabla u\+D\\nabla^\{2\}u\.\(15\)The coefficient vector𝝀=\(𝐯,D\)\\boldsymbol\{\\lambda\}=\(\\mathbf\{v\},D\)is withheld and estimated from training trajectories before rollout training\. The identification procedure and its accuracy are presented together with the partial\-physics results in Section[4\.5](https://arxiv.org/html/2608.24049#S4.SS5)\. The identified coefficients are then fixed for the corresponding trajectory and used during rollout training and evaluation\. The trajectory\-dependent coefficients are also supplied to the neural\-operator component through the coefficient\-conditioning mechanism\. The coefficient ranges, calibration window, and sampling details are reported in[D](https://arxiv.org/html/2608.24049#A4)\.

### 3\.2Model instantiation and training protocol

The proposed surrogate is trained in two stages\. In Stage 1, the FFGS encoder is trained on individual snapshots using the reconstruction objective in Eq\. \([10](https://arxiv.org/html/2608.24049#S2.E10)\) and is then frozen\. In Stage 2, the frozen representation is used inside the embedded\-physics component, and only the neural\-operator parameters are optimized using the composite rollout objective in Eq\. \([12](https://arxiv.org/html/2608.24049#S2.E12)\)\. All experiments use a rollout length ofH=5H=5\.

Unless otherwise stated, the FFGS representation uses a regular lattice withA=20A=20sites per spatial dimension and local rendering with half\-widthWloc=4W\_\{\\mathrm\{loc\}\}=4anchor cells\. Implementation details of the FFGS representation, including Gaussian parameterization, normalization, periodic boundary treatment, local rendering, and truncation, are specified in[A](https://arxiv.org/html/2608.24049#A1)\. The spectral filter and time\-integration details of the embedded\-physics component are specified in[B](https://arxiv.org/html/2608.24049#A2)\.

The neural\-operator component is instantiated using theClassicFNOimplementation ofpdequinox\. Fixed\-coefficient experiments use an unconditioned neural operator, while coefficient\-varying or coefficient\-unknown experiments use a coefficient\-conditioned variant\. The neural\-operator architecture, conditioning mechanism, optimizer, learning\-rate schedule, batch sizes, training lengths, and architecture\-specific hyperparameters are listed in[C](https://arxiv.org/html/2608.24049#A3)\.

The train and test splits are fixed within each benchmark and shared by all compared methods\. Unless otherwise stated, reported standard deviations are computed across the 30 held\-out test trajectories\.

### 3\.3Baselines

The proposed model is evaluated as a data\-driven next\-step surrogate: it is trained from trajectory data, advances the state one frame at a time, and is deployed autoregressively\. Its natural baselines are therefore purely data\-driven neural time steppers with the same input–output signature and training protocol\. We compare against FNO\[[31](https://arxiv.org/html/2608.24049#bib.bib31)\], U\-Net\[[32](https://arxiv.org/html/2608.24049#bib.bib32)\], and ResNet\[[33](https://arxiv.org/html/2608.24049#bib.bib33)\]in their APEBench configurations\[[30](https://arxiv.org/html/2608.24049#bib.bib30)\]\.

Parameter budgets are matched at the level of the trainable dynamics map\. The proposed model uses a neural\-operator component with a trainable budget comparable to the purely data\-driven baselines\. The FFGS encoder is trained separately for representation reconstruction, frozen before rollout training, and accounted for separately\.

### 3\.4Evaluation metrics

The primary field\-level metric is the relativeℓ2\\ell\_\{2\}error at lead timekk,

rL2⁡\(k\)=‖𝐮pred​\(t\+k​Δ​t\)−𝐮⁡\(t\+k​Δ​t\)‖2‖𝐮⁡\(t\+k​Δ​t\)‖2\.\\mathrm\{rL2\}\(k\)=\\frac\{\\left\\\|\\mathbf\{u\}^\{\\mathrm\{pred\}\}\(t\+k\\Delta t\)\-\\mathbf\{u\}\(t\+k\\Delta t\)\\right\\\|\_\{2\}\}\{\\left\\\|\\mathbf\{u\}\(t\+k\\Delta t\)\\right\\\|\_\{2\}\}\.\(16\)For rollout curves,rL2⁡\(k\)\\mathrm\{rL2\}\(k\)is averaged over the 30 held\-out test trajectories at each lead time\. For the summary tables, the error is first averaged over the 200\-step rollout of each trajectory and is then reported as the mean±\\pmone standard deviation across the 30 trajectory\-level averages\.

To assess spectral fidelity, we compute the radial energy spectrum of each rollout frame\. For a fieldu⁡\(𝐱\)u\(\\mathbf\{x\}\)with Fourier coefficientsu^c​\(𝐪\)\\widehat\{u\}\_\{c\}\(\\mathbf\{q\}\), the radial spectrum at integer wavenumberκ\\kappais defined as

E⁡\(κ\)=∑κ−12≤\|𝐪\|<κ\+12∑c=1C\|u^c​\(𝐪\)\|2,E\(\\kappa\)=\\sum\_\{\\kappa\-\\frac\{1\}\{2\}\\leq\|\\mathbf\{q\}\|<\\kappa\+\\frac\{1\}\{2\}\}\\;\\sum\_\{c=1\}^\{C\}\\left\|\\widehat\{u\}\_\{c\}\(\\mathbf\{q\}\)\\right\|^\{2\},\(17\)where𝐪\\mathbf\{q\}denotes the spatial Fourier mode, and the channel sum is omitted for scalar fields\.

The radial spectrum is then averaged over all rollout frames and held\-out test trajectories,

E¯​\(κ\)=1Ntest​Troll​∑n=1Ntest∑τ=1TrollEn,τ​\(κ\),Ntest=30,Troll=200,\\overline\{E\}\(\\kappa\)=\\frac\{1\}\{N\_\{\\mathrm\{test\}\}T\_\{\\mathrm\{roll\}\}\}\\sum\_\{n=1\}^\{N\_\{\\mathrm\{test\}\}\}\\sum\_\{\\tau=1\}^\{T\_\{\\mathrm\{roll\}\}\}E\_\{n,\\tau\}\(\\kappa\),\\qquad N\_\{\\mathrm\{test\}\}=30,\\quad T\_\{\\mathrm\{roll\}\}=200,\(18\)and the spectral error is defined as

Spec\.err=‖E¯pred−E¯GT‖2‖E¯GT‖2\.\\mathrm\{Spec\.\\,err\}=\\frac\{\\left\\\|\\overline\{E\}^\{\\mathrm\{pred\}\}\-\\overline\{E\}^\{\\mathrm\{GT\}\}\\right\\\|\_\{2\}\}\{\\left\\\|\\overline\{E\}^\{\\mathrm\{GT\}\}\\right\\\|\_\{2\}\}\.\(19\)
Theκ=0\\kappa=0mode represents the domain mean rather than nonzero\-wavenumber spectral content\. For the advection, advection–diffusion, and Burgers benchmarks, the reference mean is identically zero and exactly conserved\. Accordingly, Eq\. \([19](https://arxiv.org/html/2608.24049#S3.E19)\) is evaluated overκ≥1\\kappa\\geq 1for these benchmarks\. For the advection–Allen–Cahn benchmark, the order parameter has an evolving, non\-conserved mean, and theκ=0\\kappa=0mode is therefore retained\.

Lower values indicate better performance for both rL2 and spectral error\. For completeness, we also report the peak signal\-to\-noise ratio,

PSNR=10​log10​\(R2MSE\),\\mathrm\{PSNR\}=10\\log\_\{10\}\\left\(\\frac\{R^\{2\}\}\{\\mathrm\{MSE\}\}\\right\),\(20\)whereRRis the global dynamic range of the reference data over the held\-out test set\. PSNR is aggregated in the same manner as rL2: it is first averaged over the 200\-step rollout of each trajectory and is then reported as the mean±\\pmone standard deviation across the 30 trajectory\-level averages\.

## 4Results

### 4\.1Aggregate rollout accuracy

Table[2](https://arxiv.org/html/2608.24049#S4.T2)summarizes aggregate test performance over 200\-step autoregressive rollouts on the five fixed\-coefficient benchmarks under the common protocol of Section[3](https://arxiv.org/html/2608.24049#S3)\. The proposed surrogate achieves the best performance on all benchmarks and all reported metrics\. Compared with the strongest purely data\-driven baseline in each benchmark, selected by mean relativeℓ2\\ell\_\{2\}error, the proposed surrogate reduces rL2 by factors ranging from1\.5×1\.5\\timeson 3D advection–diffusion to2\.2×2\.2\\timeson Burgers, with corresponding PSNR gains of3\.33\.3to6\.66\.6dB\.

Reporting conventions follow Section[3\.4](https://arxiv.org/html/2608.24049#S3.SS4)\. \#Params counts the trainable parameters of the learned dynamics map\. The FFGS encoder is trained separately for representation reconstruction and frozen during rollout training\. Bold entries mark the best value within each benchmark\.

Table 2:Aggregate test\-rollout performance on the five PDE benchmarks\.The strongest purely data\-driven baseline varies across benchmarks\. Specifically, FNO performs best on the 2D linear transport cases, ResNet on the nonlinear cases, and U\-Net on 3D advection–diffusion\. In contrast, the proposed surrogate is the most accurate method on all five benchmarks\. This improvement is not due to a larger learned model: its neural\-operator component has the same input–output signature and a comparable trainable dynamics\-map budget as the purely data\-driven FNO stepper\. Instead, the comparison points to the FFGS\-based embedded\-physics component as the source of the improvement\. The spectral\-energy metric supports the same overall conclusion, with the proposed surrogate achieving the lowest spectral error on all five benchmarks\. Section[4\.4](https://arxiv.org/html/2608.24049#S4.SS4)examines this metric in more detail\.

### 4\.2Long\-horizon autoregressive rollout

Figure 2:Long\-horizon autoregressive rollout on the Burgers benchmark\. Relativeℓ2\\ell\_\{2\}error is plotted on a logarithmic scale for the proposed surrogate and three purely data\-driven baselines\. Panel \(a\) shows rollouts on the training trajectories, whereas panel \(b\) reports rollouts on held\-out trajectories with unseen initial conditions\. The shaded region in panel \(a\) indicates extrapolation beyond the 50\-step training horizon \(t\>50t\>50\)\.Figure[2](https://arxiv.org/html/2608.24049#S4.F2)compares the long\-horizon autoregressive rollout performance of the proposed surrogate with purely data\-driven baselines on the Burgers benchmark\. The 200\-step rollout simultaneously evaluates extrapolation beyond the 50\-step training horizon and generalization to unseen initial conditions\.

Across both the training and held\-out trajectories, the proposed surrogate consistently achieves the lowest relativeℓ2\\ell\_\{2\}error throughout the rollout\. While all methods perform similarly during the early rollout, the performance gap becomes increasingly pronounced beyond the training horizon, indicating substantially slower autoregressive error accumulation and improved long\-term stability under repeated prediction\. This behavior is retained for unseen initial conditions over the full 200\-step horizon\.

Since all methods employ comparable trainable dynamics\-map budgets, the improved rollout performance cannot be attributed to increased learned\-model capacity\. Instead, the results indicate that embedding the FFGS\-based embedded\-physics component yields more stable autoregressive prediction and better generalization to unseen trajectories\.

### 4\.3Field\-level comparison

![Refer to caption](https://arxiv.org/html/2608.24049v1/3Dadvdiff_snapshot.png)Figure 3:Three\-dimensional advection–diffusion rollout for a representative held\-out trajectory at large Péclet number\. Positive and negative scalar structures are visualized using the isosurfacesu=\+0\.4u=\+0\.4andu=−0\.4u=\-0\.4, respectively\. Columns correspond to the ground truth, the proposed method, FNO, UNet, and ResNet, while rows show the initial condition and autoregressive predictions att=50t=50,100100, and200200\. All methods are initialized from the same ground\-truth field att=0t=0, and the viewing direction and isosurface levels are fixed across all panels\.Figure[3](https://arxiv.org/html/2608.24049#S4.F3)provides a qualitative comparison of the three\-dimensional scalar\-field evolution over a long autoregressive rollout\. The paired isosurfaces expose both the spatial transport of the field and the deformation of its positive and negative structures, thereby revealing discrepancies that may not be apparent from domain\-averaged error metrics alone\. Since all models start from the same initial condition, the differences observed at later times arise from the accumulation of rollout errors\.

Att=50t=50, all methods reproduce the dominant transported structures, although deviations in the geometry and placement of smaller isosurface components have already begun to emerge for the purely data\-driven baselines\. These discrepancies become more pronounced att=100t=100, where FNO, UNet, and particularly ResNet exhibit increasing deformation, merging, or disappearance of individual structures relative to the ground truth\. The proposed method maintains closer agreement in the spatial organization, signed structure, and characteristic length scales of the evolving field\.

Byt=200t=200, the rollout has accumulated substantial temporal error for all surrogate models, as expected for this advection\-dominated three\-dimensional problem\. Nevertheless, the proposed method retains the principal positive and negative structures more faithfully than the baselines and exhibits less severe distortion of their overall topology and spatial distribution\. FNO and UNet preserve portions of the large\-scale organization but show greater displacement and morphological deviation, whereas ResNet displays the strongest departure from the reference solution\. This qualitative behavior is consistent with the quantitative performance reported in Table[2](https://arxiv.org/html/2608.24049#S4.T2), supporting the conclusion that embedding the available governing dynamics improves the structural fidelity and long\-horizon stability of the learned rollout\.

### 4\.4Spectral fidelity

Figure 4:Spectral fidelity on the 2D Burgers test set \(ν=10−4\\nu=10^\{\-4\}; 30 held\-out trajectories×\\times200 rollout steps\)\. \(a\) Time\-averaged radial energy spectrumE¯​\(κ\)\\overline\{E\}\(\\kappa\)for the reference solution and each predictor\. \(b\) RatioE¯​\(κ\)/E¯GT​\(κ\)\\overline\{E\}\(\\kappa\)/\\overline\{E\}\_\{\\mathrm\{GT\}\}\(\\kappa\)over the energy\-containing range\. \(c\) Relative error of the time\-averaged spectrum,‖E¯−E¯GT‖2/‖E¯GT‖2\\\|\\overline\{E\}\-\\overline\{E\}\_\{\\mathrm\{GT\}\}\\\|\_\{2\}/\\\|\\overline\{E\}\_\{\\mathrm\{GT\}\}\\\|\_\{2\}, evaluated overκ≥1\\kappa\\geq 1\(Section[3\.4](https://arxiv.org/html/2608.24049#S3.SS4)\)\.Figure[4](https://arxiv.org/html/2608.24049#S4.F4)compares the time\-averaged radial energy spectra of the 2D Burgers rollouts\. As shown in Fig\.[4](https://arxiv.org/html/2608.24049#S4.F4)\(a\), all methods reproduce the overall spectral decay, while the proposed surrogate provides the closest agreement with the reference spectrum over the resolved wavenumber range\. The largest discrepancy occurs at the first nonzero wavenumber, where the purely data\-driven baselines systematically overestimate the spectral energy, whereas the proposed surrogate remains nearly indistinguishable from the reference, as highlighted by the spectral ratio in Fig\.[4](https://arxiv.org/html/2608.24049#S4.F4)\(b\)\. Consequently, the proposed surrogate achieves the smallest spectral error across the entire resolved spectrum, as summarized in Fig\.[4](https://arxiv.org/html/2608.24049#S4.F4)\(c\)\.

The improved spectral fidelity is consistent with the proposed physics\-integrated formulation\. By coupling the neural operator with the FFGS\-based embedded\-physics component, the model better preserves the evolution of resolved\-scale spectral energy during long autoregressive rollouts, thereby suppressing the accumulation of spurious spectral errors\.

### 4\.5Partial\-physics regime

We now evaluate the coefficient\-unknown regime, in which the governing operator form is available but the advection–diffusion coefficients are withheld\. We first identify the coefficients from trajectory data and then examine whether the resulting physics\-integrated surrogate retains its long\-horizon rollout accuracy and spectral fidelity\.

#### 4\.5\.1Coefficient identification

In the coefficient\-unknown setting, the governing operator is available whereas its coefficients are withheld\. Before rollout training, the unknown coefficients are identified from trajectory data using the frozen FFGS representation\. The closed\-form spatial derivatives supplied by FFGS enable the prescribed governing terms to be evaluated directly, yielding a linear regression system for the unknown coefficients\. The coefficient vector is recovered using the least\-squares formulation in Eq\.\([14](https://arxiv.org/html/2608.24049#S2.E14)\)\. Before solving the regression system, the columns of the design matrix are normalized to improve numerical conditioning; the complete calibration procedure is detailed in[D](https://arxiv.org/html/2608.24049#A4)\.

For the advection–diffusion benchmark, the recovered coefficients are highly accurate, with relative errors of0\.34%0\.34\\%for the velocity and2\.84%2\.84\\%for the diffusivity\. The slightly larger diffusivity error is expected in the large\-Péclet\-number regime, where advection dominates the dynamics and the diffusive contribution provides a weaker regression signal\. Nevertheless, both coefficients are recovered with sufficient accuracy to instantiate the embedded\-physics component for the subsequent rollout experiments\.

#### 4\.5\.2Rollout accuracy and spectral fidelity

Figure 5:Partial\-physics regime on the advection–diffusion benchmark\. \(a\) Test\-rollout relativeℓ2\\ell\_\{2\}error for the proposed surrogate and the purely data\-driven baselines at matched and tenfold \(10×10\\times\) trainable dynamics\-map budgets\. \(b\) Time\-averaged energy spectra; the dotted line marks the initial\-condition band edge\.Figures[5](https://arxiv.org/html/2608.24049#S4.F5)\(a,b\) compare the proposed surrogate in the partial\-physics regime with the purely data\-driven steppers\. The rollout error of each baseline grows toward𝒪⁡\(1\)\\mathcal\{O\}\(1\)within the 200\-step horizon; the10×10\\timesvariants delay this degradation but do not eliminate it\. In contrast, the proposed surrogate maintains slower error growth throughout the rollout\. The time\-averaged spectra provide a consistent diagnosis: the baselines accumulate spurious high\-wavenumber energy beyond the reference band, most severely for the10×10\\timesFNO, whereas the proposed surrogate follows the reference spectrum down to the numerical floor\. Thus, even when the physical coefficients are estimated rather than prescribed, the embedded\-physics component continues to improve long\-horizon accuracy and spectral fidelity\. The resulting surrogate remains data\-driven end to end, including the coefficient identification step, while outperforming the purely data\-driven baselines\.

### 4\.6Inference cost

Figure 6:Per\-step inference cost measured as wall\-clock time per autoregressive update\. Values are the minimum over five timed rollouts after warm\-up and are shown on a logarithmic scale\.Figure[6](https://arxiv.org/html/2608.24049#S4.F6)reports wall\-clock inference time per autoregressive step for 2D and 3D advection–diffusion\. The proposed surrogate is slower than FNO and U\-Net because it evaluates the FFGS\-based embedded\-physics component during each update\. However, its cost remains comparable to the slower neural baselines: it is close to ResNet in 2D and faster than ResNet in 3D\. Together with the accuracy results in Table[2](https://arxiv.org/html/2608.24049#S4.T2), this shows an accuracy–cost trade\-off\. The additional computation comes from evaluating the physics\-coupled Gaussian representation inside the update map, not from increasing the trainable neural\-operator component\.

## 5Discussion

### 5\.1Component analysis

The results above show consistent gains over purely data\-driven neural steppers\. We now use ablation studies to identify which parts of the composite update are responsible for these gains\. In the fixed\-coefficient regime, we test the role of the FFGS\-based embedded\-physics component and the low\-pass filter \(Fig\.[7](https://arxiv.org/html/2608.24049#S5.F7)\)\. In the partial\-physics regime, after the coefficients have been identified, we compare an FFGS\-only variant containing only the embedded\-physics component, a coefficient\-conditioned FNO without the embedded\-physics component, and the proposed composite surrogate \(Fig\.[8](https://arxiv.org/html/2608.24049#S5.F8)\)\. This comparison tests whether the embedded\-physics component and the coefficient\-conditioned neural\-operator component are complementary\.

Figure 7:Ablation in the full\-physics regime\. \(a\) Removing the FFGS\-based embedded\-physics component leaves the neural\-operator\-only variant with the same trainable dynamics\-map budget\. \(b\) Removing the low\-pass filter from the embedded\-physics component on the Burgers benchmark\. Both panels report rollout relativeℓ2\\ell\_\{2\}error on a logarithmic scale\.Figure[7](https://arxiv.org/html/2608.24049#S5.F7)shows that both the embedded\-physics component and the low\-pass filter contribute to rollout accuracy\. Removing the FFGS\-based embedded\-physics component reduces the model to a neural\-operator\-only stepper and increases the rollout error on all five fixed\-coefficient benchmarks\. The effect is especially pronounced on the advection–Allen–Cahn benchmark, where the available governing terms carry substantial information about the reaction\-interface dynamics\. This indicates that the gains in Table[2](https://arxiv.org/html/2608.24049#S4.T2)are not produced by the neural\-operator component alone\.

The low\-pass filter has a different role\. On Burgers, removing the filter from the embedded\-physics component leads to a growing error gap over the rollout horizon\. This is consistent with the fact that spatial differentiation amplifies high\-wavenumber reconstruction errors in the Gaussian representation\. The filter therefore acts as a stabilization device for physics evaluation: it suppresses unreliable high\-wavenumber derivative content before the governing terms are evaluated\. The filtered and unfiltered variants use the same neural\-operator component, trainable budget, and rollout loss, so the difference isolates the effect of the filter\.

Figure 8:Component ablation in the partial\-physics regime\. All variants use the same recovered coefficients\. FFGS\-only uses only the embedded\-physics component; FiLM\-FNO supplies the recovered coefficients to a coefficient\-conditioned neural operator without the embedded\-physics component; Ours combines the coefficient\-conditioned neural\-operator component with the FFGS\-based embedded\-physics component\. \(a\) Test\-rollout relativeℓ2\\ell\_\{2\}error\. \(b\) Time\-averaged radial energy spectrum\. \(c\) Mean rollout error over the 200\-step horizon\.The partial\-physics ablation tests whether the embedded\-physics component and the coefficient\-conditioned neural\-operator component are complementary after coefficient identification\. The FFGS\-only variant uses the recovered coefficients inside the governing operator but contains no neural\-operator component\. Its rollout error grows slowly, but it starts from a large error offset because the available physics and the filtered Gaussian representation do not fully reproduce the reference dynamics\. Its spectrum also shows excessive damping of high\-wavenumber content\. The FiLM\-FNO variant uses the same recovered coefficients only as neural\-network conditioning\. It improves short\-time prediction, but without the embedded\-physics component, its rollout drifts over long horizons and accumulates high\-wavenumber energy\. The proposed surrogate combines the two components\. The embedded\-physics component uses the recovered coefficients structurally, while the coefficient\-conditioned neural\-operator component accounts for the part of the update not captured by the available physics\. This combination yields lower rollout error and better spectral fidelity than either component alone\. The comparison shows that exposing the coefficient values to a neural operator is not equivalent to evaluating the corresponding governing terms on the differentiable FFGS representation\.

Taken together, the fixed\-coefficient and partial\-physics ablations show that the embedded\-physics component and the neural\-operator component are complementary\. The embedded\-physics component improves long\-horizon behavior by injecting the available governing structure into the update map, while the neural\-operator component compensates for representation error, filtering effects, and dynamics not fully captured by the available physics\. Removing either component degrades performance\. The recurrence of this pattern across prescribed\-coefficient and coefficient\-unknown settings suggests that the gain comes from the composite update itself, rather than from a benchmark\-specific tuning or a particular coefficient choice\.

### 5\.2Limitations and future directions

The results demonstrate the benefit of embedding available governing structure through a differentiable Gaussian representation, but the present implementation has several limitations\. These limitations concern the fixed capacity of the representation, its specialization to a given solution family, and the manual selection of the spectral cutoff\.

A first limitation is the fixed capacity of the Gaussian representation\. The feed\-forward encoder predicts a fixed number of anisotropic kernels on a regular anchor lattice\. This design makes the representation efficient and amortized, but the spectral range over which it can support accurate physics evaluation is limited by the anchor density, kernel scales, and derivative accuracy\. Sharp fronts or persistent high\-wavenumber content may therefore be reconstructed only approximately, and derivative quantities in those spectral ranges can become unreliable\. The limitation is therefore not only pointwise reconstruction error, but also the accuracy of the spatial derivatives used by the embedded operator\. This motivates adaptive Gaussian representations, such as variable kernel budgets or anchor refinement in regions of high local complexity\. Such adaptivity could increase the usable spectral bandwidth of the representation and reduce the reliance on low\-pass filtering during embedded\-physics evaluation\.

A second limitation is specialization to a given equation family\. In the present experiments, the FFGS encoder is trained on trajectories from the same equation family used for evaluation\. Its reconstruction quality, derivative accuracy, and filter choice are therefore tied to that solution family\. Extending the method across multiple PDEs would require a shared representation whose rendered fields and closed\-form derivatives remain reliable across different operators, parameter regimes, and solution morphologies\. This requirement is stronger than interpolation accuracy alone, because the derivatives are evaluated inside the governing operator\. A possible direction is to condition the encoder and the neural\-operator component on the governing coefficients, nondimensional parameters, or an operator descriptor while retaining a common FFGS interface\.

A third limitation is the manual selection of the spectral cutoff\. In this work, the cutoff is chosen per benchmark to retain the spectral range in which the FFGS reconstruction and its derivatives are reliable, while suppressing high\-wavenumber derivative errors\. This choice is effective but remains external to the model\. A more systematic strategy would estimate the reliable spectral range directly from diagnostics of the reconstructed field and its derivatives\. Since the FFGS reconstruction is explicit, quantities such as spectral decay, reconstruction error, first\-derivative error, and second\-derivative error can be measured against reference fields \([F](https://arxiv.org/html/2608.24049#A6)\)\. These diagnostics could be used to select the cutoff automatically and, in future extensions, to adapt it to the field, governing operator, and time step\.

## 6Conclusion

This work introduced a physics\-integrated neural\-operator framework that incorporates available governing\-equation knowledge at the level of the field representation\. Rather than imposing the physics through an auxiliary residual loss or redesigning the neural\-operator architecture for a particular PDE, the proposed approach retains a generic neural\-operator backbone and uses feed\-forward Gaussian splatting \(FFGS\) as a continuous interface between sampled solution fields and governing operators\. FFGS maps each discrete state to a compact continuous Gaussian representation in a single forward pass, whose analytic structure provides direct access to the spatial differential quantities required by the governing equations\. The resulting surrogate combines a neural\-operator component with an embedded\-physics component evaluated directly on this representation\.

Across five two\- and three\-dimensional PDE benchmarks spanning transport, diffusion, nonlinear self\-advection, and reaction dynamics, the proposed framework consistently improves long\-horizon autoregressive accuracy and spectral fidelity over purely data\-driven neural steppers with comparable trainable dynamics\-map budgets\. The gains persist beyond the temporal window used for training and on held\-out initial conditions, while the ablation studies show that both the embedded\-physics component and the spectral filtering of reconstructed differential quantities contribute to the improved rollout behavior\. The same representation also extends naturally to the partial\-physics setting: when the governing form is known but its coefficients are unavailable, the analytic derivatives supplied by FFGS enable the coefficients to be identified from trajectory data and subsequently used within the embedded\-physics component\. The resulting model retains most of the accuracy and spectral\-fidelity gains observed when the coefficients are prescribed\.

These results support representation\-level physics integration as a practical route for coupling neural operators with fully or partially known PDE structure\. More broadly, the role of the learned representation extends beyond reconstructing the state: it provides the continuous physical field on which governing operators can be evaluated before their contribution is incorporated into the learned time\-evolution map\. The present implementation remains limited by the fixed capacity of the Gaussian representation, its specialization to a given solution family, and the use of manually selected spectral cutoffs\. Addressing these limitations through adaptive representations, cross\-equation conditioning, and data\-driven selection of the reliable spectral range provides a natural path toward applying the framework to more complex multiscale and multiphysics systems\.

## References

- \[1\]L\. Lu, P\. Jin, G\. Pang, Z\. Zhang, G\. E\. Karniadakis,[Learning nonlinear operators via DeepONet based on the universal approximation theorem of operators](https://www.nature.com/articles/s42256-021-00302-5)3 \(3\) 218–229\.[doi:10\.1038/s42256\-021\-00302\-5](https://doi.org/10.1038/s42256-021-00302-5)\. URL[https://www\.nature\.com/articles/s42256\-021\-00302\-5](https://www.nature.com/articles/s42256-021-00302-5)
- \[2\]N\. Kovachki, Z\. Li, B\. Liu, K\. Azizzadenesheli, K\. Bhattacharya, A\. Stuart, A\. Anandkumar,[Neural Operator: Learning Maps Between Function Spaces](http://arxiv.org/abs/2108.08481), arXiv:2108\.08481 \[cs\.LG\] \(May 2024\)\.[doi:10\.5555/3648699\.3648788](https://doi.org/10.5555/3648699.3648788)\. URL[http://arxiv\.org/abs/2108\.08481](http://arxiv.org/abs/2108.08481)
- \[3\]J\. Pathak, S\. Subramanian, P\. Harrington, S\. Raja, A\. Chattopadhyay, M\. Mardani, T\. Kurth, D\. Hall, Z\. Li, K\. Azizzadenesheli, et al\., Fourcastnet: A global data\-driven high\-resolution weather model using adaptive fourier neural operators, arXiv preprint arXiv:2202\.11214 \(2022\)\.
- \[4\]K\. Bi, L\. Xie, H\. Zhang, X\. Chen, X\. Gu, Q\. Tian, Accurate medium\-range global weather forecasting with 3d neural networks, Nature 619 \(7970\) \(2023\) 533–538\.
- \[5\]Z\. Li, W\. Peng, Z\. Yuan, J\. Wang,[Fourier neural operator approach to large eddy simulation of three\-dimensional turbulence](https://www.sciencedirect.com/science/article/pii/S2095034922000691)12 \(6\) 100389\.[doi:10\.1016/j\.taml\.2022\.100389](https://doi.org/10.1016/j.taml.2022.100389)\. URL[https://www\.sciencedirect\.com/science/article/pii/S2095034922000691](https://www.sciencedirect.com/science/article/pii/S2095034922000691)
- \[6\]Y\. Wang, Z\. Li, Z\. Yuan, W\. Peng, T\. Liu, J\. Wang,[Prediction of turbulent channel flow using Fourier neural operator\-based machine\-learning strategy](https://link.aps.org/doi/10.1103/PhysRevFluids.9.084604)9 \(8\) 084604\.[doi:10\.1103/PhysRevFluids\.9\.084604](https://doi.org/10.1103/PhysRevFluids.9.084604)\. URL[https://link\.aps\.org/doi/10\.1103/PhysRevFluids\.9\.084604](https://link.aps.org/doi/10.1103/PhysRevFluids.9.084604)
- \[7\]P\. Lippe, B\. Veeling, P\. Perdikaris, R\. Turner, J\. Brandstetter,[PDE\-Refiner: Achieving Accurate Long Rollouts with Neural PDE Solvers](https://proceedings.neurips.cc/paper_files/paper/2023/hash/d529b943af3dba734f8a7d49efcb6d09-Abstract-Conference.html), in: Advances in Neural Information Processing Systems, Vol\. 36, Curran Associates, Inc\., pp\. 67398–67433\. URL[https://proceedings\.neurips\.cc/paper\_files/paper/2023/hash/d529b943af3dba734f8a7d49efcb6d09\-Abstract\-Conference\.html](https://proceedings.neurips.cc/paper_files/paper/2023/hash/d529b943af3dba734f8a7d49efcb6d09-Abstract-Conference.html)
- \[8\]M\. Takamoto, T\. Praditia, R\. Leiteritz, D\. MacKinlay, F\. Alesiani, D\. Pflüger, M\. Niepert,[PDEBENCH: An Extensive Benchmark for Scientific Machine Learning](http://arxiv.org/abs/2210.07182)\.[arXiv:2210\.07182](http://arxiv.org/abs/2210.07182),[doi:10\.48550/arXiv\.2210\.07182](https://doi.org/10.48550/arXiv.2210.07182)\. URL[http://arxiv\.org/abs/2210\.07182](http://arxiv.org/abs/2210.07182)
- \[9\]J\. A\. L\. Benitez, T\. Furuya, F\. Faucher, A\. Kratsios, X\. Tricoche, M\. V\. de Hoop, Out\-of\-distributional risk bounds for neural operators with applications to the helmholtz equation, Journal of Computational Physics 513 \(2024\) 113168\.
- \[10\]S\. Wang, H\. Wang, P\. Perdikaris,[Learning the solution operator of parametric partial differential equations with physics\-informed DeepONets](https://www.science.org/doi/10.1126/sciadv.abi8605)7 \(40\) eabi8605\.[doi:10\.1126/sciadv\.abi8605](https://doi.org/10.1126/sciadv.abi8605)\. URL[https://www\.science\.org/doi/10\.1126/sciadv\.abi8605](https://www.science.org/doi/10.1126/sciadv.abi8605)
- \[11\]Z\. Li, H\. Zheng, N\. Kovachki, D\. Jin, H\. Chen, B\. Liu, K\. Azizzadenesheli, A\. Anandkumar,[Physics\-Informed Neural Operator for Learning Partial Differential Equations](http://arxiv.org/abs/2111.03794), arXiv:2111\.03794 \[cs\.LG\] \(Jul\. 2023\)\.[doi:10\.48550/arXiv\.2111\.03794](https://doi.org/10.48550/arXiv.2111.03794)\. URL[http://arxiv\.org/abs/2111\.03794](http://arxiv.org/abs/2111.03794)
- \[12\]M\. Raissi, P\. Perdikaris, G\. E\. Karniadakis,[Physics\-informed neural networks: A deep learning framework for solving forward and inverse problems involving nonlinear partial differential equations](https://www.sciencedirect.com/science/article/pii/S0021999118307125)378 686–707\.[doi:10\.1016/j\.jcp\.2018\.10\.045](https://doi.org/10.1016/j.jcp.2018.10.045)\. URL[https://www\.sciencedirect\.com/science/article/pii/S0021999118307125](https://www.sciencedirect.com/science/article/pii/S0021999118307125)
- \[13\]S\. Wang, X\. Yu, P\. Perdikaris,[When and why PINNs fail to train: A neural tangent kernel perspective](https://www.sciencedirect.com/science/article/pii/S002199912100663X)449 110768\.[doi:10\.1016/j\.jcp\.2021\.110768](https://doi.org/10.1016/j.jcp.2021.110768)\. URL[https://www\.sciencedirect\.com/science/article/pii/S002199912100663X](https://www.sciencedirect.com/science/article/pii/S002199912100663X)
- \[14\]A\. Krishnapriyan, A\. Gholami, S\. Zhe, R\. Kirby, M\. Mahoney,[Characterizing possible failure modes in physics\-informed neural networks](https://proceedings.neurips.cc/paper/2021/hash/df438e5206f31600e6ae4af72f2725f1-Abstract.html), in: Advances in Neural Information Processing Systems, Vol\. 34, Curran Associates, Inc\., pp\. 26548–26560\. URL[https://proceedings\.neurips\.cc/paper/2021/hash/df438e5206f31600e6ae4af72f2725f1\-Abstract\.html](https://proceedings.neurips.cc/paper/2021/hash/df438e5206f31600e6ae4af72f2725f1-Abstract.html)
- \[15\]N\. Liu, Y\. Fan, X\. Zeng, M\. Klöwer, L\. Zhang, Y\. Yu,[Harnessing the power of neural operators with automatically encoded conservation laws](https://proceedings.mlr.press/v235/liu24p.html), in: R\. Salakhutdinov, Z\. Kolter, K\. Heller, A\. Weller, N\. Oliver, J\. Scarlett, F\. Berkenkamp \(Eds\.\), Proceedings of the 41st International Conference on Machine Learning, Vol\. 235 of Proceedings of Machine Learning Research, PMLR, 2024, pp\. 30965–30997\. URL[https://proceedings\.mlr\.press/v235/liu24p\.html](https://proceedings.mlr.press/v235/liu24p.html)
- \[16\]N\. Liu, Y\. Yu, H\. You, N\. Tatikola,[INO: Invariant Neural Operators for Learning Complex Physical Systems with Momentum Conservation](http://arxiv.org/abs/2212.14365)\.[arXiv:2212\.14365](http://arxiv.org/abs/2212.14365),[doi:10\.48550/arXiv\.2212\.14365](https://doi.org/10.48550/arXiv.2212.14365)\. URL[http://arxiv\.org/abs/2212\.14365](http://arxiv.org/abs/2212.14365)
- \[17\]J\. Helwig, X\. Zhang, C\. Fu, J\. Kurtin, S\. Wojtowytsch, S\. Ji,[Group Equivariant Fourier Neural Operators for Partial Differential Equations](http://arxiv.org/abs/2306.05697)\.[arXiv:2306\.05697](http://arxiv.org/abs/2306.05697),[doi:10\.48550/arXiv\.2306\.05697](https://doi.org/10.48550/arXiv.2306.05697)\. URL[http://arxiv\.org/abs/2306\.05697](http://arxiv.org/abs/2306.05697)
- \[18\]X\. Hu, Q\. Ma, P\. Zhao, X\. Wang, Physics\-aware neural operator for high\-fidelity fluid dynamics modeling with geometric and spectral priors, Physics of Fluids 37 \(11\) \(2025\)\.
- \[19\]M\. Liu\-Schiaffini, J\. Berner, B\. Bonev, T\. Kurth, K\. Azizzadenesheli, A\. Anandkumar,[Neural Operators with Localized Integral and Differential Kernels](http://arxiv.org/abs/2402.16845)\.[arXiv:2402\.16845](http://arxiv.org/abs/2402.16845),[doi:10\.48550/arXiv\.2402\.16845](https://doi.org/10.48550/arXiv.2402.16845)\. URL[http://arxiv\.org/abs/2402\.16845](http://arxiv.org/abs/2402.16845)
- \[20\]E\. J\. Kansa,[Multiquadrics—A scattered data approximation scheme with applications to computational fluid\-dynamics—II solutions to parabolic, hyperbolic and elliptic partial differential equations](https://www.sciencedirect.com/science/article/pii/089812219090271K)19 \(8\) 147–161\.[doi:10\.1016/0898\-1221\(90\)90271\-K](https://doi.org/10.1016/0898-1221(90)90271-K)\. URL[https://www\.sciencedirect\.com/science/article/pii/089812219090271K](https://www.sciencedirect.com/science/article/pii/089812219090271K)
- \[21\]C\. Franke, R\. Schaback,[Solving partial differential equations by collocation using radial basis functions](https://www.sciencedirect.com/science/article/pii/S0096300397101047)93 \(1\) 73–82\.[doi:10\.1016/S0096\-3003\(97\)10104\-7](https://doi.org/10.1016/S0096-3003(97)10104-7)\. URL[https://www\.sciencedirect\.com/science/article/pii/S0096300397101047](https://www.sciencedirect.com/science/article/pii/S0096300397101047)
- \[22\]B\. Kerbl, G\. Kopanas, T\. Leimkühler, G\. Drettakis,[3d Gaussian Splatting for Real\-Time Radiance Field Rendering](https://dl.acm.org/doi/10.1145/3592433), ACM Transactions on Graphics 42 \(4\) \(2023\) 139:1–139:14\.[doi:10\.1145/3592433](https://doi.org/10.1145/3592433)\. URL[https://dl\.acm\.org/doi/10\.1145/3592433](https://dl.acm.org/doi/10.1145/3592433)
- \[23\]D\. V\. Shenoy, S\. H\. Frankel, Gaussian field representations for turbulent flow: Compression, scale separation, and physical fidelity, Computers & Fluids \(2026\) 107202\.
- \[24\]Z\. Li, Y\. Feng, Z\. Lai, W\. Wang, From basis to basis: Gaussian particle representation for interpretable pde operators, arXiv preprint arXiv:2602\.21551 \(2026\)\.
- \[25\]J\. Xing, B\. Wang, M\. Chu, B\. Chen,[Gaussian Fluids: A Grid\-Free Fluid Solver based on Gaussian Spatial Representation](http://arxiv.org/abs/2405.18133), arXiv:2405\.18133 \[cs\.GR\] \(Jul\. 2025\)\.[doi:10\.48550/arXiv\.2405\.18133](https://doi.org/10.48550/arXiv.2405.18133)\. URL[http://arxiv\.org/abs/2405\.18133](http://arxiv.org/abs/2405.18133)
- \[26\]D\. Charatan, S\. Li, A\. Tagliasacchi, V\. Sitzmann,[pixelSplat: 3D Gaussian Splats from Image Pairs for Scalable Generalizable 3D Reconstruction](http://arxiv.org/abs/2312.12337)\.[arXiv:2312\.12337](http://arxiv.org/abs/2312.12337),[doi:10\.48550/arXiv\.2312\.12337](https://doi.org/10.48550/arXiv.2312.12337)\. URL[http://arxiv\.org/abs/2312\.12337](http://arxiv.org/abs/2312.12337)
- \[27\]W\. Wang, Y\. Chen, Z\. Zhang, H\. Liu, H\. Wang, Z\. Feng, W\. Qin, F\. Chen, Z\. Zhu, D\. Y\. Chen, B\. Zhuang,[VolSplat: Rethinking Feed\-Forward 3D Gaussian Splatting with Voxel\-Aligned Prediction](http://arxiv.org/abs/2509.19297)\.[arXiv:2509\.19297](http://arxiv.org/abs/2509.19297),[doi:10\.48550/arXiv\.2509\.19297](https://doi.org/10.48550/arXiv.2509.19297)\. URL[http://arxiv\.org/abs/2509\.19297](http://arxiv.org/abs/2509.19297)
- \[28\]Z\. Zhang, X\. Meng, K\. Wu, W\. Ding,[Sparsesplat: Towards applicable feed\-forward 3d gaussian splatting with pixel\-unaligned prediction](https://arxiv.org/abs/2604.03069)\(2026\)\.[arXiv:2604\.03069](http://arxiv.org/abs/2604.03069)\. URL[https://arxiv\.org/abs/2604\.03069](https://arxiv.org/abs/2604.03069)
- \[29\]H\. Huang, M\. Li, Z\. Gao, X\. Zhou, X\. Huang, X\. Sun, Fluidsplat: Reconstructing physical fields from sparse sensors via gaussian primitives, arXiv preprint arXiv:2605\.18866 \(2026\)\.
- \[30\]F\. Koehler, S\. Niedermayr, R\. Westermann, N\. Thuerey,[APEBench: A benchmark for autoregressive neural emulators of pdes](https://proceedings.neurips.cc/paper_files/paper/2024/file/d9875ebcf74bccdc5076acab0dbee62c-Paper-Datasets_and_Benchmarks_Track.pdf), in: A\. Globerson, L\. Mackey, D\. Belgrave, A\. Fan, U\. Paquet, J\. Tomczak, C\. Zhang \(Eds\.\), Advances in Neural Information Processing Systems, Vol\. 37, Curran Associates, Inc\., pp\. 120252–120310\.[doi:10\.52202/079017\-3822](https://doi.org/10.52202/079017-3822)\. URL[https://proceedings\.neurips\.cc/paper\_files/paper/2024/file/d9875ebcf74bccdc5076acab0dbee62c\-Paper\-Datasets\_and\_Benchmarks\_Track\.pdf](https://proceedings.neurips.cc/paper_files/paper/2024/file/d9875ebcf74bccdc5076acab0dbee62c-Paper-Datasets_and_Benchmarks_Track.pdf)
- \[31\]Z\. Li, N\. Kovachki, K\. Azizzadenesheli, B\. Liu, K\. Bhattacharya, A\. Stuart, A\. Anandkumar,[Fourier Neural Operator for Parametric Partial Differential Equations](http://arxiv.org/abs/2010.08895), arXiv:2010\.08895 \[cs, math\] \(May 2021\)\. URL[http://arxiv\.org/abs/2010\.08895](http://arxiv.org/abs/2010.08895)
- \[32\]O\. Ronneberger, P\. Fischer, T\. Brox,[U\-Net: Convolutional Networks for Biomedical Image Segmentation](https://arxiv.org/abs/1505.04597v1)\. URL[https://arxiv\.org/abs/1505\.04597v1](https://arxiv.org/abs/1505.04597v1)
- \[33\]K\. He, X\. Zhang, S\. Ren, J\. Sun,[Deep Residual Learning for Image Recognition](http://arxiv.org/abs/1512.03385)\.[arXiv:1512\.03385](http://arxiv.org/abs/1512.03385),[doi:10\.48550/arXiv\.1512\.03385](https://doi.org/10.48550/arXiv.1512.03385)\. URL[http://arxiv\.org/abs/1512\.03385](http://arxiv.org/abs/1512.03385)

## Appendix AFFGS architecture and implementation

This appendix summarizes the implementation details of the feed\-forward Gaussian splatting \(FFGS\) representation used throughout this work, including the encoder architecture, Gaussian parameterization, normalization, periodic rendering, and local rendering\.

### A\.1Encoder architecture

The FFGS encoder maps an input field of sizeNd×CN^\{d\}\\times Cto a Gaussian representation defined on an auxiliary anchor lattice containingAdA^\{d\}cells\. Unless otherwise stated, all experiments useA=20A=20anchors per spatial dimension\. The encoder is implemented as a four\-level U\-Net with periodic padding\. A1×11\\times 1prediction head produces one Gaussian primitive for each anchor cell by predicting the center offset, logarithmic scales, orientation, and channel\-wise amplitude\. The covariance matrix is subsequently assembled from the predicted scales and orientation\.

The encoder contains approximately 451k trainable parameters for the two\-dimensional benchmarks and 1\.33M parameters for the three\-dimensional benchmark\.

### A\.2Field normalization

Each input frame is normalized independently for every physical channel\. For each physical channel, the spatial mean and standard deviation are

mt=⟨ut⟩Ω,st=stdΩ⁡\(ut\),m\_\{t\}=\\langle u\_\{t\}\\rangle\_\{\\Omega\},\\qquad s\_\{t\}=\\operatorname\{std\}\_\{\\Omega\}\(u\_\{t\}\),\(21\)and the normalized field is

u¯t=ut−mtst\+ε,\\bar\{u\}\_\{t\}=\\frac\{u\_\{t\}\-m\_\{t\}\}\{s\_\{t\}\+\\varepsilon\},\(22\)whereε=10−6\\varepsilon=10^\{\-6\}\.

The encoder operates on the normalized field\. After rendering, the reconstructed field is restored to the original physical scale,

u^t=\(st\+ε\)​u¯^t\+mt\.\\hat\{u\}\_\{t\}=\(s\_\{t\}\+\\varepsilon\)\\hat\{\\bar\{u\}\}\_\{t\}\+m\_\{t\}\.\(23\)
The reconstructed derivatives are transformed accordingly,

∇u^t=\(st\+ε\)∇u¯^t,∇2u^t=\(st\+ε\)∇2u¯^t\.\\nabla\\hat\{u\}\_\{t\}=\(s\_\{t\}\+\\varepsilon\)\\nabla\\hat\{\\bar\{u\}\}\_\{t\},\\qquad\\nabla^\{2\}\\hat\{u\}\_\{t\}=\(s\_\{t\}\+\\varepsilon\)\\nabla^\{2\}\\hat\{\\bar\{u\}\}\_\{t\}\.\(24\)

### A\.3Gaussian parameterization

Each Gaussian primitive is associated with one anchor cell\. Leth=1/Ah=1/Adenote the anchor spacing\. The Gaussian center is represented as

Δ​𝝁i=h2​tanh⁡\(𝝁~i\),𝝁i=\(𝐱ci\+Δ​𝝁i\)mod1,\\Delta\\boldsymbol\{\\mu\}\_\{i\}=\\frac\{h\}\{2\}\\tanh\(\\widetilde\{\\boldsymbol\{\\mu\}\}\_\{i\}\),\\qquad\\boldsymbol\{\\mu\}\_\{i\}=\\left\(\\mathbf\{x\}\_\{c\_\{i\}\}\+\\Delta\\boldsymbol\{\\mu\}\_\{i\}\\right\)\\bmod 1,\(25\)where𝐱ci\\mathbf\{x\}\_\{c\_\{i\}\}denotes the anchor position\. The bounded offset keeps every primitive associated with its anchor cell while allowing sub\-cell positioning\.

The covariance matrix is parameterized as

Σi=Ri​diag⁡\(𝝈i2\)​Ri⊤,\\Sigma\_\{i\}=R\_\{i\}\\operatorname\{diag\}\(\\boldsymbol\{\\sigma\}\_\{i\}^\{2\}\)R\_\{i\}^\{\\top\},\(26\)where the principal scales are obtained by exponentiating the predicted logarithmic scales and clipping them to the interval\[ρmin,ρmax\]\[\\rho\_\{\\min\},\\rho\_\{\\max\}\], withρmin=0\.008\\rho\_\{\\min\}=0\.008andρmax=0\.25\\rho\_\{\\max\}=0\.25\. Rotations are represented by one angle in two dimensions and by normalized unit quaternions in three dimensions\. Each primitive additionally predicts a channel\-wise amplitude vector𝜶i∈ℝC\\boldsymbol\{\\alpha\}\_\{i\}\\in\\mathbb\{R\}^\{C\}\.

### A\.4Periodic and local rendering

All experiments are performed on the periodic domain\[0,1\)d\[0,1\)^\{d\}\. Periodicity is enforced using the minimum\-image convention,

𝐫i​\(𝐱\)=𝐱−𝝁i−round⁡\(𝐱−𝝁i\),\\mathbf\{r\}\_\{i\}\(\\mathbf\{x\}\)=\\mathbf\{x\}\-\\boldsymbol\{\\mu\}\_\{i\}\-\\operatorname\{round\}\\left\(\\mathbf\{x\}\-\\boldsymbol\{\\mu\}\_\{i\}\\right\),\(27\)where the rounding operator is applied independently to each spatial coordinate\. This periodic displacement is used consistently in Gaussian evaluation and analytical differentiation\.

To reduce computational cost, rendering is restricted to a fixed local candidate window containing\(2​Wloc\+1\)d\(2W\_\{\\mathrm\{loc\}\}\+1\)^\{d\}neighboring anchor cells, withWloc=4W\_\{\\mathrm\{loc\}\}=4in all experiments\. Thus, each query uses a9d9^\{d\}local neighborhood\. Every Gaussian primitive within this window is included in the rendering sum, with no additional distance\-based culling\. This fixed\-window strategy exploits the rapid spatial decay of the Gaussian kernels and is applied consistently during both representation training and embedded\-physics evaluation\.

To quantify the approximation introduced by the finite rendering window, Table[3](https://arxiv.org/html/2608.24049#A1.T3)compares theWloc=4W\_\{\\mathrm\{loc\}\}=4local rendering used in this work against dense summation over all Gaussian primitives\.

Table 3:Comparison of theWloc=4W\_\{\\mathrm\{loc\}\}=4local rendering used in this work with dense summation over all Gaussian primitives\. RelativeL2L\_\{2\}discrepancies are reported for the reconstructed field and its spatial derivatives, together with the rendering speedup\.Across the considered benchmarks, theWloc=4W\_\{\\mathrm\{loc\}\}=4local rendering closely reproduces the dense evaluation for both the reconstructed field and its spatial derivatives, while providing a13\.913\.9–14\.6×14\.6\\timesrendering speedup\.

## Appendix BFiltering and time integration

The frozen FFGS representation supplies the reconstructed quantities and analytical spatial derivatives required by the available governing operator\. State\-dependent algebraic terms are evaluated directly on the current intermediate state, whereas spatial differential terms are evaluated from the closed\-form FFGS derivatives\. Before entering the embedded\-physics operator, the reconstructed derivative quantities are filtered onto a retained spectral band using the low\-pass operator introduced in Section[2\.2](https://arxiv.org/html/2608.24049#S2.SS2)\.

For a grid quantityq⁡\(𝐱\)q\(\\mathbf\{x\}\), the filter is applied in Fourier space as

Πκ​q^​\(𝐤\)=wκ​\(‖𝐤‖2\)​q^​\(𝐤\),\\widehat\{\\Pi\_\{\\kappa\}q\}\(\\mathbf\{k\}\)=w\_\{\\kappa\}\\\!\\left\(\\\|\\mathbf\{k\}\\\|\_\{2\}\\right\)\\hat\{q\}\(\\mathbf\{k\}\),\(28\)whereq^\\hat\{q\}denotes the Fourier coefficients,‖𝐤‖2\\\|\\mathbf\{k\}\\\|\_\{2\}is the radial wavenumber, andwκw\_\{\\kappa\}is the spectral filter\. For sharp filtering,

wκ​\(k\)=\{1,k≤κ,0,k\>κ,w\_\{\\kappa\}\(k\)=\\begin\{cases\}1,&k\\leq\\kappa,\\\\ 0,&k\>\\kappa,\\end\{cases\}\(29\)whereas smooth filtering uses

wκ​\(k\)=\{1,k≤κ,exp⁡\[−\(k−κWf\)2\],k\>κ\.w\_\{\\kappa\}\(k\)=\\begin\{cases\}1,&k\\leq\\kappa,\\\\\[5\.69054pt\] \\exp\\\!\\left\[\-\\left\(\\dfrac\{k\-\\kappa\}\{W\_\{f\}\}\\right\)^\{2\}\\right\],&k\>\\kappa\.\\end\{cases\}\(30\)
The filter suppresses high\-wavenumber derivative errors from the finite\-capacity Gaussian representation, restricts the bandwidth entering the explicit time integrator, and removes out\-of\-band modes generated by nonlinear products\. For the Burgers and advection–Allen–Cahn equations, the assembled right\-hand side is therefore spectrally filtered again after the nonlinear terms have been evaluated\. This additional filter is redundant for the linear benchmarks because the spectral filter commutes with the corresponding linear combinations\.

The cutoff parameters used in the fixed\-coefficient benchmarks are summarized in Table[4](https://arxiv.org/html/2608.24049#A2.T4)\. The cutoff is chosen to balance retention of the dynamically relevant modes, the reliability of the FFGS derivatives, and the stability of the explicit time integration\. Since all initial conditions are generated from truncated Fourier series satisfying\|ki\|≤5\|k\_\{i\}\|\\leq 5, the linear transport cases use cutoffs close to the initial\-condition band, whereas the nonlinear benchmarks retain a wider range to accommodate dynamically generated harmonics\.

Table 4:Spectral filtering and time\-step parameters used by the embedded physics component\. HerekNyq=N/2k\_\{\\mathrm\{Nyq\}\}=N/2denotes the maximum one\-dimensional Fourier index\.The available governing operator is advanced using a classical fourth\-order Runge–Kutta integrator\. At every RK stage, the intermediate state is passed through the frozen FFGS encoder, the required spatial derivatives are rendered analytically and filtered, and the governing right\-hand side is assembled\. For the nonlinear benchmarks, the assembled right\-hand side is additionally filtered as described above\. The FFGS normalization statistics are recomputed for each intermediate state\. One RK4 step is used per dataset interval, with no temporal substepping, so the integration step equals the frame spacingΔ​t\\Delta tlisted in Table[4](https://arxiv.org/html/2608.24049#A2.T4)\.

After the four RK stages are combined, the updated state is filtered once more,

ΦΔ​tFFGS​\(𝐮⁡\(t\),𝝀\)=Πκ​\[𝐮⁡\(t\)\+Δ​t6​\(K1\+2​K2\+2​K3\+K4\)\]\.\\Phi\_\{\\Delta t\}^\{\\mathrm\{FFGS\}\}\\\!\\left\(\\mathbf\{u\}\(t\);\\boldsymbol\{\\lambda\}\\right\)=\\Pi\_\{\\kappa\}\\left\[\\mathbf\{u\}\(t\)\+\\frac\{\\Delta t\}\{6\}\\left\(K\_\{1\}\+2K\_\{2\}\+2K\_\{3\}\+K\_\{4\}\\right\)\\right\]\.\(31\)whereK1,…,K4K\_\{1\},\\ldots,K\_\{4\}denote the standard RK4 stage evaluations of the FFGS\-based governing right\-hand side\. This final filtering step restricts the output of the embedded\-physics component to the retained spectral band and prevents unresolved high\-wavenumber content from accumulating within the embedded\-physics component over repeated updates\. All spectral filters described in this section are confined to the embedded\-physics component; the neural\-operator component acts directly on the current grid state without spectral filter\.

## Appendix CNeural operator implementation and training protocol

For the five fixed\-coefficient benchmarks, the neural\-operator component is implemented using theClassicFNOarchitecture provided by thepdequinoxlibrary\. The network contains four Fourier blocks, ten retained Fourier modes per spatial dimension, six hidden channels, and GELU activations\. This configuration contains 57,787 trainable parameters for the two\-dimensional scalar benchmarks, 57,800 parameters for the two\-component Burgers benchmark, and 1,152,187 parameters for the three\-dimensional benchmark\.

In both the fixed\-coefficient and coefficient\-conditioned implementations, the final output projection layer is initialized to zero\. Consequently, before Stage 2 optimization, the neural\-operator contribution vanishes and the composite surrogate reduces exactly to the frozen embedded\-physics operator\.

### C\.1Coefficient conditioning

The coefficient\-varying and coefficient\-unknown experiments use a separate Flax implementation of a coefficient\-conditioned FNO\. Its retained Fourier modes, hidden width, number of Fourier blocks, and overall parameter budget are matched to the unconditioned FNO configuration, while Feature\-wise Linear Modulation \(FiLM\) is introduced within every Fourier block\.

The coefficient vector𝝀\\boldsymbol\{\\lambda\}is standardized using statistics computed from the training coefficients\. For each Fourier blockℓ\\ell, an independent affine layer, with no hidden layer or activation, maps the standardized coefficient vector to a scale vector𝐬\(ℓ\)\\mathbf\{s\}^\{\(\\ell\)\}and bias vector𝐛\(ℓ\)\\mathbf\{b\}^\{\(\\ell\)\}, with one entry per hidden channel\. The block features are modulated as

𝐳\(ℓ\)←𝐳\(ℓ\)⊙\(𝟏\+𝐬\(ℓ\)\)\+𝐛\(ℓ\),\\mathbf\{z\}^\{\(\\ell\)\}\\leftarrow\\mathbf\{z\}^\{\(\\ell\)\}\\odot\\left\(\\mathbf\{1\}\+\\mathbf\{s\}^\{\(\\ell\)\}\\right\)\+\\mathbf\{b\}^\{\(\\ell\)\},\(32\)where𝐳\(ℓ\)\\mathbf\{z\}^\{\(\\ell\)\}denotes the hidden feature field of theℓ\\ell\-th Fourier block\. The factor𝟏\+𝐬\(ℓ\)\\mathbf\{1\}\+\\mathbf\{s\}^\{\(\\ell\)\}parameterizes the multiplicative modulation around the identity map\.

### C\.2Training protocol

Training follows the two\-stage procedure described in Section[2\.3](https://arxiv.org/html/2608.24049#S2.SS3)\. During Stage 1, the FFGS encoder is optimized using the frame\-reconstruction objective in Eq\. \([10](https://arxiv.org/html/2608.24049#S2.E10)\)\. After the prescribed number of epochs, the final encoder checkpoint is frozen and used for the corresponding Stage 2 training and evaluation\. During Stage 2, only the neural\-operator parameters are optimized using the five\-step rollout loss in Eq\. \([12](https://arxiv.org/html/2608.24049#S2.E12)\)\. The complete embedded\-physics component, including the frozen FFGS evaluation and RK4 integration, is enclosed by the stop\-gradient operation\.

The principal optimization settings are summarized in Table[5](https://arxiv.org/html/2608.24049#A3.T5)\. Stage 1 uses cosine learning\-rate decay from1\.5×10−31\.5\\times 10^\{\-3\}to10−610^\{\-6\}\. Stage 2 uses a 2,000\-step linear warmup to a peak learning rate of10−310^\{\-3\}, followed by cosine decay to zero over 10,000 optimization steps\.

Table 5:Training hyperparameters for the two\-stage optimization procedure\. No validation\-based checkpoint selection or early stopping is used\. For the three\-dimensional benchmark, Stage 1 uses 40 epochs and batch size 8, while Stage 2 uses a global batch size of 12\. The coefficient\-conditioned experiments use 100 Stage 1 epochs and batch size 64\.

## Appendix DCoefficient identification and coefficient\-conditioned implementation

This appendix describes the implementation of the coefficient\-unknown advection–diffusion experiments presented in Section[4\.5](https://arxiv.org/html/2608.24049#S4.SS5), including coefficient sampling, trajectory\-wise identification, coefficient conditioning, and the training and evaluation protocol\.

### D\.1Coefficient\-unknown dataset

The benchmark is governed by

∂tu=−v⁡\(∂xu\+∂yu\)\+D​∇2u,\\partial\_\{t\}u=\-v\\left\(\\partial\_\{x\}u\+\\partial\_\{y\}u\\right\)\+D\\nabla^\{2\}u,\(33\)where the advection speedvvand diffusivityDDare unknown\. The velocity is isotropic, with𝐯=v​𝟏\\mathbf\{v\}=v\\mathbf\{1\}, so the trajectory\-specific coefficient vector is

𝝀=\(v,D\)∈ℝ2\.\\boldsymbol\{\\lambda\}=\(v,D\)\\in\\mathbb\{R\}^\{2\}\.For each trajectory,log⁡v\\log vandlog⁡D\\log Dare sampled independently and uniformly over the logarithms of the intervals

v∈\[0\.0125,0\.025\],D∈\[10−6,4×10−6\],v\\in\[0\.0125,\\,0\.025\],\\qquad D\\in\[10^\{\-6\},\\,4\\times 10^\{\-6\}\],respectively\. The sampled coefficients remain constant within one trajectory and vary across trajectories\. Training and test trajectories are sampled independently from the same distributions\. The training trajectories contain 50 evolution steps, while evaluation uses 200\-step autoregressive rollouts\.

### D\.2Coefficient identification

Coefficient identification is performed independently for every trajectory before Stage 2 rollout training or test\-time rollout\. Identification uses seven consecutive observed frames centered on one calibration time\. The temporal derivative at the center frame is approximated using the sixth\-order centered finite\-difference scheme,

∂tu^k=uk\+3−9​uk\+2\+45​uk\+1−45​uk−1\+9​uk−2−uk−360​Δ​t\.\\widehat\{\\partial\_\{t\}u\}\_\{k\}=\\frac\{u\_\{k\+3\}\-9u\_\{k\+2\}\+45u\_\{k\+1\}\-45u\_\{k\-1\}\+9u\_\{k\-2\}\-u\_\{k\-3\}\}\{60\\Delta t\}\.\(34\)
The frozen FFGS map is applied to the center frame to obtain the required closed\-form spatial derivatives\. For coefficient identification, these derivatives are used without the rollout low\-pass filter\. We introduce the auxiliary regression coefficient vector

𝜷=\(vx,vy,D\),\\boldsymbol\{\\beta\}=\(v\_\{x\},v\_\{y\},D\),\(35\)where the directional advection coefficientsvxv\_\{x\}andvyv\_\{y\}are estimated independently in the least\-squares solve and subsequently combined to recover the isotropic coefficientvv\.

The regression system is then written as

𝐀𝜷≈𝐛,𝐀=\[−∂xu^,−∂yu^,Δu^\],𝐛=∂tu^\.\\mathbf\{A\}\\boldsymbol\{\\beta\}\\approx\\mathbf\{b\},\\qquad\\mathbf\{A\}=\\left\[\-\\partial\_\{x\}\\hat\{u\},\\,\-\\partial\_\{y\}\\hat\{u\},\\,\\Delta\\hat\{u\}\\right\],\\qquad\\mathbf\{b\}=\\widehat\{\\partial\_\{t\}u\}\.\(36\)where the spatial grid points are flattened into the rows of𝐀\\mathbf\{A\}and𝐛\\mathbf\{b\}\.

Before solving, each column of𝐀\\mathbf\{A\}is normalized by its Euclidean norm,

𝐀~=𝐀𝐒−1,𝐒=diag⁡\(s1,s2,s3\),\\widetilde\{\\mathbf\{A\}\}=\\mathbf\{A\}\\mathbf\{S\}^\{\-1\},\\qquad\\mathbf\{S\}=\\operatorname\{diag\}\(s\_\{1\},s\_\{2\},s\_\{3\}\),\(37\)wheresj=∥𝐀:,j∥2s\_\{j\}=\\\|\\mathbf\{A\}\_\{:,j\}\\\|\_\{2\}denotes the Euclidean norm of thejj\-th column of𝐀\\mathbf\{A\}\. The normalized system is then solved by ordinary least squares,

𝜷^=arg​min𝜷⁡‖𝐀~​𝜷−𝐛‖22\.\\hat\{\\boldsymbol\{\\beta\}\}=\\operatorname\*\{arg\\,min\}\_\{\\boldsymbol\{\\beta\}\}\\left\\\|\\widetilde\{\\mathbf\{A\}\}\\boldsymbol\{\\beta\}\-\\mathbf\{b\}\\right\\\|\_\{2\}^\{2\}\.\(38\)
The solution is transformed back to the original physical scaling\. Since the benchmark assumes an isotropic advection velocity, the two directional advection coefficients are averaged,

𝜷^phys=𝐒−1​𝜷^,v^=β^phys,1\+β^phys,22,D^=β^phys,3\.\\hat\{\\boldsymbol\{\\beta\}\}\_\{\\rm phys\}=\\mathbf\{S\}^\{\-1\}\\hat\{\\boldsymbol\{\\beta\}\},\\qquad\\hat\{v\}=\\frac\{\\hat\{\\beta\}\_\{\{\\rm phys\},1\}\+\\hat\{\\beta\}\_\{\{\\rm phys\},2\}\}\{2\},\\qquad\\hat\{D\}=\\hat\{\\beta\}\_\{\{\\rm phys\},3\}\.\(39\)
The final coefficient estimate is therefore𝝀^=\(v^,D^\)\\hat\{\\boldsymbol\{\\lambda\}\}=\(\\hat\{v\},\\hat\{D\}\)\. No explicit regularization is used; numerical conditioning is provided by column normalization before solving the least\-squares system\.

### D\.3Coefficient\-conditioned implementation

The same identified coefficient vector is used by both branches of the composite surrogate\. The embedded\-physics component substitutes\(v^,D^\)\(\\hat\{v\},\\hat\{D\}\)directly into the advection–diffusion operator, while the neural\-operator component receives the same coefficient estimate through the FiLM\-conditioned neural operator described in[C](https://arxiv.org/html/2608.24049#A3)\.

The FiLM normalization statistics are computed from the identified coefficients of the training trajectories\. During Stage 2 training and evaluation, both the embedded\-physics component and the neural operator use only the identified coefficients; the ground\-truth coefficient values are never supplied to either branch\.

Coefficient identification uses the unfiltered FFGS derivatives, whereas rollout prediction employs the spectrally filtered embedded\-physics component described in[B](https://arxiv.org/html/2608.24049#A2)\. The coefficient\-unknown experiments use a sharp spectral cutoff ofκ=7\\kappa=7\.

### D\.4Training and evaluation protocol

For each training trajectory, coefficient identification is performed once before Stage 2 optimization\. The recovered coefficients remain fixed throughout rollout training and are never updated online\.

For each held\-out test trajectory, the same seven\-frame calibration procedure is performed before autoregressive rollout\. The recovered coefficients remain fixed throughout prediction, and no reference states beyond the calibration window are used during coefficient estimation\.

Relative to the fixed\-coefficient experiments, the coefficient\-unknown experiments differ in three respects\. First, the FFGS encoder is trained for 100 epochs with batch size 64\. Second, the neural\-operator component is replaced by the FiLM\-conditioned FNO described in[C](https://arxiv.org/html/2608.24049#A3)\. Third, the embedded\-physics component uses a sharp spectral cutoff ofκ=7\\kappa=7\. Otherwise, the rollout horizon \(50 training steps and 200 evaluation steps\), rollout loss, optimizer settings, and checkpoint policy remain unchanged\. All reported rollout and spectral metrics are computed using the identified coefficients rather than the ground\-truth values\.

## Appendix ESupplementary experimental results

This appendix provides additional experimental results that complement the quantitative comparisons presented in the main text\. Unless otherwise stated, all experimental settings are identical to those described in Section[3](https://arxiv.org/html/2608.24049#S3)\.

### E\.1Additional long\-horizon rollout results

Figure 9:Long\-horizon test\-rollout accuracy on the remaining fixed\-coefficient benchmarks\. Relativeℓ2\\ell\_\{2\}error is averaged over the 30 held\-out trajectories and plotted as a function of autoregressive lead time on a logarithmic scale\. \(a\) Two\-dimensional advection, \(b\) two\-dimensional advection–diffusion, \(c\) two\-dimensional advection–Allen–Cahn, and \(d\) three\-dimensional advection–diffusion\. All models are trained on 50\-step trajectories and evaluated over 200\-step autoregressive rollouts\.Figure[9](https://arxiv.org/html/2608.24049#A5.F9)extends the long\-horizon comparison in Fig\.[2](https://arxiv.org/html/2608.24049#S4.F2)to the remaining fixed\-coefficient benchmarks\. Across all four cases, the proposed surrogate maintains the lowest relativeℓ2\\ell\_\{2\}error over the full 200\-step rollout\. The advantage persists across linear transport, transport–diffusion, reaction\-driven interface dynamics, and three\-dimensional evolution, indicating that the improved rollout stability is not specific to the Burgers benchmark\.

For the two\-dimensional advection and advection–diffusion cases, the error gap between the proposed surrogate and the purely data\-driven baselines increases progressively with lead time\. A similar trend is observed for three\-dimensional advection–diffusion, where the proposed surrogate exhibits the slowest long\-horizon error growth\. On the advection–Allen–Cahn benchmark, the proposed surrogate and ResNet remain substantially more accurate than FNO and U\-Net, with the proposed surrogate achieving the lowest error over most of the rollout horizon\.These results are consistent with the aggregate trajectory\-averaged metrics reported in Table[2](https://arxiv.org/html/2608.24049#S4.T2)\.

![Refer to caption](https://arxiv.org/html/2608.24049v1/snapshot_AC.png)Figure 10:Representative long\-horizon rollout of the two\-dimensional advection–Allen–Cahn benchmark\. Columns show the initial condition and autoregressive predictions att=33t=33,6666,100100,133133,166166, and199199\. Rows correspond to the reference solution, ResNet, FNO, U\-Net, and the proposed surrogate\. A common color scale is used for the scalar order parameteruuacross all methods and lead times\.Figure[10](https://arxiv.org/html/2608.24049#A5.F10)shows the evolution of the order parameter in the advection–Allen–Cahn benchmark\. All methods reproduce the initial coarsening of the phase\-field structures, but their long\-horizon behavior differs substantially\. FNO develops pronounced oscillatory artifacts after the intermediate rollout times and eventually departs qualitatively from the reference dynamics\. ResNet and U\-Net retain the dominant interfaces but exhibit increasing discrepancies in their positions and shapes\. The proposed surrogate most closely preserves the reference interface geometry and large\-scale phase organization throughout the rollout, consistent with its lower trajectory\-averaged error in Fig\.[9](https://arxiv.org/html/2608.24049#A5.F9)\(c\)\.

![Refer to caption](https://arxiv.org/html/2608.24049v1/snapshot_AD.png)Figure 11:Representative long\-horizon rollout of the two\-dimensional advection–diffusion benchmark\. Columns show the initial condition and autoregressive predictions att=33t=33,6666,100100,133133,166166, and199199\. Rows correspond to the reference solution, ResNet, FNO, U\-Net, and the proposed surrogate\. The scalar fielduuis displayed using a common color scale across all methods and lead times\.Figure[11](https://arxiv.org/html/2608.24049#A5.F11)compares the transported and diffused scalar structures in the two\-dimensional advection–diffusion benchmark\. The methods remain visually similar at early lead times, but spatial and amplitude errors become increasingly apparent during the later rollout\. ResNet exhibits the strongest long\-time deformation, while U\-Net shows noticeable changes in the morphology and intensity of individual structures\. FNO preserves much of the reference organization but accumulates visible local discrepancies\. The proposed surrogate maintains the closest agreement with the reference field over the full rollout, including the locations, signs, and characteristic scales of the dominant scalar structures\.

![Refer to caption](https://arxiv.org/html/2608.24049v1/snapshot_burgers.png)Figure 12:Representative long\-horizon rollout of the two\-dimensional Burgers benchmark\. Columns show the initial condition and autoregressive predictions att=40t=40,8080,120120,160160, and200200\. Rows correspond to the reference solution, ResNet, FNO, U\-Net, and the proposed surrogate\. The displayed quantity is the vorticity normalized by the instantaneous reference standard deviation,ω/σω​\(t\)\\omega/\\sigma\_\{\\omega\}\(t\); the corresponding values ofσω​\(t\)\\sigma\_\{\\omega\}\(t\)are reported below each column\.Figure[12](https://arxiv.org/html/2608.24049#A5.F12)provides a field\-level view of the Burgers rollout through the normalized vorticity\. As the flow evolves, the initial small\-scale structures merge into progressively smoother and more elongated vortical regions\. All surrogate models reproduce the dominant early\-time evolution, but the purely data\-driven baselines accumulate increasing phase, amplitude, and morphological errors at longer lead times\. ResNet exhibits substantial displacement and deformation of the principal structures, while U\-Net increasingly overpredicts their spatial extent\. FNO remains closer to the reference but still develops visible long\-time discrepancies\. The proposed surrogate most faithfully tracks the location, orientation, and width of the dominant vorticity structures throught=200t=200, in agreement with the rollout\-error comparison in Fig\.[2](https://arxiv.org/html/2608.24049#S4.F2)\.

Taken together, the snapshot comparisons show that the lower aggregate rollout errors of the proposed surrogate correspond to improved preservation of physically relevant spatial structures rather than only smaller pointwise discrepancies\. The advantage appears across transport–diffusion, reaction\-interface, and nonlinear self\-advection dynamics, and becomes most visible at lead times well beyond the training horizon\.

### E\.2Additional spectral comparisons

Figure 13:Time\-averaged radial energy spectra for the remaining fixed\-coefficient benchmarks\. Panels show \(a\) advection–Allen–Cahn \(2D\), \(b\) advection \(2D\), \(c\) advection–diffusion \(2D\), and \(d\) advection–diffusion \(3D\)\. Spectra are averaged over the 30 held\-out trajectories and 200\-step autoregressive rollouts\. The zero\-wavenumber mode is retained only for the advection–Allen–Cahn benchmark, following the evaluation protocol in Section[3\.4](https://arxiv.org/html/2608.24049#S3.SS4)\.Figure[13](https://arxiv.org/html/2608.24049#A5.F13)extends the spectral comparison of Fig\.[4](https://arxiv.org/html/2608.24049#S4.F4)to the remaining fixed\-coefficient benchmarks\. Consistent with the Burgers results, the proposed surrogate provides the closest agreement with the reference time\-averaged spectrum across all four cases\. In particular, it maintains substantially lower high\-wavenumber spectral distortion than the purely data\-driven baselines throughout the 200\-step autoregressive rollouts\. These qualitative observations are consistent with the spectral errors summarized in Table[2](https://arxiv.org/html/2608.24049#S4.T2)\.

### E\.3Model\-capacity control on the Burgers benchmark

Figure 14:Model\-capacity control on the Burgers benchmark\. \(a\) Long\-horizon test\-rollout relativeℓ2\\ell\_\{2\}error for the proposed surrogate and purely data\-driven baselines with matched and tenfold trainable dynamics\-map budgets\. Solid lines denote the matched\-budget models, and dashed lines denote the enlarged baselines\. \(b\) Relative error of the time\-averaged radial spectrum,‖E¯−E¯GT‖2/‖E¯GT‖2\\\|\\overline\{E\}\-\\overline\{E\}\_\{\\mathrm\{GT\}\}\\\|\_\{2\}/\\\|\\overline\{E\}\_\{\\mathrm\{GT\}\}\\\|\_\{2\}, evaluated overκ≥1\\kappa\\geq 1following Section[3\.4](https://arxiv.org/html/2608.24049#S3.SS4)\.Figure[14](https://arxiv.org/html/2608.24049#A5.F14)examines whether the performance gain of the proposed surrogate can be explained solely by increased model capacity\. The trainable dynamics\-map budgets of FNO, U\-Net, and ResNet are increased by a factor of ten, whereas the proposed surrogate uses the same configuration as in the main experiments\.

As shown in Fig\.[14](https://arxiv.org/html/2608.24049#A5.F14)\(a\), increasing the capacity of the purely data\-driven baselines consistently improves long\-horizon rollout accuracy\. The enlarged FNO achieves the best rollout performance among the purely data\-driven models and slightly outperforms the proposed surrogate over part of the rollout horizon\. Nevertheless, the proposed surrogate remains more accurate than both the matched\-budget baselines and the enlarged U\-Net and ResNet models\.

The time\-averaged spectral error exhibits a different ordering, as shown in Fig\.[14](https://arxiv.org/html/2608.24049#A5.F14)\(b\)\. Despite using the same trainable dynamics\-map budget as the matched FNO baseline, the proposed surrogate achieves the lowest spectral error among all compared models\. Increasing the capacity of the purely data\-driven baselines reduces, but does not eliminate, the spectral discrepancy\. This observation suggests that the improved spectral fidelity cannot be explained solely by increased neural\-model capacity\.

## Appendix FRepresentation diagnostics

To assess the suitability of the frozen FFGS representation for physics evaluation, we evaluate field\- and derivative\-level reconstruction errors on the periodic test problems\. Reference spatial derivatives are computed by pseudo\-spectral differentiation\. For the FFGS reconstructionu^\\hat\{u\}, we define

eu=‖u^−u‖2‖u‖2,e∇u=‖∇u^−∇u‖2‖∇u‖2,eΔ​u=‖Δ​u^−Δ​u‖2‖Δ​u‖2\.e\_\{u\}=\\frac\{\\\|\\hat\{u\}\-u\\\|\_\{2\}\}\{\\\|u\\\|\_\{2\}\},\\qquad e\_\{\\nabla u\}=\\frac\{\\\|\\nabla\\hat\{u\}\-\\nabla u\\\|\_\{2\}\}\{\\\|\\nabla u\\\|\_\{2\}\},\\qquad e\_\{\\Delta u\}=\\frac\{\\\|\\Delta\\hat\{u\}\-\\Delta u\\\|\_\{2\}\}\{\\\|\\Delta u\\\|\_\{2\}\}\.\(40\)
Table[6](https://arxiv.org/html/2608.24049#A6.T6)reports the resulting errors for each benchmark\.

Table 6:RelativeL2L\_\{2\}reconstruction errors of the frozen FFGS representation for the field, gradient, and Laplacian\. Reference spatial derivatives are obtained by pseudo\-spectral differentiation\.The frozen FFGS representation maintains low field\-level reconstruction errors across all benchmarks, while the errors increase under spatial differentiation and are largest for the Laplacian\. This increased sensitivity of higher\-order derivatives to reconstruction error motivates the spectral filtering in Eq\. \([9](https://arxiv.org/html/2608.24049#S2.E9)\) before the differential quantities are used to evaluate the available governing terms\.