Geometry-aware Incremental Neural Operator for Long-Horizon PDE prediction

# Geometry-aware Incremental Neural Operator for Long-Horizon PDE prediction
Source: [https://arxiv.org/html/2608.11237](https://arxiv.org/html/2608.11237)
,Shuxu ChenElectronics and Information Convergence Engineering, KHUYongin\-siKorea,Haifan MengSchool of Computer Science and Engineering, UESTCChengduChina,Yi LuDepartment of Mathematical Sciences, UOLLiverpoolEngland,Zhihan LyuSchool of Computer Science and Technology, XDUShanxiChina,Fan MoSchool of Mechanical and Electrical Engineering, UESTCChengduChina,Wei DongCollege of Computer and Information Engineering, XAUATXi’anChina,Yang YangSchool of Computer Science and Engineering, UESTCChengduChinaandChaoning ZhangSchool of Computer Science and Engineering, UESTCChengduChina

\(5 June 2009\)

###### Abstract\.

Neural operators have shown strong potential for learning solution operators of partial differential equations \(PDEs\)\. However, long\-horizon autoregressive prediction remains challenging: local errors accumulate as spectral inconsistency, phase misalignment, or mean drift\. Existing methods mainly improve state representations and operator backbones, while leaving the repeatedly applied latent transition increment weakly structured, allowing spectral errors and unstable channel couplings to accumulate during rollout\. To address these issues, we propose a geometry\-aware incremental neural operator \(GeoIncNO\) for stable long\-horizon PDE prediction\. GeoIncNO predicts latent increments for residual advancement and uses lightweight low\-rank projectors to regulate channel coupling within active frequency bands derived from the increment spectral energy distribution\. To reduce physical\-space reconstruction errors, GeoIncNO further introduces a mean–fluctuation decoupled reconstruction mechanism, where stable mean structures and dynamic fluctuations are fused separately, and phase correction is applied only to the zero\-mean fluctuation component\. Extensive experiments on six PDE benchmarks, covering 1D, 2D, and 3D dynamical systems, show that GeoIncNO achieves consistently strong prediction accuracy, improved rollout stability, and better spectral fidelity compared with competitive neural\-operator baselines\.

††copyright:acmlicensed††journalyear:2018††doi:XXXXXXX\.XXXXXXX††conference:Make sure to enter the correct conference title from your rights confirmation email; June 03–05, 2018; Woodstock, NY††ccs:Computing methodologies Neural networks## 1\.Introduction

Partial differential equations \(PDEs\) describe the spatiotemporal evolution of many complex physical systems, such as fluid flows, heat transfer and multiscale materials\(Yip et al\.,[2022](https://arxiv.org/html/2608.11237#bib.bib44); Brunton and Kutz,[2024](https://arxiv.org/html/2608.11237#bib.bib6); Aarts and Van Der Veer,[2001](https://arxiv.org/html/2608.11237#bib.bib2)\)\. Although traditional numerical solvers are accurate, they often require fine\-grid discretization and iterative computation, which becomes costly in high\-dimensional and long\-horizon settings\(Karniadakis et al\.,[2021](https://arxiv.org/html/2608.11237#bib.bib17)\)\. Neural operators offer a data\-driven paradigm for PDE modeling by directly learning mappings between function spaces\(Azizzadenesheli et al\.,[2024](https://arxiv.org/html/2608.11237#bib.bib3); Liu et al\.,[2025](https://arxiv.org/html/2608.11237#bib.bib25); Dai et al\.,[2026](https://arxiv.org/html/2608.11237#bib.bib9); Ren et al\.,[2026](https://arxiv.org/html/2608.11237#bib.bib33)\)\. This enables fast prediction of future physical fields from initial conditions, making neural operators an important tool for scientific machine learning\(Li et al\.,[2021a](https://arxiv.org/html/2608.11237#bib.bib22)\)\. Despite the significant progress of existing neural operators in one\-step and short\-term prediction, long\-horizon autoregressive prediction remains challenging due to error accumulation\(Brandstetter et al\.,[2022](https://arxiv.org/html/2608.11237#bib.bib4); McCabe et al\.,[2023](https://arxiv.org/html/2608.11237#bib.bib28); Worrall et al\.,[2024](https://arxiv.org/html/2608.11237#bib.bib39)\)\. Since the prediction at each step is repeatedly fed back as the input for subsequent steps, local errors can be amplified during rollout, leading to degraded spectral consistency and unstable temporal evolution\(Wu et al\.,[2025](https://arxiv.org/html/2608.11237#bib.bib40)\)\.

Existing neural PDE predictors can be broadly grouped according to where temporal evolution is modeled\. One class learns the solution evolution in the physical grid space or in lifted operator features, without explicitly formulating a compact latent\-state transition process\. These methods typically improve prediction by introducing Fourier parameterization, multiscale structures, local convolutions, or stronger spatiotemporal mixing modules\(Rahman et al\.,[2023](https://arxiv.org/html/2608.11237#bib.bib30); Hu et al\.,[2024](https://arxiv.org/html/2608.11237#bib.bib16)\)\. Another class first maps physical fields into a latent space and then learns the evolution of latent states, thereby reducing the modeling complexity in the original physical space\(Wang and Wang,[2024](https://arxiv.org/html/2608.11237#bib.bib37); Tiwari et al\.,[2025](https://arxiv.org/html/2608.11237#bib.bib35); Chen and Wu,[2025](https://arxiv.org/html/2608.11237#bib.bib8)\)\. Although latent\-space neural operators provide a compact representation space for PDE dynamics, existing formulations usually model the evolution by mapping one latent state to the next\. The latent change between consecutive states remains entangled with the state representation, so its frequency\-domain structure and channel\-wise coupling are difficult to characterize or regulate\. This entanglement is particularly problematic under autoregressive rollout, where the same transition update is applied repeatedly and small per\-step errors accumulate over long horizons\. Therefore, stable long\-horizon prediction requires not only an effective latent\-state representation, but also a structured increment formulation that directly regularizes the repeated transition process\. Motivated by this limitation, we examine latent\-state propagation from the perspective of the transition increment\. As shown in Section[3](https://arxiv.org/html/2608.11237#S3), latent states and latent increments exhibit different dominant spatial responses\. Specifically, the latent stateztz\_\{t\}retains broad spatial organization, whereas the transition incrementδ​zt\\delta z\_\{t\}exhibits a distinct response associated with the change between consecutive latent states \(Figure[1](https://arxiv.org/html/2608.11237#S3.F1)\)\. The covariance spectrum ofδ​zt\\delta z\_\{t\}also decays faster than that ofztz\_\{t\}, with more of its variation captured by the leading channel directions \(Figure[2](https://arxiv.org/html/2608.11237#S3.F2)\)\. In addition, its spectral energy is distributed non\-uniformly across spatial frequencies \(Figure[3](https://arxiv.org/html/2608.11237#S3.F3)\)\. Together, these observations motivate modeling the latent increment explicitly in the frequency–channel domain\. We refer to this modeling perspective as increment geometry\. GeoIncNO accordingly applies structured modulation to the increment before residual state advancement\.

Building on this increment\-centered perspective, we propose a Geometry\-aware Incremental Neural Operator \(GeoIncNO\) for long\-horizon PDE prediction\. GeoIncNO first encodes the input physical field into a latent state and uses a latent\-space backbone to predict a raw transition increment\. Instead of directly applying this increment for state update, GeoIncNO introduces Active\-Band Projection \(ABP\), which constructs active frequency bands from the spectral structure of the increment and uses lightweight low\-rank geometric projectors to modulate band\-wise channel coupling\. The resulting geometrically shaped increment is then used for residual latent\-state advancement, allowing the frequency\-band structure to directly participate in forward propagation and mitigate unstable spectral components and redundant channel coupling that may otherwise accumulate during rollout\. GeoIncNO further introduces Mean–Fluctuation Decoupled Reconstruction \(MFDR\) to address physical\-space reconstruction errors\. Since the candidate physical predictions may contain slow mean shifts, dynamic fluctuations, and phase misalignment, directly fusing complete predicted fields can entangle these errors\. MFDR therefore uses two complementary prediction branches, decomposes their outputs into temporal mean and zero\-mean fluctuation fields, and fuses the two components separately\. Phase correction is applied only to the zero\-mean fluctuation field, reducing interference with the mean structure and alleviating mean drift during long\-horizon prediction\.

The main contributions of this work are summarized as follows\.\(i\)We introduce an increment\-centered perspective for neural PDE operators, where the latent transition increment is explicitly modeled as the object governing autoregressive state propagation\.\(ii\)We propose an active\-band increment projection mechanism that constructs frequency bands from the spectral energy of latent increments and applies lightweight low\-rank projectors to shape band\-wise channel geometry before latent\-state advancement\.\(iii\)We design a mean–fluctuation decoupled reconstruction strategy that separately fuses stable mean structures and zero\-mean fluctuations, while restricting phase correction to the fluctuation field to reduce mean\-field drift\.\(iv\)We validate GeoIncNO on six PDE benchmarks spanning 1D, 2D, and 3D systems, demonstrating improved prediction accuracy, spectral fidelity, and rollout stability over competitive neural\-operator baselines\. Ablation studies further confirm the contribution of each component\.

## 2\.Related Work

Neural operators\(Kovachki et al\.,[2023](https://arxiv.org/html/2608.11237#bib.bib20)\)learn mappings between function spaces, unlike physics\-informed neural networks\(Raissi et al\.,[2019](https://arxiv.org/html/2608.11237#bib.bib31)\)that embed governing equations into the loss and retrain per instance; FNO\(Li et al\.,[2021a](https://arxiv.org/html/2608.11237#bib.bib22)\)and DeepONet\(Lu et al\.,[2021](https://arxiv.org/html/2608.11237#bib.bib26)\)are representative\. Subsequent work either strengthens the operator backbone, via adaptive Fourier mixing\(Guibas et al\.,[2022](https://arxiv.org/html/2608.11237#bib.bib10)\), spectral\-basis enrichment\(Gupta et al\.,[2021](https://arxiv.org/html/2608.11237#bib.bib11); Hu et al\.,[2025](https://arxiv.org/html/2608.11237#bib.bib15)\), attention\-based mixing\(Hao et al\.,[2023](https://arxiv.org/html/2608.11237#bib.bib13)\), and Koopman\-linearized operators\(Xiong et al\.,[2024](https://arxiv.org/html/2608.11237#bib.bib43)\), or compresses PDE dynamics into a latent space\(Lusch et al\.,[2018](https://arxiv.org/html/2608.11237#bib.bib27)\), as in LSM\(Wu et al\.,[2023](https://arxiv.org/html/2608.11237#bib.bib41)\), LNO\(Wang and Wang,[2024](https://arxiv.org/html/2608.11237#bib.bib37)\), and PI\-Latent\-NO\(Karumuri et al\.,[2026](https://arxiv.org/html/2608.11237#bib.bib18)\)\. Both routes model the state\-to\-state or latent\-state mapping, however, leaving the latent increment that drives temporal evolution entangled with the state representation rather than modeled as a structured object in its own right\. We provide the complete related work in Appendix[A](https://arxiv.org/html/2608.11237#A1)\.

## 3\.Motivation

This section presents the core empirical observations that motivate GeoIncNO’s latent\-increment modeling and active\-band construction\. We first show that latent increments differ from full latent states in spatial response and covariance geometry\. We then analyze their non\-uniform spectral energy distribution, which motivates energy\-adaptive frequency partitioning\. Additional analyses of frequency\-dependent band geometry and mean–fluctuation separation are provided in Appendix[B](https://arxiv.org/html/2608.11237#A2)\.

![Refer to caption](https://arxiv.org/html/2608.11237v1/Figure/Figure1.png)Figure 1\.First principal\-component projections ofztz\_\{t\}andδ​zt\\delta z\_\{t\}\. PCA is applied independently along the latent channel dimension, so the two panels compare dominant spatial patterns rather than cross\-panel magnitudes\.##### Latent Increment Dynamics

To examine the roles of latent states and latent increments, we perform a forward pass using a trained GeoIncNO model on a held\-out 2D Navier–Stokes \(NS\) test sample\. We visualize the latent stateztz\_\{t\}and the corresponding latent incrementδ​zt\\delta z\_\{t\}in Figure[1](https://arxiv.org/html/2608.11237#S3.F1)\. Since both are multi\-channel latent fields, we apply principal component analysis \(PCA\) along the latent channel dimension and visualize the spatial response of the first principal component\. The PCA projections are computed independently, so the two panels compare dominant spatial patterns rather than cross\-panel magnitudes\.

As shown in Figure[1](https://arxiv.org/html/2608.11237#S3.F1),ztz\_\{t\}exhibits broad spatial organization, indicating that it preserves the current field configuration in the latent space\. In contrast,δ​zt\\delta z\_\{t\}presents a distinct dominant response pattern\. This difference suggests that the latent state and latent increment encode different information rather than repeating the same spatial structure\. Explicitly modeling the latent increment is therefore more targeted than imposing structure only on the full latent state, which may entangle the current field configuration with the transition information required for subsequent propagation\.

![Refer to caption](https://arxiv.org/html/2608.11237v1/Figure/Figure2.png)Figure 2\.Covariance spectral analysis ofztz\_\{t\}andδ​zt\\delta z\_\{t\}, showing faster eigenvalue decay and higher principal\-energy concentration for latent increments\.![Refer to caption](https://arxiv.org/html/2608.11237v1/Figure/Figure3.png)Figure 3\.Radial spectral energy distribution ofδ​zt\\delta z\_\{t\}\.
##### Covariance Geometry of Latent Increments

We analyze the covariance geometry of latent increments by applying a trained GeoIncNO model to held\-out 2D NS test samples and collectingztz\_\{t\}andδ​zt\\delta z\_\{t\}during inference\. We compute their covariance matrices across the latent channels and compare the resulting spectra in Figure[2](https://arxiv.org/html/2608.11237#S3.F2)\. Specifically, we perform eigenvalue decomposition and visualize the normalized eigenvalue spectra and cumulative energy captured by the top\-kkprincipal directions\. As shown in Figure[2](https://arxiv.org/html/2608.11237#S3.F2), the covariance eigenvalues ofδ​zt\\delta z\_\{t\}decay more rapidly than those ofztz\_\{t\}, and its cumulative principal energy rises faster within the first few directions\. This indicates that the variation ofδ​zt\\delta z\_\{t\}is concentrated in fewer principal directions, implying a more compact covariance structure\. Therefore, its dominant transition directions can be captured by a small number of principal modes, supporting the increment\-centered geometric design of GeoIncNO\.

##### Non\-uniform Increment Spectra

Next, we analyze the spectral distribution of latent increments\. We use a trained FNO checkpoint on the 2D shallow\-water equation \(SWE\) benchmark to verify that the observed spectral non\-uniformity is not specific to our model\. For each test sample, two consecutive physical contexts and their coordinates\(x,y\)\(x,y\)are passed through the same FNO lifting layer and GELU activation to obtainztz\_\{t\}andzt\+1z\_\{t\+1\}\. The latent increment is then computed asδ​zt=zt\+1−zt\\delta z\_\{t\}=z\_\{t\+1\}\-z\_\{t\}\. We apply a 2D real FFT toδ​zt\\delta z\_\{t\}over the spatial dimensions and average the spectral power over samples, time steps, and latent channels\. For each Fourier mode, letkxk\_\{x\}andkyk\_\{y\}denote the frequency indices along the two spatial axes, and define its radial frequency asρ​\(kx,ky\)=kx2\+ky2\\rho\(k\_\{x\},k\_\{y\}\)=\\sqrt\{k\_\{x\}^\{2\}\+k\_\{y\}^\{2\}\}\. As illustrated in Figure[3](https://arxiv.org/html/2608.11237#S3.F3), the black curve shows the normalized radial spectral energy ofδ​zt\\delta z\_\{t\}, obtained by grouping Fourier modes according toρ​\(kx,ky\)\\rho\(k\_\{x\},k\_\{y\}\)\. The energy is non\-uniformly distributed and is concentrated in low\-frequency regions and a few effective frequency intervals, decaying rapidly as the radial frequency increases\. The gray dashed lines denote a fixed uniform frequency partition, which assigns several bands to low\-energy or nearly inactive regions\. In contrast, the blue solid lines show active\-band boundaries that adapt to the observed spectral support ofδ​zt\\delta z\_\{t\}\. This comparison motivates constructing active bands according to the increment spectral energy, so that geometric modeling focuses on frequency regions carrying meaningful transition dynamics\.

![Refer to caption](https://arxiv.org/html/2608.11237v1/Figure/Figure7.png)Figure 4\.Overview of GeoIncNO\. ABP shapes the spectral–channel geometry of the latent increment before residual state advancement, while MFDR performs dual\-branch mean–fluctuation reconstruction with fluctuation\-only phase correction\.

## 4\.Method

GeoIncNO is built on a guiding principle: it explicitly treats the latent transition increment as the object that governs autoregressive evolution\. As shown in Figure[4](https://arxiv.org/html/2608.11237#S3.F4), GeoIncNO follows a geometry\-aware incremental pipeline: i\) the input physical field is lifted into a latent state, from which a raw latent transition increment is predicted; ii\) active spectral bands are constructed from the increment, and low\-rank band\-wise projection is applied to shape its spectral\-channel geometry; iii\) the refined increment advances the latent state through residual propagation; iv\) two complementary physical predictions are reconstructed from the updated state and the refined increment; and v\) temporal mean and zero\-mean fluctuation components are decoupled and fused separately, after which phase correction is applied only to the fluctuation component\.

### 4\.1\.Problem Formulation

We consider long\-horizon prediction for time\-dependent PDE systems\. Letuinu\_\{\\mathrm\{in\}\}denote the observed physical fields over an input time window, whereuin∈ℝB×Tin×N1in×⋯×Ndsin×Cinu\_\{\\mathrm\{in\}\}\\in\\mathbb\{R\}^\{B\\times T\_\{\\mathrm\{in\}\}\\times N^\{\\mathrm\{in\}\}\_\{1\}\\times\\cdots\\times N^\{\\mathrm\{in\}\}\_\{d\_\{s\}\}\\times C\_\{\\mathrm\{in\}\}\}\. Here,BBis the batch size,TinT\_\{\\mathrm\{in\}\}is the number of input time steps,dsd\_\{s\}is the spatial dimension,N1in,…,NdsinN^\{\\mathrm\{in\}\}\_\{1\},\\ldots,N^\{\\mathrm\{in\}\}\_\{d\_\{s\}\}are the input spatial grid sizes, andCinC\_\{\\mathrm\{in\}\}is the number of input physical variables\. The goal is to predict the future physical fieldsu^∈ℝB×Tout×N1out×⋯×Ndsout×Cout\\hat\{u\}\\in\\mathbb\{R\}^\{B\\times T\_\{\\mathrm\{out\}\}\\times N^\{\\mathrm\{out\}\}\_\{1\}\\times\\cdots\\times N^\{\\mathrm\{out\}\}\_\{d\_\{s\}\}\\times C\_\{\\mathrm\{out\}\}\}, whereToutT\_\{\\mathrm\{out\}\}is the prediction horizon,N1out,…,NdsoutN^\{\\mathrm\{out\}\}\_\{1\},\\ldots,N^\{\\mathrm\{out\}\}\_\{d\_\{s\}\}are the output spatial grid sizes, andCoutC\_\{\\mathrm\{out\}\}is the number of output physical variables\. Letuoutu\_\{\\mathrm\{out\}\}denote the corresponding ground\-truth future fields with the same shape\. The task is to learn a neural operator𝒢θ:uin↦u^\\mathcal\{G\}\_\{\\theta\}:u\_\{\\mathrm\{in\}\}\\mapsto\\hat\{u\}that approximates the PDE solution operator over the target window\.

### 4\.2\.Latent State Encoding and Increment Prediction

#### 4\.2\.1\.Encoding

GeoIncNO first encodes the input physical field into a grid\-structured latent representation, where subsequent dynamics are modeled through latent increments\. Specifically, the input field is mapped aszin=Eϕ​\(uin\)z\_\{\\mathrm\{in\}\}=E\_\{\\phi\}\(u\_\{\\mathrm\{in\}\}\), wherezin∈ℝB×d×Tin×N1in×⋯×Ndsinz\_\{\\mathrm\{in\}\}\\in\\mathbb\{R\}^\{B\\times d\\times T\_\{\\mathrm\{in\}\}\\times N^\{\\mathrm\{in\}\}\_\{1\}\\times\\cdots\\times N^\{\\mathrm\{in\}\}\_\{d\_\{s\}\}\}, anddddenotes the latent channel width\. The encoder preserves the spatiotemporal grid by concatenating the physical variables with normalized temporal and spatial coordinates, followed by a latent projection and several dimension\-specific spatial convolutional layers\. The resultingzinz\_\{\\mathrm\{in\}\}serves as the latent representation from which GeoIncNO predicts the transition direction\.

#### 4\.2\.2\.Increment Prediction

For notational simplicity, we denotezinz\_\{\\mathrm\{in\}\}byztz\_\{t\}when describing a single latent transition, wherettindexes the transition itself rather than the intra\-window temporal coordinateτ\\tau\. Givenztz\_\{t\}, GeoIncNO extracts transition features asht=𝒯θ​\(zt\)h\_\{t\}=\\mathcal\{T\}\_\{\\theta\}\(z\_\{t\}\)\. The backbone𝒯θ\\mathcal\{T\}\_\{\\theta\}consists of mixed latent blocks that combine spectral\-domain modeling, local convolution, pointwise channel mixing, and a channel MLP, allowing the model to capture both global modes and local spatiotemporal patterns\. An increment prediction headGψG\_\{\\psi\}then predicts the raw latent increment asδ​zraw=Gψ​\(ht\)\\delta z\_\{\\mathrm\{raw\}\}=G\_\{\\psi\}\(h\_\{t\}\)\. The raw increment provides an unconstrained estimate of the latent transition direction and is subsequently refined through the active\-band mechanism before latent\-state advancement\.

### 4\.3\.Active\-Band Increment Projection

The predicted incrementδ​zraw\\delta z\_\{\\mathrm\{raw\}\}encodes the transition direction of the latent dynamics\. Compared with the full latent stateztz\_\{t\}, which also contains static context, boundary\-related information, and sample\-specific background structures, the increment more directly captures the dynamic change that drives state evolution\. GeoIncNO therefore applies geometric shaping to the increment rather than to the full state\. Since the freely predicted increment may contain unstructured spectral components or redundant channel couplings, this shaping is performed within active spectral bands before residual state advancement\.

#### 4\.3\.1\.Active\-Band Construction

To focus geometric shaping on dynamically relevant frequency regions, GeoIncNO constructs active bands according to the spectral energy distribution of the raw latent increment\. The construction is defined on a generaldsd\_\{s\}\-dimensional spatial domain and therefore applies to one\-, two\-, and three\-dimensional PDE systems\. Givenδ​zraw∈ℝB×d×Tin×N1in×⋯×Ndsin\\delta z\_\{\\mathrm\{raw\}\}\\in\\mathbb\{R\}^\{B\\times d\\times T\_\{\\mathrm\{in\}\}\\times N^\{\\mathrm\{in\}\}\_\{1\}\\times\\cdots\\times N^\{\\mathrm\{in\}\}\_\{d\_\{s\}\}\}, whereN1in,…,NdsinN^\{\\mathrm\{in\}\}\_\{1\},\\ldots,N^\{\\mathrm\{in\}\}\_\{d\_\{s\}\}denote the spatial resolutions of the latent grid, we compute its spectral representation asδ​z^raw=ℱ𝐱​\(δ​zraw\)\\widehat\{\\delta z\}\_\{\\mathrm\{raw\}\}=\\mathcal\{F\}\_\{\\mathbf\{x\}\}\(\\delta z\_\{\\mathrm\{raw\}\}\), where𝐱=\(x1,…,xds\)\\mathbf\{x\}=\(x\_\{1\},\\ldots,x\_\{d\_\{s\}\}\)\. For periodic domains,ℱ𝐱\\mathcal\{F\}\_\{\\mathbf\{x\}\}is implemented as adsd\_\{s\}\-dimensional fast Fourier transform\. For non\-periodic discretizations, it can be replaced with a boundary\-compatible orthogonal spectral transform\. Let𝐤=\(k1,…,kds\)∈ℤds\\mathbf\{k\}=\(k\_\{1\},\\ldots,k\_\{d\_\{s\}\}\)\\in\\mathbb\{Z\}^\{d\_\{s\}\}denote a spatial frequency index\. The average spectral energy of the raw latent increment at𝐤\\mathbf\{k\}is defined as

\(1\)E​\(𝐤\)\\displaystyle E\(\\mathbf\{k\}\)=1B​Tin​d​∑b,t,c\|δ​z^raw​\(b,c,t,𝐤\)\|2\.\\displaystyle=\\frac\{1\}\{BT\_\{\\mathrm\{in\}\}d\}\\sum\_\{b,t,c\}\\left\|\\widehat\{\\delta z\}\_\{\\mathrm\{raw\}\}\(b,c,t,\\mathbf\{k\}\)\\right\|^\{2\}\.The active\-band partition is constructed within the effective spectral supportΩeff⊆ℤds\\Omega\_\{\\mathrm\{eff\}\}\\subseteq\\mathbb\{Z\}^\{d\_\{s\}\}, preventing inactive frequency locations from affecting the band boundaries\. To obtain a dimension\-independent radial partition, we define

\(2\)ρ​\(𝐤\)\\displaystyle\\rho\(\\mathbf\{k\}\)=\(∑j=1ds\(\|kj\|kjmax\+ϵ\)2\)1/2,\\displaystyle=\\left\(\\sum\_\{j=1\}^\{d\_\{s\}\}\\left\(\\frac\{\|k\_\{j\}\|\}\{k\_\{j\}^\{\\max\}\+\\epsilon\}\\right\)^\{2\}\\right\)^\{1/2\},wherekjmaxk\_\{j\}^\{\\max\}denotes the maximum retained frequency index along thejj\-th spatial direction andϵ\\epsilonis a numerical stability constant\.

GeoIncNO determines the active\-band boundaries from the cumulative increment spectral energy\. The cumulative energy function overΩeff\\Omega\_\{\\mathrm\{eff\}\}is defined as

\(3\)C​\(r\)\\displaystyle C\(r\)=∑𝐤∈Ωeffρ​\(𝐤\)≤rE​\(𝐤\)∑𝐤∈ΩeffE​\(𝐤\)\+ϵ\.\\displaystyle=\\frac\{\\sum\_\{\\begin\{subarray\}\{c\}\\mathbf\{k\}\\in\\Omega\_\{\\mathrm\{eff\}\}\\\\ \\rho\(\\mathbf\{k\}\)\\leq r\\end\{subarray\}\}E\(\\mathbf\{k\}\)\}\{\\sum\_\{\\mathbf\{k\}\\in\\Omega\_\{\\mathrm\{eff\}\}\}E\(\\mathbf\{k\}\)\+\\epsilon\}\.GivenMMactive bands, the boundary of themm\-th band is determined by the energy quantile

\(4\)em\\displaystyle e\_\{m\}=inf\{r∣C​\(r\)≥mM\},m=1,…,M\.\\displaystyle=\\inf\\left\\\{r\\mid C\(r\)\\geq\\frac\{m\}\{M\}\\right\\\},\\qquad m=1,\\ldots,M\.The active band isℬm=\{𝐤∈Ωeff∣em−1≤ρ​\(𝐤\)<em\}\\mathcal\{B\}\_\{m\}=\\left\\\{\\mathbf\{k\}\\in\\Omega\_\{\\mathrm\{eff\}\}\\mid e\_\{m\-1\}\\leq\\rho\(\\mathbf\{k\}\)<e\_\{m\}\\right\\\}\. The last band includes the right endpoint so that all frequency locations inΩeff\\Omega\_\{\\mathrm\{eff\}\}are covered\. For 1D PDEs, the active bands are intervals on the frequency axis\. For 2D PDEs, they correspond to radial annuli in the frequency plane, and for 3D PDEs, they correspond to radial shells in the frequency volume\. Thus, ABP constructs active bands in a unifieddsd\_\{s\}\-dimensional frequency space rather than relying on a specifically two\-dimensional formulation\.

To maintain stable band semantics, GeoIncNO calibrates the active\-band boundaries during the early training stage\. LetEb​\(𝐤\)E\_\{b\}\(\\mathbf\{k\}\)denote the increment spectral energy estimated from thebb\-th mini\-batch\. The training\-set estimate isE¯​\(𝐤\)=Nbatch−1​∑b=1NbatchEb​\(𝐤\)\\bar\{E\}\(\\mathbf\{k\}\)=N\_\{\\mathrm\{batch\}\}^\{\-1\}\\sum\_\{b=1\}^\{N\_\{\\mathrm\{batch\}\}\}E\_\{b\}\(\\mathbf\{k\}\)\. The cumulative energy function is then computed by replacingE​\(𝐤\)E\(\\mathbf\{k\}\)withE¯​\(𝐤\)\\bar\{E\}\(\\mathbf\{k\}\), and the boundaries\{em\}m=1M\\\{e\_\{m\}\\\}\_\{m=1\}^\{M\}are obtained using the same energy\-quantile rule\. Once calibrated, the boundaries remain fixed during subsequent training and inference\. If the effective spectral energy is insufficient or too few frequency locations are available, GeoIncNO falls back to a uniform radial partition\.

#### 4\.3\.2\.Low\-Rank Band Projection

After constructing the active bands, GeoIncNO applies a lightweight low\-rank projection within each band to reshape the channel coupling of the spectral increment\. Here, projection refers to residual low\-rank geometric shaping rather than a strict orthogonal projection\. For a frequency location𝐤∈ℬm\\mathbf\{k\}\\in\\mathcal\{B\}\_\{m\}, letv=a\+i​b∈ℂdv=a\+ib\\in\\mathbb\{C\}^\{d\}, wherea,b∈ℝda,b\\in\\mathbb\{R\}^\{d\}, denote the complex channel vector of the spectral increment\.

For each band, we introduce a learnable low\-rank channel basisUm∈ℝd×rU\_\{m\}\\in\\mathbb\{R\}^\{d\\times r\}, wherer≪dr\\ll d, and define

\(5\)Pmℝ​\(x\)\\displaystyle P\_\{m\}^\{\\mathbb\{R\}\}\(x\)=x\+αm​Um​Um⊤​x,\\displaystyle=x\+\\alpha\_\{m\}U\_\{m\}U\_\{m\}^\{\\top\}x,whereαm\\alpha\_\{m\}is a learnable projection strength andx∈ℝdx\\in\\mathbb\{R\}^\{d\}\. The low\-rank basis restricts band\-wise channel modulation to a few learnable directions, avoiding unconstrained full\-channel mixing\. The same real\-valued projector is applied to the real and imaginary parts asPm​\(v\)=Pmℝ​\(a\)\+i​Pmℝ​\(b\)P\_\{m\}\(v\)=P\_\{m\}^\{\\mathbb\{R\}\}\(a\)\+iP\_\{m\}^\{\\mathbb\{R\}\}\(b\)\. This shared projector preserves the complex spectral representation while restricting channel transformation to a low\-rank subspace\.

The projectorPmP\_\{m\}is applied independently to the channel vector at each batch index, time step, and frequency location inℬm\\mathcal\{B\}\_\{m\}\. The projected bands are placed back into their corresponding spectral locations to formδ​z^proj\\widehat\{\\delta z\}\_\{\\mathrm\{proj\}\}, and the projected increment is obtained asδ​zproj=ℱ𝐱−1​\(δ​z^proj\)\\delta z\_\{\\mathrm\{proj\}\}=\\mathcal\{F\}\_\{\\mathbf\{x\}\}^\{\-1\}\(\\widehat\{\\delta z\}\_\{\\mathrm\{proj\}\}\)\. We denote the complete active\-band projection operator byPABP\_\{\\mathrm\{AB\}\}, such thatδ​zproj=PAB​\(δ​zraw\)\\delta z\_\{\\mathrm\{proj\}\}=P\_\{\\mathrm\{AB\}\}\(\\delta z\_\{\\mathrm\{raw\}\}\)\. The projected increment is used as a residual geometric correction:

\(6\)δ​z\\displaystyle\\delta z=δ​zraw\+η​\(PAB​\(δ​zraw\)−δ​zraw\),\\displaystyle=\\delta z\_\{\\mathrm\{raw\}\}\+\\eta\\left\(P\_\{\\mathrm\{AB\}\}\(\\delta z\_\{\\mathrm\{raw\}\}\)\-\\delta z\_\{\\mathrm\{raw\}\}\\right\),\(7\)zt\+1\\displaystyle z\_\{t\+1\}=zt\+δ​z\.\\displaystyle=z\_\{t\}\+\\delta z\.Here,η\\etacontrols the correction strength\. Since this correction is applied before residual state advancement, the active\-band geometry directly participates in the latent transition rather than serving only as an auxiliary regularizer\. The frequency\-dependent latent\-channel geometry underlying this band\-wise design is empirically analyzed in Appendix[B\.2](https://arxiv.org/html/2608.11237#A2.SS2)\.

### 4\.4\.Dual\-Branch Physical\-Space Reconstruction

After the latent update, GeoIncNO obtains the refined incrementδ​z\\delta zand the updated latent statezt\+1z\_\{t\+1\}\. Sinceδ​z\\delta zcaptures the transition\-induced correction whilezt\+1z\_\{t\+1\}represents the evolved latent state, GeoIncNO reconstructs the physical output through two complementary branches\. The input\-anchored correction branch preserves stable input structures with increment\-induced correction, whereas the state\-decoding branch predicts the output from the updated latent state\.

#### 4\.4\.1\.Input\-Anchored Correction Branch

The input\-anchored correction branch maps the input physical field to the target output window and grid usingubase=Aphys​\(uin\)u\_\{\\mathrm\{base\}\}=A\_\{\\mathrm\{phys\}\}\(u\_\{\\mathrm\{in\}\}\), whereAphys=Rsp∘MT∘MCA\_\{\\mathrm\{phys\}\}=R\_\{\\mathrm\{sp\}\}\\circ M\_\{T\}\\circ M\_\{C\}\. Here,MC:ℝCin→ℝCoutM\_\{C\}:\\mathbb\{R\}^\{C\_\{\\mathrm\{in\}\}\}\\rightarrow\\mathbb\{R\}^\{C\_\{\\mathrm\{out\}\}\}maps physical variables,MT:ℝTin→ℝToutM\_\{T\}:\\mathbb\{R\}^\{T\_\{\\mathrm\{in\}\}\}\\rightarrow\\mathbb\{R\}^\{T\_\{\\mathrm\{out\}\}\}maps the temporal window, andRspR\_\{\\mathrm\{sp\}\}resamples the spatial grid from\(N1in,…,Ndsin\)\(N^\{\\mathrm\{in\}\}\_\{1\},\\ldots,N^\{\\mathrm\{in\}\}\_\{d\_\{s\}\}\)to\(N1out,…,Ndsout\)\(N^\{\\mathrm\{out\}\}\_\{1\},\\ldots,N^\{\\mathrm\{out\}\}\_\{d\_\{s\}\}\)\. The refined increment is decoded into a physical correction asδ​u=DΔ​\(δ​z\)\\delta u=D\_\{\\Delta\}\(\\delta z\), whereDΔD\_\{\\Delta\}consists of a pointwise spatiotemporal convolutional head and a temporal linear projection\. The correction branch is then given by

\(8\)ucorr\\displaystyle u\_\{\\mathrm\{corr\}\}=α​ubase\+δ​u,\\displaystyle=\\alpha u\_\{\\mathrm\{base\}\}\+\\delta u,whereα\\alphais a learnable anchor mixing coefficient\.

#### 4\.4\.2\.State\-Decoding Branch

The state\-decoding branch directly predictsustate=Dstate​\(zt\+1\)u\_\{\\mathrm\{state\}\}=D\_\{\\mathrm\{state\}\}\(z\_\{t\+1\}\), whereDstateD\_\{\\mathrm\{state\}\}is independent ofDΔD\_\{\\Delta\}\. WhileDΔ​\(δ​z\)D\_\{\\Delta\}\(\\delta z\)predicts a correction relative to the input\-anchored base field,Dstate​\(zt\+1\)D\_\{\\mathrm\{state\}\}\(z\_\{t\+1\}\)predicts the complete physical field from the evolved latent representation\. Thus, GeoIncNO obtains two complementary predictions,ucorru\_\{\\mathrm\{corr\}\}andustateu\_\{\\mathrm\{state\}\}\. Both decoders map to the physical output space, i\.e\.,DΔ​\(δ​z\),Dstate​\(zt\+1\)∈ℝB×Tout×N1out×⋯×Ndsout×CoutD\_\{\\Delta\}\(\\delta z\),\\,D\_\{\\mathrm\{state\}\}\(z\_\{t\+1\}\)\\in\\mathbb\{R\}^\{B\\times T\_\{\\mathrm\{out\}\}\\times N\_\{1\}^\{\\mathrm\{out\}\}\\times\\cdots\\times N\_\{d\_\{s\}\}^\{\\mathrm\{out\}\}\\times C\_\{\\mathrm\{out\}\}\}\.

### 4\.5\.Mean–Fluctuation Decoupled Reconstruction

Directly fusingucorru\_\{\\mathrm\{corr\}\}andustateu\_\{\\mathrm\{state\}\}mixes mean structure, fluctuation dynamics, and phase errors in the same prediction space\. Motivated by the distinct temporal roles of the mean and zero\-mean fluctuation components illustrated in Appendix[B\.1](https://arxiv.org/html/2608.11237#A2.SS1), GeoIncNO decomposes each candidate prediction into these two components, fuses them separately, and applies phase correction only to the fluctuation component\.

#### 4\.5\.1\.Mean–Fluctuation Decomposition and Fusion

For an output sequenceu∈ℝB×Tout×N1out×⋯×Ndsout×Coutu\\in\\mathbb\{R\}^\{B\\times T\_\{\\mathrm\{out\}\}\\times N^\{\\mathrm\{out\}\}\_\{1\}\\times\\cdots\\times N^\{\\mathrm\{out\}\}\_\{d\_\{s\}\}\\times C\_\{\\mathrm\{out\}\}\}, we define its temporal mean and zero\-mean fluctuation as

\(9\)μ​\(u\)\\displaystyle\\mu\(u\)=1Tout​∑τ=1Toutuτ,\\displaystyle=\\frac\{1\}\{T\_\{\\mathrm\{out\}\}\}\\sum\_\{\\tau=1\}^\{T\_\{\\mathrm\{out\}\}\}u\_\{\\tau\},\(10\)r​\(u\)\\displaystyle r\(u\)=u−repeatt⁡\(μ​\(u\)\)\.\\displaystyle=u\-\\operatorname\{repeat\}\_\{t\}\(\\mu\(u\)\)\.Here,μ​\(u\)∈ℝB×1×N1out×⋯×Ndsout×Cout\\mu\(u\)\\in\\mathbb\{R\}^\{B\\times 1\\times N^\{\\mathrm\{out\}\}\_\{1\}\\times\\cdots\\times N^\{\\mathrm\{out\}\}\_\{d\_\{s\}\}\\times C\_\{\\mathrm\{out\}\}\}, andrepeatt⁡\(μ​\(u\)\)\\operatorname\{repeat\}\_\{t\}\(\\mu\(u\)\)copies the temporal mean along the output time dimension\. By construction,meant⁡\(r​\(u\)\)=0\\operatorname\{mean\}\_\{t\}\(r\(u\)\)=0, andu=repeatt⁡\(μ​\(u\)\)\+r​\(u\)u=\\operatorname\{repeat\}\_\{t\}\(\\mu\(u\)\)\+r\(u\)\.

Applying this decomposition to the two candidate predictions gives

\(11\)μcorr\\displaystyle\\mu\_\{\\mathrm\{corr\}\}=μ​\(ucorr\),\\displaystyle=\\mu\(u\_\{\\mathrm\{corr\}\}\),rcorr\\displaystyle r\_\{\\mathrm\{corr\}\}=ucorr−repeatt⁡\(μcorr\),\\displaystyle=u\_\{\\mathrm\{corr\}\}\-\\operatorname\{repeat\}\_\{t\}\(\\mu\_\{\\mathrm\{corr\}\}\),\(12\)μstate\\displaystyle\\mu\_\{\\mathrm\{state\}\}=μ​\(ustate\),\\displaystyle=\\mu\(u\_\{\\mathrm\{state\}\}\),rstate\\displaystyle r\_\{\\mathrm\{state\}\}=ustate−repeatt⁡\(μstate\)\.\\displaystyle=u\_\{\\mathrm\{state\}\}\-\\operatorname\{repeat\}\_\{t\}\(\\mu\_\{\\mathrm\{state\}\}\)\.GeoIncNO uses separate gates for mean\-field and fluctuation\-field fusion:

\(13\)μ^\\displaystyle\\hat\{\\mu\}=Γμ​\(zt\+1\)⊙μcorr\+\[1−Γμ​\(zt\+1\)\]⊙μstate,\\displaystyle=\\Gamma\_\{\\mu\}\(z\_\{t\+1\}\)\\odot\\mu\_\{\\mathrm\{corr\}\}\+\\left\[1\-\\Gamma\_\{\\mu\}\(z\_\{t\+1\}\)\\right\]\\odot\\mu\_\{\\mathrm\{state\}\},\(14\)r^\\displaystyle\\hat\{r\}=centert⁡\(Γr​\(zt\+1\)⊙rcorr\+\[1−Γr​\(zt\+1\)\]⊙rstate\)\.\\displaystyle=\\operatorname\{center\}\_\{t\}\\left\(\\Gamma\_\{r\}\(z\_\{t\+1\}\)\\odot r\_\{\\mathrm\{corr\}\}\+\\left\[1\-\\Gamma\_\{r\}\(z\_\{t\+1\}\)\\right\]\\odot r\_\{\\mathrm\{state\}\}\\right\)\.Here,Γμ\\Gamma\_\{\\mu\}andΓr\\Gamma\_\{r\}are generated from the updated latent state, andcentert⁡\(⋅\)\\operatorname\{center\}\_\{t\}\(\\cdot\)removes the temporal mean to keepr^\\hat\{r\}zero\-mean\. The separate gates allow the mean structure and fluctuation dynamics to use different branch preferences\.

#### 4\.5\.2\.Fluctuation\-Only Phase Correction

Phase shift is meaningful for non\-zero temporal fluctuations rather than for the temporal mean\. GeoIncNO therefore applies phase correction only to the fused fluctuationr^\\hat\{r\}\. Let𝒫phase\\mathcal\{P\}\_\{\\mathrm\{phase\}\}denote a lightweight temporal advancement operator\. The phase residual and gated correction are

\(15\)Δ​rphase\\displaystyle\\Delta r\_\{\\mathrm\{phase\}\}=𝒫phase​\(r^\)−r^,\\displaystyle=\\mathcal\{P\}\_\{\\mathrm\{phase\}\}\(\\hat\{r\}\)\-\\hat\{r\},\(16\)rphase\\displaystyle r\_\{\\mathrm\{phase\}\}=centert⁡\(Γphase​\(zt\+1\)⊙Δ​rphase\)\.\\displaystyle=\\operatorname\{center\}\_\{t\}\\left\(\\Gamma\_\{\\mathrm\{phase\}\}\(z\_\{t\+1\}\)\\odot\\Delta r\_\{\\mathrm\{phase\}\}\\right\)\.The centering operation keeps the phase correction zero\-mean and prevents direct perturbation of the temporal mean\.

The final prediction is reconstructed as

\(17\)u^\\displaystyle\\hat\{u\}=repeatt⁡\(μ^\)\+r^\+rphase\.\\displaystyle=\\operatorname\{repeat\}\_\{t\}\(\\hat\{\\mu\}\)\+\\hat\{r\}\+r\_\{\\mathrm\{phase\}\}\.Thus, the reconstruction separates stable mean structures, non\-stationary fluctuations, and phase correction within the fluctuation subspace\.

### 4\.6\.Training Objective

GeoIncNO is trained to jointly promote accurate physical\-field prediction, structured active\-band increment geometry, and limited deviation between the raw and geometrically shaped increments\. The primary prediction loss isℒpred=MSE⁡\(u^,uout\)\\mathcal\{L\}\_\{\\mathrm\{pred\}\}=\\operatorname\{MSE\}\(\\hat\{u\},u\_\{\\mathrm\{out\}\}\)\.

The frequency–channel structure of the latent increment is regulated through an active\-band geometry loss\. For themm\-th active bandℬm\\mathcal\{B\}\_\{m\}, we flatten the batch, temporal, and frequency dimensions ofδ​z^raw\(m\)\\widehat\{\\delta z\}\_\{\\mathrm\{raw\}\}^\{\(m\)\}and concatenate its real and imaginary components, yieldingXm∈ℝNm×dX\_\{m\}\\in\\mathbb\{R\}^\{N\_\{m\}\\times d\}\. Its centered covariance and normalized correlation matrices are

\(18\)Cm\\displaystyle C\_\{m\}=1Nm−1​X¯m⊤​X¯m,\\displaystyle=\\frac\{1\}\{N\_\{m\}\-1\}\\bar\{X\}\_\{m\}^\{\\top\}\\bar\{X\}\_\{m\},\(19\)Rm\\displaystyle R\_\{m\}=Dm−1/2​Cm​Dm−1/2,\\displaystyle=D\_\{m\}^\{\-1/2\}C\_\{m\}D\_\{m\}^\{\-1/2\},\(20\)Dm\\displaystyle D\_\{m\}=diag⁡\(Cm\)\+ϵ​Id,\\displaystyle=\\operatorname\{diag\}\(C\_\{m\}\)\+\\epsilon I\_\{d\},whereX¯m\\bar\{X\}\_\{m\}is centered along the sample dimension\. Redundant channel coupling is penalized by

\(21\)ℒcorr\\displaystyle\\mathcal\{L\}\_\{\\mathrm\{corr\}\}=1M​∑m=1M‖offdiag⁡\(Rm\)‖F2d​\(d−1\)\.\\displaystyle=\\frac\{1\}\{M\}\\sum\_\{m=1\}^\{M\}\\frac\{\\left\\\|\\operatorname\{offdiag\}\(R\_\{m\}\)\\right\\\|\_\{F\}^\{2\}\}\{d\(d\-1\)\}\.This term suppresses strong off\-diagonal correlations and encourages different channels to represent complementary transition directions\. A variance support term is further introduced to prevent feature collapse:

\(22\)ℒvar\\displaystyle\\mathcal\{L\}\_\{\\mathrm\{var\}\}=1M​∑m=1M1d​∑c=1d\[max⁡\(0,γ−Cm​\(c,c\)\+ϵ\)\]2\.\\displaystyle=\\frac\{1\}\{M\}\\sum\_\{m=1\}^\{M\}\\frac\{1\}\{d\}\\sum\_\{c=1\}^\{d\}\\left\[\\max\\left\(0,\\gamma\-\\sqrt\{C\_\{m\}\(c,c\)\+\\epsilon\}\\right\)\\right\]^\{2\}\.This term maintains sufficient channel variation within each active band\. The complete geometry loss isℒgeom=ℒcorr\+ℒvar\\mathcal\{L\}\_\{\\mathrm\{geom\}\}=\\mathcal\{L\}\_\{\\mathrm\{corr\}\}\+\\mathcal\{L\}\_\{\\mathrm\{var\}\}\. Projection consistency is enforced by keeping the refined increment close to the transition direction predicted by the latent backbone:

\(23\)ℒproj\\displaystyle\\mathcal\{L\}\_\{\\mathrm\{proj\}\}=‖δ​z−δ​zraw‖22‖δ​zraw‖22\+ϵ\.\\displaystyle=\\frac\{\\left\\\|\\delta z\-\\delta z\_\{\\mathrm\{raw\}\}\\right\\\|\_\{2\}^\{2\}\}\{\\left\\\|\\delta z\_\{\\mathrm\{raw\}\}\\right\\\|\_\{2\}^\{2\}\+\\epsilon\}\.
The final objective is

\(24\)ℒ\\displaystyle\\mathcal\{L\}=ℒpred\+λgeom​\(s\)​ℒgeom\+λproj​\(s\)​ℒproj,\\displaystyle=\\mathcal\{L\}\_\{\\mathrm\{pred\}\}\+\\lambda\_\{\\mathrm\{geom\}\}\(s\)\\mathcal\{L\}\_\{\\mathrm\{geom\}\}\+\\lambda\_\{\\mathrm\{proj\}\}\(s\)\\mathcal\{L\}\_\{\\mathrm\{proj\}\},wheressdenotes the training progress\. A warm\-up schedule gradually increasesλgeom​\(s\)\\lambda\_\{\\mathrm\{geom\}\}\(s\)andλproj​\(s\)\\lambda\_\{\\mathrm\{proj\}\}\(s\), allowing GeoIncNO to first learn a reliable predictive transition and then progressively enforce increment geometry and projection consistency\.

The complete inference procedure is summarized in Algorithm[1](https://arxiv.org/html/2608.11237#alg1)\. Appendix[C\.1](https://arxiv.org/html/2608.11237#A3.SS1)provides a conditional interpretation of GeoIncNO\. Under bounded latent reconstruction error and Lipschitz decoding, controlling the latent increment error also controls the physical prediction error; under spectral\-support and band\-wise low\-rank assumptions, ABP admits a geometric interpretation\.

## 5\.Experiments

### 5\.1\.Experimental Setup

#### 5\.1\.1\.Baseline Models

We select seven baselines to cover neural\-operator paradigms for PDE prediction\. DeepONet\(Lu et al\.,[2021](https://arxiv.org/html/2608.11237#bib.bib26)\)is included as a classical operator\-learning baseline based on the branch–trunk formulation\. FNO\(Li et al\.,[2021a](https://arxiv.org/html/2608.11237#bib.bib22)\), UNO\(Rahman et al\.,[2023](https://arxiv.org/html/2608.11237#bib.bib30)\), and WNO\(Tripura and Chakraborty,[2023](https://arxiv.org/html/2608.11237#bib.bib36)\)represent transform\-domain and multi\-scale neural operators, covering Fourier\-based global spectral modeling, U\-shaped multi\-resolution operator learning, and wavelet\-based localized multi\-scale modeling\. We also include PINO\(Li et al\.,[2021b](https://arxiv.org/html/2608.11237#bib.bib24)\)as a physics\-informed neural operator that incorporates PDE residual constraints during training\. We further include two recent latent\-space neural operators, LNO\(Wang and Wang,[2024](https://arxiv.org/html/2608.11237#bib.bib37)\)and LaMO\(Tiwari et al\.,[2025](https://arxiv.org/html/2608.11237#bib.bib35)\), as the most direct baselines to GeoIncNO\. While these methods also reduce PDE prediction to latent\-space operator learning, they primarily focus on learning compact latent representations or efficient latent dynamics\. In contrast, GeoIncNO explicitly formulates the latent transition as an increment and shapes its active spectral components before state advancement\.

#### 5\.1\.2\.PDE Benchmarks

We evaluate GeoIncNO on six PDE benchmarks generated under a PDEBench\-style protocol: Burgers, Kuramoto–Sivashinsky \(KS\), Navier–Stokes \(NS\), shallow\-water \(SW\), 3D compressible Euler \(CE\), and Maxwell\. These benchmarks span one\-, two\-, and three\-dimensional systems and cover convection–diffusion, spatiotemporally chaotic, vortical, wave\-propagation, compressible\-flow, and electromagnetic dynamics\. Table[3](https://arxiv.org/html/2608.11237#A4.T3)summarizes their spatial domains, grid resolutions, temporal ranges, numbers of snapshots, and numbers of physical state variables\. In the Time Range column, the notation\[tmin,tmax\]/Nt\[t\_\{\\min\},t\_\{\\max\}\]/N\_\{t\}denotes the simulated time interval and the number of temporal snapshots\. The corresponding governing equations and variable definitions are provided in Appendix[D\.1\.1](https://arxiv.org/html/2608.11237#A4.SS1.SSS1)\.

#### 5\.1\.3\.Evaluation Metrics

We report metrics averaged over the test set to evaluate prediction accuracy and rollout stability\. Field\-level accuracy is measured by the relativeL2L\_\{2\}errorRel​\-​L2=‖u−u^‖2/‖u‖2\\mathrm\{Rel\}\\text\{\-\}L\_\{2\}=\\\|u\-\\hat\{u\}\\\|\_\{2\}/\\\|u\\\|\_\{2\}and MSE, where Rel\-L2L\_\{2\}evaluates scale\-normalized errors and MSE measures absolute deviations\. We further report the relative Sobolev errorRel​\-​H1=‖u−u^‖H1/‖u‖H1\\mathrm\{Rel\}\\text\{\-\}H^\{1\}=\\\|u\-\\hat\{u\}\\\|\_\{H^\{1\}\}/\\\|u\\\|\_\{H^\{1\}\}, which accounts for both field and gradient discrepancies\. For spectral accuracy, we use the weighted log ratio \(WLR\) to measure discrepancies between predicted and reference spectral energies\. For autoregressive rollouts, we report the final\-step RMSE \(F\-RMSE\) and cumulative RMSE \(C\-RMSE\), which measure terminal and accumulated prediction errors, respectively\. Boundary RMSE \(B\-RMSE\) is used when boundary\-region accuracy is evaluated\. We also report the relative temporal\-mean drift \(Mean Drift\) to quantify long\-term bias in the rollout trajectory\. The full definitions of WLR and Mean Drift are provided in Appendix[D\.1\.2](https://arxiv.org/html/2608.11237#A4.SS1.SSS2)\.

#### 5\.1\.4\.Implementation Details

All experiments are implemented in PyTorch on a single NVIDIA RTX 4090 GPU \(24 GB\)\. For each benchmark, we use 1,000 trajectories for training and 200 independently generated trajectories for testing\. All non\-physics\-informed baselines are trained with the MSE prediction loss only, whereas PINO additionally uses the PDE residual loss following its original formulation, and GeoIncNO is optimized with the objective in Section[4\.6](https://arxiv.org/html/2608.11237#S4.SS6), which augments the prediction loss with active\-band increment geometry regularization and projection consistency\. All models are trained for 100 epochs using Adam\(Kingma and Ba,[2015](https://arxiv.org/html/2608.11237#bib.bib19)\)with an initial learning rate of1×10−41\\times 10^\{\-4\}and cosine annealing, together with the pushforward strategy at a rollout length ofT=5T=5for long\-horizon stability\. For GeoIncNO, we useM=4M=4active frequency bands, a low\-rank projector dimension ofr=8r=8, and an ABP correction strength ofη=0\.05\\eta=0\.05unless otherwise specified\. Full training and dimension\-specific operator details are provided in Appendix[D\.1\.3](https://arxiv.org/html/2608.11237#A4.SS1.SSS3)\.

### 5\.2\.Overall Prediction Performance

Table[1](https://arxiv.org/html/2608.11237#S5.T1)compares GeoIncNO with competitive neural\-operator baselines across six PDE benchmarks in field accuracy, local\-structure preservation, spectral fidelity, and cumulative rollout error\. The complete results, including MSE, F\-RMSE, parameter counts, and standard deviations, are provided in Appendix Table[4](https://arxiv.org/html/2608.11237#A4.T4)\. GeoIncNO achieves the best results on all four primary metrics across all benchmarks\. Compared with the FNO backbone, GeoIncNO substantially reduces Rel\-L2L\_\{2\}, from3\.22​E\-​23\.22\\text\{E\-\}2to1\.64​E\-​21\.64\\text\{E\-\}2on 1D KS, from2\.65​E\-​12\.65\\text\{E\-\}1to7\.54​E\-​37\.54\\text\{E\-\}3on 2D NS, and from2\.75​E\-​22\.75\\text\{E\-\}2to5\.43​E\-​35\.43\\text\{E\-\}3on 3D CE\. The gains are particularly evident on the chaotic 1D KS benchmark, where several baselines exhibit severe error accumulation, whereas GeoIncNO achieves a C\-RMSE of2\.02​E\-​12\.02\\text\{E\-\}1and a Rel\-H1H^\{1\}of2\.12​E\-​22\.12\\text\{E\-\}2\. GeoIncNO also consistently outperforms the latent\-space baselines LNO and LaMO\. On 2D NS, for example, it reduces LaMO’s Rel\-L2L\_\{2\}, Rel\-H1H^\{1\}, WLR, and C\-RMSE from1\.37​E\-​21\.37\\text\{E\-\}2,2\.42​E\-​22\.42\\text\{E\-\}2,2\.05​E\-​22\.05\\text\{E\-\}2, and1\.43​E\-​11\.43\\text\{E\-\}1to7\.54​E\-​37\.54\\text\{E\-\}3,1\.33​E\-​21\.33\\text\{E\-\}2,1\.13​E\-​21\.13\\text\{E\-\}2, and7\.87​E\-​27\.87\\text\{E\-\}2, respectively\. The simultaneous improvements in Rel\-H1H^\{1\}, WLR, and C\-RMSE indicate that GeoIncNO better preserves local and spectral structures throughout autoregressive rollout, supporting the effectiveness of explicitly structuring the latent transition increment\.

Table 1\.Main prediction results across six PDE benchmarks
### 5\.3\.Long\-Horizon Rollout Stability

To evaluate the long\-horizon stability of GeoIncNO, we conduct autoregressive rollout experiments across all six PDE benchmarks\. Starting from the initial condition, each model is applied autoregressively, where the prediction at the current step is fed back as the input for the next step\. We report MSE, Rel\-L2L\_\{2\}, and Rel\-H1H^\{1\}at five checkpoints,0\.2​T0\.2T–1\.0​T1\.0T, to measure how errors accumulate over the horizon\. All models are evaluated under the same rollout protocol, initial conditions, prediction horizon, and evaluation metrics\. As shown in Figure[7](https://arxiv.org/html/2608.11237#A4.F7)\(in Appendix[D\.2](https://arxiv.org/html/2608.11237#A4.SS2)\), GeoIncNO attains the lowest error at every checkpoint on all three metrics across the 1D, 2D, and 3D benchmarks, and its advantage widens with the horizon\. It stays nearly flat on the chaotic 1D KS benchmark where several baselines diverge, and preserves gradient\-level structure on 2D NS and SW where competing operators inflate Rel\-H1H^\{1\}\. Compared with the FNO backbone and the closest latent\-space baseline LaMO, GeoIncNO shows a markedly slower error\-accumulation slope rather than merely a lower one\-step error, so the gap is largest at the late checkpoints where long\-horizon stability matters most\. Overall, these results indicate that explicitly structuring the latent transition increment, rather than repeatedly applying an unconstrained state\-to\-state update, suppresses the high\-frequency error accumulation that destabilizes autoregressive rollout\.

### 5\.4\.Ablation Study

To examine whether each component contributes to the long\-horizon prediction behavior of GeoIncNO, we conduct component ablations on the 2D Navier–Stokes \(NS\) benchmark\. We define the variants as follows: i\)Base LatentNOcontains the encoder, latent backbone, residual latent\-state update, and decoder, without geometric shaping or structured reconstruction; ii\)\+ Fixed\-band LIGadds latent increment geometry \(LIG\), where low\-rank projectors shape the latent increment over fixed frequency bands; iii\)\+ Active Bandsreplaces fixed bands with spectrum\-adaptive active bands constructed from the increment spectral distribution; iv\)\+ MFDRadds mean–fluctuation decoupled reconstruction, which separately fuses temporal mean and zero\-mean fluctuation components; v\)\+ FFPCapplies full\-field phase correction to the reconstructed prediction; and vi\)GeoIncNOkeeps the same components but restricts phase correction to the zero\-mean fluctuation component\. Each row in Table[2](https://arxiv.org/html/2608.11237#S5.T2)is cumulative with respect to the previous row unless otherwise specified\.

Table 2\.Cumulative component ablation on the 2D NSTable[2](https://arxiv.org/html/2608.11237#S5.T2)shows monotone gains from each component\. Fixed increment geometry already improves over the unconstrained baseline \(Rel\-L2L\_\{2\}: 2\.18E\-02 to 1\.78E\-02\), and active bands reduce all four metrics further, supporting spectrum\-adaptive shaping\. MFDR produces the largest Mean Drift drop, from 1\.32E\-03 to 1\.05E\-04, while also lowering prediction and spectral error\. Full\-field phase correction lowers Rel\-L2L\_\{2\}but raises Mean Drift back to 3\.42E\-04, whereas restricting it to the zero\-mean fluctuation \(GeoIncNO\) reaches the best value on every metric and drives Mean Drift to 6\.21E\-05\. Overall, active\-band increment geometry and mean–fluctuation decoupling jointly improve spectral fidelity and rollout stability while confining phase correction to the fluctuation subspace preserves the temporal mean\.

### 5\.5\.Additional Results

Additional experiments further evaluate GeoIncNO beyond the standard benchmark setting\. Cross\-parameter generalization to an unseen Reynolds number is reported in Appendix[D\.3](https://arxiv.org/html/2608.11237#A4.SS3), while cross\-resolution transfer from coarse to finer grids is examined in Appendix[D\.5](https://arxiv.org/html/2608.11237#A4.SS5)\. Extended rollout behavior under the unseen parameter regime is analyzed in Appendix[D\.4](https://arxiv.org/html/2608.11237#A4.SS4)\. Results on the RealPDEBench fluid–structure interaction task are provided in Appendix[D\.6](https://arxiv.org/html/2608.11237#A4.SS6), and qualitative comparisons across competitive 1D, 2D, and 3D systems are shown in Appendix[D\.7](https://arxiv.org/html/2608.11237#A4.SS7)\.

## 6\.Conclusion

We proposed GeoIncNO for stable long\-horizon PDE prediction\. GeoIncNO treats the latent transition increment as the object governing autoregressive evolution, rather than only improving the latent state representation or direct state\-to\-state mapping\. This formulation allows the repeated dynamical update to be regularized before residual state advancement\. Specifically, active\-band low\-rank projection shapes the spectral\-channel geometry of latent increments, and mean–fluctuation decoupled reconstruction separates stable mean structures, dynamic fluctuations, and fluctuation\-level phase correction in the physical output space\. Experiments on six one\-, two\-, and three\-dimensional PDE benchmarks have shown that GeoIncNO improves prediction accuracy, spectral fidelity, and rollout stability over competitive neural\-operator baselines\. Ablation, cross\-parameter, cross\-resolution, and qualitative results have further confirmed that structured latent increments and mean–fluctuation decoupling reduce error accumulation during autoregressive prediction\.

## 7\.Limitations and Ethical Considerations

GeoIncNO currently calibrates active bands from training\-set spectral statistics and keeps them fixed during inference\. Future work may explore adaptive band construction and broader validation under more diverse geometries, boundary conditions, and distribution shifts\. In addition, incorporating physical constraints or uncertainty estimation could further improve reliability for very long rollouts\. This work uses no human subjects, personal information, or privacy\-sensitive data\.

## 8\.Generative AI Usage

Generative AI was used only for language editing\.

## References

- \(1\)
- Aarts and Van Der Veer \(2001\)Lucie P Aarts and Peter Van Der Veer\. 2001\.Neural network method for solving partial differential equations\.*Neural processing letters*14, 3 \(2001\), 261–271\.
- Azizzadenesheli et al\.\(2024\)Kamyar Azizzadenesheli, Nikola Kovachki, Zongyi Li, Miguel Liu\-Schiaffini, Jean Kossaifi, and Anima Anandkumar\. 2024\.Neural operators for accelerating scientific simulations and design\.*Nature Reviews Physics*6, 5 \(2024\), 320–328\.
- Brandstetter et al\.\(2022\)Johannes Brandstetter, Daniel Worrall, and Max Welling\. 2022\.Message passing neural PDE solvers\.*arXiv preprint arXiv:2202\.03376*\(2022\)\.
- Brunton et al\.\(2022\)S\. L\. Brunton, M\. Budišić, E\. Kaiser, and J\. N\. Kutz\. 2022\.Modern Koopman theory for dynamical systems\.*SIAM Rev\.*64, 2 \(2022\), 229–340\.
- Brunton and Kutz \(2024\)Steven L Brunton and J Nathan Kutz\. 2024\.Promising directions of machine learning for partial differential equations\.*Nature Computational Science*4, 7 \(2024\), 483–494\.
- Cao \(2021\)Shuhao Cao\. 2021\.Choose a transformer: Fourier or galerkin\.*Advances in neural information processing systems*34 \(2021\), 24924–24940\.
- Chen and Wu \(2025\)Chuanqi Chen and Jin\-Long Wu\. 2025\.Neural dynamical operator: Continuous spatial\-temporal model with gradient\-based and derivative\-free optimization methods\.*J\. Comput\. Phys\.*520 \(2025\), 113480\.
- Dai et al\.\(2026\)Yilong Dai, Shengyu Chen, Ziyi Wang, Xiaowei Jia, Yiqun Xie, Vipin Kumar, and Runlong Yu\. 2026\.Learning PDE Solvers with Physics and Data: A Unifying View of Physics\-Informed Neural Networks and Neural Operators\.*arXiv preprint arXiv:2601\.14517*\(2026\)\.
- Guibas et al\.\(2022\)J\. Guibas, M\. Mardani, Z\. Li, A\. Tao, A\. Anandkumar, and B\. Catanzaro\. 2022\.Adaptive Fourier neural operators: Efficient token mixers for transformers\. In*Advances in Neural Information Processing Systems \(NeurIPS\)*\.
- Gupta et al\.\(2021\)G\. Gupta, X\. Xiao, and P\. Bogdan\. 2021\.Multiwavelet\-based operator learning for differential equations\. In*Advances in Neural Information Processing Systems \(NeurIPS\)*\.
- Hao et al\.\(2024\)Zhongkai Hao, Chang Su, Songming Liu, Julius Berner, Chengyang Ying, Hang Su, Anima Anandkumar, Jian Song, and Jun Zhu\. 2024\.Dpot: Auto\-regressive denoising operator transformer for large\-scale pde pre\-training\.*arXiv preprint arXiv:2403\.03542*\(2024\)\.
- Hao et al\.\(2023\)Z\. Hao, C\. Ying, Z\. Wang, H\. Su, Y\. Dong, S\. Liu, Z\. Cheng, J\. Zhu, and J\. Song\. 2023\.GNOT: A general neural operator transformer for operator learning\. In*International Conference on Machine Learning \(ICML\)*\.
- Hu et al\.\(2026\)Peiyan Hu, Haodong Feng, Hongyuan Liu, Tongtong Yan, Wenhao Deng, Tianrun Gao, Rong Zheng, Haoren Zheng, Chenglei Yu, Chuanrui Wang, et al\.2026\.RealPDEBench: A Benchmark for Complex Physical Systems with Real\-World Data\.*arXiv preprint arXiv:2601\.01829*\(2026\)\.
- Hu et al\.\(2025\)Peiyan Hu, Rui Wang, Xiang Zheng, Tao Zhang, Haodong Feng, Ruiqi Feng, Long Wei, Yue Wang, Zhi\-Ming Ma, and Tailin Wu\. 2025\.Wavelet diffusion neural operator\. In*International Conference on Learning Representations*, Vol\. 2025\. 12291–12333\.
- Hu et al\.\(2024\)Peiyan Hu, Yue Wang, and Zhi\-Ming Ma\. 2024\.Better neural PDE solvers through data\-free mesh movers\. In*International Conference on Learning Representations*, Vol\. 2024\. 4550–4576\.
- Karniadakis et al\.\(2021\)George Em Karniadakis, Ioannis G Kevrekidis, Lu Lu, Paris Perdikaris, Sifan Wang, and Liu Yang\. 2021\.Physics\-informed machine learning\.*Nature Reviews Physics*3, 6 \(2021\), 422–440\.
- Karumuri et al\.\(2026\)Sharmila Karumuri, Lori Graham\-Brady, and Somdatta Goswami\. 2026\.Physics\-informed latent neural operator for real\-time predictions of time\-dependent parametric PDEs\.*Computer Methods in Applied Mechanics and Engineering*450 \(2026\), 118599\.
- Kingma and Ba \(2015\)Diederik P\. Kingma and Jimmy Ba\. 2015\.Adam: A Method for Stochastic Optimization\. In*3rd International Conference on Learning Representations, ICLR 2015, San Diego, CA, USA, May 7\-9, 2015, Conference Track Proceedings*, Yoshua Bengio and Yann LeCun \(Eds\.\)\.[http://arxiv\.org/abs/1412\.6980](http://arxiv.org/abs/1412.6980)
- Kovachki et al\.\(2023\)N\. Kovachki, Z\. Li, B\. Liu, K\. Azizzadenesheli, K\. Bhattacharya, A\. Stuart, and A\. Anandkumar\. 2023\.Neural operator: Learning maps between function spaces with applications to PDEs\.*Journal of Machine Learning Research*24, 89 \(2023\), 1–97\.
- Kutz et al\.\(2016\)J Nathan Kutz, Steven L Brunton, Bingni W Brunton, and Joshua L Proctor\. 2016\.*Dynamic mode decomposition: data\-driven modeling of complex systems*\.SIAM\.
- Li et al\.\(2021a\)Zongyi Li, Nikola Borislavov Kovachki, Kamyar Azizzadenesheli, Burigede Liu, Kaushik Bhattacharya, Andrew M\. Stuart, and Anima Anandkumar\. 2021a\.Fourier Neural Operator for Parametric Partial Differential Equations\. In*9th International Conference on Learning Representations, ICLR 2021, Virtual Event, Austria, May 3\-7, 2021*\.
- Li et al\.\(2023\)Z\. Li, K\. Meidani, and A\. B\. Farimani\. 2023\.Transformer for partial differential equations’ operator learning\.*Transactions on Machine Learning Research \(TMLR\)*\(2023\)\.
- Li et al\.\(2021b\)Zongyi Li, Hongkai Zheng, Nikola B\. Kovachki, David Jin, Haoxuan Chen, Burigede Liu, Kamyar Azizzadenesheli, and Anima Anandkumar\. 2021b\.Physics\-Informed Neural Operator for Learning Partial Differential Equations\.*CoRR*abs/2111\.03794 \(2021\)\.
- Liu et al\.\(2025\)Shengjun Liu, Yu Yu, Ting Zhang, Hanchao Liu, Xinru Liu, and Deyu Meng\. 2025\.Architectures, variants, and performance of neural operators: A comparative review\.*Neurocomputing*648 \(2025\), 130518\.
- Lu et al\.\(2021\)L\. Lu, P\. Jin, G\. Pang, Z\. Zhang, and G\. E\. Karniadakis\. 2021\.Learning nonlinear operators via DeepONet based on the universal approximation theorem of operators\.*Nature Machine Intelligence*3 \(2021\), 218–229\.
- Lusch et al\.\(2018\)Bethany Lusch, J Nathan Kutz, and Steven L Brunton\. 2018\.Deep learning for universal linear embeddings of nonlinear dynamics\.*Nature communications*9, 1 \(2018\), 4950\.
- McCabe et al\.\(2023\)Michael McCabe, Peter Harrington, Shashank Subramanian, and Jed Brown\. 2023\.Towards stability of autoregressive neural operators\.*arXiv preprint arXiv:2306\.10619*\(2023\)\.
- Morton et al\.\(2018\)J\. Morton, A\. Jameson, M\. J\. Kochenderfer, and F\. D\. Witherden\. 2018\.Deep dynamical modeling and control of unsteady fluid flows\.*Advances in Neural Information Processing Systems \(NeurIPS\)*\(2018\)\.
- Rahman et al\.\(2023\)Md Ashiqur Rahman, Zachary E\. Ross, and Kamyar Azizzadenesheli\. 2023\.U\-NO: U\-shaped Neural Operators\.*Trans\. Mach\. Learn\. Res\.*2023 \(2023\)\.
- Raissi et al\.\(2019\)M\. Raissi, P\. Perdikaris, and G\. E\. Karniadakis\. 2019\.Physics\-informed neural networks: A deep learning framework for solving forward and inverse problems involving nonlinear partial differential equations\.*J\. Comput\. Phys\.*378 \(2019\), 686–707\.
- Raonic et al\.\(2023\)Bogdan Raonic, Roberto Molinaro, Tim De Ryck, Tobias Rohner, Francesca Bartolucci, Rima Alaifari, Siddhartha Mishra, and Emmanuel De Bézenac\. 2023\.Convolutional neural operators for robust and accurate learning of pdes\.*Advances in Neural Information Processing Systems*36 \(2023\), 77187–77200\.
- Ren et al\.\(2026\)Xingyu Ren, Pengkai Wang, Pengwei Liu, Xihang Yue, Huanshuo Dong, Zhenxin Huang, Zhongkai Hao, Ziqian Hu, Zhen Huang, Yian Wang, et al\.2026\.Foundation neural operators: A survey on pretraining methods, the data ecosystem, and efficient adaptation\.\(2026\)\.
- Ronneberger et al\.\(2015\)Olaf Ronneberger, Philipp Fischer, and Thomas Brox\. 2015\.U\-net: Convolutional networks for biomedical image segmentation\. In*International Conference on Medical image computing and computer\-assisted intervention*\. Springer, 234–241\.
- Tiwari et al\.\(2025\)Karn Tiwari, Niladri Dutta, NM Krishnan, et al\.2025\.Latent mamba operator for partial differential equations\.*arXiv preprint arXiv:2505\.19105*\(2025\)\.
- Tripura and Chakraborty \(2023\)Tapas Tripura and Souvik Chakraborty\. 2023\.Wavelet neural operator for solving parametric partial differential equations in computational mechanics problems\.*Computer Methods in Applied Mechanics and Engineering*404 \(2023\), 115783\.
- Wang and Wang \(2024\)Tian Wang and Chuang Wang\. 2024\.Latent neural operator for solving forward and inverse pde problems\.*Advances in Neural Information Processing Systems*37 \(2024\), 33085–33107\.
- Wen et al\.\(2022\)G\. Wen, Z\. Li, K\. Azizzadenesheli, A\. Anandkumar, and S\. M\. Benson\. 2022\.U\-FNO—An enhanced Fourier neural operator\-based deep\-learning model for multiphase flow\.*Advances in Water Resources*163 \(2022\), 104180\.
- Worrall et al\.\(2024\)Daniel E Worrall, Miles Cranmer, J Nathan Kutz, and Peter Battaglia\. 2024\.Spectral Shaping for Neural PDE Surrogates\.\(2024\)\.
- Wu et al\.\(2025\)Hao Wu, Yuan Gao, Fan Xu, Fan Zhang, Qingsong Wen, Kun Wang, Xiaomeng Huang, and Xian Wu\. 2025\.Differential\-Integral Neural Operator for Long\-Term Turbulence Forecasting\.*arXiv preprint arXiv:2509\.21196*\(2025\)\.
- Wu et al\.\(2023\)H\. Wu, T\. Hu, H\. Luo, J\. Wang, and M\. Long\. 2023\.Solving high\-dimensional PDEs with latent spectral models\. In*International Conference on Machine Learning \(ICML\)*\.
- Wu et al\.\(2024\)Haixu Wu, Huakun Luo, Haowen Wang, Jianmin Wang, and Mingsheng Long\. 2024\.Transolver: A fast transformer solver for pdes on general geometries\.*arXiv preprint arXiv:2402\.02366*\(2024\)\.
- Xiong et al\.\(2024\)W\. Xiong, X\. Huang, Z\. Zhang, R\. Deng, P\. Sun, and Y\. Tian\. 2024\.Koopman neural operator as a mesh\-free solver of non\-linear partial differential equations\.*J\. Comput\. Phys\.*513 \(2024\), 113194\.
- Yip et al\.\(2022\)Chung Yip, Lokyi Seol, and Xu Zhu Hon\. 2022\.A Step\-by\-Step Approach to Partial Differential Equations\.*Fusion of Multidisciplinary Research, An International Journal*3, 1 \(2022\), 302–315\.
- You et al\.\(2025\)Z\. You, Z\. Xu, and W\. Cai\. 2025\.MscaleFNO: Multi\-scale Fourier neural operator learning for oscillatory functions and wave scattering problems\.*J\. Comput\. Phys\.*\(2025\)\.in press\.

## Appendix ARelated Work

### A\.1\.Neural Operators and Latent Dynamics

Unlike physics\-informed neural networks\(Raissi et al\.,[2019](https://arxiv.org/html/2608.11237#bib.bib31); Karniadakis et al\.,[2021](https://arxiv.org/html/2608.11237#bib.bib17)\), which embed governing equations into the loss and are typically retrained per instance, neural operators\(Kovachki et al\.,[2023](https://arxiv.org/html/2608.11237#bib.bib20)\)learn mappings between function spaces, with FNO\(Li et al\.,[2021a](https://arxiv.org/html/2608.11237#bib.bib22)\)and DeepONet\(Lu et al\.,[2021](https://arxiv.org/html/2608.11237#bib.bib26)\)as the two representative architectures\. Subsequent work follows two routes\. The first route strengthens the operator backbone: multi\-scale and adaptive Fourier mixing\(Wen et al\.,[2022](https://arxiv.org/html/2608.11237#bib.bib38); Guibas et al\.,[2022](https://arxiv.org/html/2608.11237#bib.bib10)\), spectral\-basis enrichment\(Gupta et al\.,[2021](https://arxiv.org/html/2608.11237#bib.bib11); Hu et al\.,[2025](https://arxiv.org/html/2608.11237#bib.bib15); You et al\.,[2025](https://arxiv.org/html/2608.11237#bib.bib45)\), attention\-based interactions in place of fixed spectral mixing\(Li et al\.,[2023](https://arxiv.org/html/2608.11237#bib.bib23); Hao et al\.,[2023](https://arxiv.org/html/2608.11237#bib.bib13)\), and Koopman\-linearized Fourier operators\(Xiong et al\.,[2024](https://arxiv.org/html/2608.11237#bib.bib43)\)\. The second route compresses PDE dynamics into a learned latent space following the encode–evolve–decode paradigm\(Morton et al\.,[2018](https://arxiv.org/html/2608.11237#bib.bib29); Lusch et al\.,[2018](https://arxiv.org/html/2608.11237#bib.bib27); Brunton et al\.,[2022](https://arxiv.org/html/2608.11237#bib.bib5)\), realized within neural operators by LSM\(Wu et al\.,[2023](https://arxiv.org/html/2608.11237#bib.bib41)\), LNO\(Wang and Wang,[2024](https://arxiv.org/html/2608.11237#bib.bib37)\), and PI\-Latent\-NO\(Karumuri et al\.,[2026](https://arxiv.org/html/2608.11237#bib.bib18)\)\. Both routes improve prediction by enhancing backbone expressiveness, broadening spectral coverage, or learning compact latent representations\. However, they mostly model the state\-to\-state mapping or latent\-state evolution, while the latent increment that drives temporal evolution is rarely modeled as a structured object in its own right\.

## Appendix BMotivation

### B\.1\.Mean–Fluctuation Separation

The physical output contains components with different temporal properties\. As shown in Figure[5](https://arxiv.org/html/2608.11237#A2.F5), the temporal meanμ​\(u\)\\mu\(u\)captures the stable background or slowly varying structure, whereas the zero\-mean fluctuationr​\(u\)r\(u\)represents non\-stationary oscillations and phase\-related dynamics\. This decomposition suggests that mean errors and fluctuation errors should not be treated identically: mean errors correspond to mean drift, while fluctuation errors are more related to oscillatory mismatch and phase misalignment\. This motivates a structured reconstruction strategy that decouples mean and fluctuation components and restricts phase correction to the zero\-mean fluctuation field\.

![Refer to caption](https://arxiv.org/html/2608.11237v1/Figure/Figure4.png)Figure 5\.Mean–fluctuation decomposition of a physical time series, separating the stable mean component from the zero\-mean fluctuation dynamics\.
### B\.2\.Frequency\-Dependent Band Geometry

Beyond the non\-uniform spectral support, different active frequency bands also exhibit distinct latent\-channel geometries\. As shown in Figure[6](https://arxiv.org/html/2608.11237#A2.F6), the left panel visualizes the latent channel correlation matrices within representative low\-, middle\-, and high\-frequency active bands\. Different bands present clearly different channel\-coupling patterns, indicating that the geometric organization of latent increments varies across frequency regions\. The right panel further reports the effective rank and top eigen energy of different active bands\. These statistics differ noticeably across bands, suggesting that the complexity and concentration of channel\-wise variation are frequency dependent\. A single global geometric constraint is therefore unlikely to characterize all frequency regions adequately\. This supports a band\-wise design in which each active band uses an independent low\-rank geometric projector, allowing incremental geometry to match local spectral and channel structures\.

![Refer to caption](https://arxiv.org/html/2608.11237v1/Figure/Figure5.png)Figure 6\.Frequency\-dependent channel geometry of latent increments\. Distinct channel\-correlation matrices and band\-dependent rank statistics indicate that a single global geometry is insufficient for all active bands\.

## Appendix CMethod

### C\.1\.Theoretical Analysis

This section provides a conditional theoretical interpretation of GeoIncNO\. The goal is not to prove that the learned model always outperforms existing neural operators, but to clarify how increment\-centered modeling and active\-band geometric shaping can reduce physical prediction error under explicit structural assumptions\.

#### C\.1\.1\.Preliminaries and Notation

LetΩ⊂ℝds\\Omega\\subset\\mathbb\{R\}^\{d\_\{s\}\}denote the spatial domain, wheredsd\_\{s\}is the spatial dimension\. In this work,ds=1,2,3d\_\{s\}=1,2,3corresponds to the one\-, two\-, and three\-dimensional PDE benchmarks considered in the experiments\. LetIinI\_\{\\mathrm\{in\}\}andIoutI\_\{\\mathrm\{out\}\}denote the input and output time windows, respectively\. We consider input and output physical fieldsuin∈𝒳in:=L2​\(Iin×Ω;ℝCin\)u\_\{\\mathrm\{in\}\}\\in\\mathcal\{X\}\_\{\\mathrm\{in\}\}:=L^\{2\}\(I\_\{\\mathrm\{in\}\}\\times\\Omega;\\mathbb\{R\}^\{C\_\{\\mathrm\{in\}\}\}\)anduout∈𝒳out:=L2​\(Iout×Ω;ℝCout\)u\_\{\\mathrm\{out\}\}\\in\\mathcal\{X\}\_\{\\mathrm\{out\}\}:=L^\{2\}\(I\_\{\\mathrm\{out\}\}\\times\\Omega;\\mathbb\{R\}^\{C\_\{\\mathrm\{out\}\}\}\), whereCinC\_\{\\mathrm\{in\}\}andCoutC\_\{\\mathrm\{out\}\}are the numbers of input and output physical variables\. Unless otherwise specified,∥⋅∥\\\|\\cdot\\\|denotes theL2L^\{2\}norm in the corresponding function space, and its discrete Euclidean or Frobenius counterpart on finite grids\.

Let𝒢τ:𝒜⊂𝒳in→𝒳out\\mathcal\{G\}\_\{\\tau\}:\\mathcal\{A\}\\subset\\mathcal\{X\}\_\{\\mathrm\{in\}\}\\rightarrow\\mathcal\{X\}\_\{\\mathrm\{out\}\}denote the ground\-truth PDE solution operator over a prediction spanτ\>0\\tau\>0, where𝒜\\mathcal\{A\}is the admissible set of input trajectories\. The target output isuout⋆=𝒢τ​\(uin\)u\_\{\\mathrm\{out\}\}^\{\\star\}=\\mathcal\{G\}\_\{\\tau\}\(u\_\{\\mathrm\{in\}\}\)\.

To compare input and output fields in a common latent space, we introduce two compatible encoding mapsEin:𝒳in→𝒵E\_\{\\mathrm\{in\}\}:\\mathcal\{X\}\_\{\\mathrm\{in\}\}\\rightarrow\\mathcal\{Z\}andEout:𝒳out→𝒵E\_\{\\mathrm\{out\}\}:\\mathcal\{X\}\_\{\\mathrm\{out\}\}\\rightarrow\\mathcal\{Z\}, where𝒵\\mathcal\{Z\}is a latent function space\. Here,EinE\_\{\\mathrm\{in\}\}corresponds to the encoderEϕE\_\{\\phi\}used in the main architecture, whileEoutE\_\{\\mathrm\{out\}\}is introduced only for analysis to embed the ground\-truth output into the same latent space\. The current latent state and the ground\-truth next latent state are defined aszt=Ein​\(uin\)z\_\{t\}=E\_\{\\mathrm\{in\}\}\(u\_\{\\mathrm\{in\}\}\)andzt\+1⋆=Eout​\(𝒢τ​\(uin\)\)z\_\{t\+1\}^\{\\star\}=E\_\{\\mathrm\{out\}\}\(\\mathcal\{G\}\_\{\\tau\}\(u\_\{\\mathrm\{in\}\}\)\)\. The ground\-truth latent increment is then

\(25\)Δτ⋆​\(uin\):=zt\+1⋆−zt=Eout​\(𝒢τ​\(uin\)\)−Ein​\(uin\)\.\\Delta\_\{\\tau\}^\{\\star\}\(u\_\{\\mathrm\{in\}\}\):=z\_\{t\+1\}^\{\\star\}\-z\_\{t\}=E\_\{\\mathrm\{out\}\}\(\\mathcal\{G\}\_\{\\tau\}\(u\_\{\\mathrm\{in\}\}\)\)\-E\_\{\\mathrm\{in\}\}\(u\_\{\\mathrm\{in\}\}\)\.
GeoIncNO learns a latent increment operatorΔτpred\\Delta\_\{\\tau\}^\{\\mathrm\{pred\}\}and predicts the next latent state by residual advancement,zt\+1pred=zt\+Δτpred​\(zt\)z\_\{t\+1\}^\{\\mathrm\{pred\}\}=z\_\{t\}\+\\Delta\_\{\\tau\}^\{\\mathrm\{pred\}\}\(z\_\{t\}\)\. The corresponding physical\-space prediction is written as𝒢^τ​\(uin\)=D​\(zt\+1pred\)\\widehat\{\\mathcal\{G\}\}\_\{\\tau\}\(u\_\{\\mathrm\{in\}\}\)=D\(z\_\{t\+1\}^\{\\mathrm\{pred\}\}\), whereD:𝒵→𝒳outD:\\mathcal\{Z\}\\rightarrow\\mathcal\{X\}\_\{\\mathrm\{out\}\}denotes the overall decoding map\. In the actual GeoIncNO architecture,DDis implemented by the increment decoder, the state decoder, and the mean\-fluctuation reconstruction module; in the theoretical analysis, we treat them as a unified decoder when deriving error bounds\.

For a time\-dependent fieldv∈L2​\(I×Ω;ℝC\)v\\in L^\{2\}\(I\\times\\Omega;\\mathbb\{R\}^\{C\}\), we define the temporal mean and zero\-mean fluctuation operators asμt​\(v\)​\(x\):=1\|I\|​∫Iv​\(t,x\)​𝑑t\\mu\_\{t\}\(v\)\(x\):=\\frac\{1\}\{\|I\|\}\\int\_\{I\}v\(t,x\)\\,dtandcentert⁡\(v\)​\(t,x\):=v​\(t,x\)−μt​\(v\)​\(x\)\\operatorname\{center\}\_\{t\}\(v\)\(t,x\):=v\(t,x\)\-\\mu\_\{t\}\(v\)\(x\)\. By construction,μt​\(centert⁡\(v\)\)=0\\mu\_\{t\}\(\\operatorname\{center\}\_\{t\}\(v\)\)=0\. This continuous definition corresponds to the discrete temporal meanμ​\(⋅\)\\mu\(\\cdot\)and centering operation used in the main method\. We useℱ𝐱\\mathcal\{F\}\_\{\\mathbf\{x\}\}to denote the spatial spectral transform over alldsd\_\{s\}spatial dimensions, with frequency variableξ∈ℝds\\xi\\in\\mathbb\{R\}^\{d\_\{s\}\}in the continuous setting or𝐤∈ℤds\\mathbf\{k\}\\in\\mathbb\{Z\}^\{d\_\{s\}\}on a discrete grid\. We useℱt\\mathcal\{F\}\_\{t\}to denote the temporal Fourier transform, with temporal frequency variableω∈ℝ\\omega\\in\\mathbb\{R\}\. A hat over a variable denotes its spectral representation unless it is explicitly used for a model prediction, such asu^\\hat\{u\}\.

#### C\.1\.2\.Increment\-Centered Latent Evolution

A standard latent\-space neural operator maps the current latent state directly to the next latent state, i\.e\.,uin→Einzt→𝒦θzt\+1→𝐷u^u\_\{\\mathrm\{in\}\}\\xrightarrow\{E\_\{\\mathrm\{in\}\}\}z\_\{t\}\\xrightarrow\{\\mathcal\{K\}\_\{\\theta\}\}z\_\{t\+1\}\\xrightarrow\{D\}\\hat\{u\}, where𝒦θ:𝒵→𝒵\\mathcal\{K\}\_\{\\theta\}:\\mathcal\{Z\}\\rightarrow\\mathcal\{Z\}\. Such a formulation focuses on learning a state\-to\-state transition\. GeoIncNO instead adopts an increment\-centered formulation,zt\+1pred=zt\+Δτpred​\(zt\)z\_\{t\+1\}^\{\\mathrm\{pred\}\}=z\_\{t\}\+\\Delta\_\{\\tau\}^\{\\mathrm\{pred\}\}\(z\_\{t\}\), which separates the current latent state from the direction that drives its evolution\.

The increment formulation is consistent with the integral form of an evolutionary PDE\. Suppose that a physical trajectory satisfies∂tu​\(t\)=𝒩​\(u​\(t\)\)\\partial\_\{t\}u\(t\)=\\mathcal\{N\}\(u\(t\)\)\. Then the exact time advancement over a spanτ\\tausatisfies

\(26\)u​\(t\+τ\)−u​\(t\)=∫tt\+τ𝒩​\(u​\(s\)\)​𝑑s=τ​𝒩¯\[t,t\+τ\],u\(t\+\\tau\)\-u\(t\)=\\int\_\{t\}^\{t\+\\tau\}\\mathcal\{N\}\(u\(s\)\)\\,ds=\\tau\\overline\{\\mathcal\{N\}\}\_\{\[t,t\+\\tau\]\},where𝒩¯\[t,t\+τ\]:=1τ​∫tt\+τ𝒩​\(u​\(s\)\)​𝑑s\\overline\{\\mathcal\{N\}\}\_\{\[t,t\+\\tau\]\}:=\\frac\{1\}\{\\tau\}\\int\_\{t\}^\{t\+\\tau\}\\mathcal\{N\}\(u\(s\)\)\\,dsis the average evolution direction over\[t,t\+τ\]\[t,t\+\\tau\]\. Thus, learning the future state can be equivalently viewed as learning the state increment induced by the underlying dynamics\.

The same interpretation applies in the latent space\. If the latent trajectoryz​\(t\)z\(t\)is differentiable, then

\(27\)Δ​z⋆​\(t,τ\):=z​\(t\+τ\)−z​\(t\)=∫tt\+τ∂sz​\(s\)​d​s\.\\Delta z^\{\\star\}\(t,\\tau\):=z\(t\+\\tau\)\-z\(t\)=\\int\_\{t\}^\{t\+\\tau\}\\partial\_\{s\}z\(s\)\\,ds\.If∂tz\\partial\_\{t\}zis locally Lipschitz in time, thenΔ​z⋆​\(t,τ\)=τ​∂tz​\(t\)\+O​\(τ2\)\\Delta z^\{\\star\}\(t,\\tau\)=\\tau\\partial\_\{t\}z\(t\)\+O\(\\tau^\{2\}\), showing that the latent increment represents the average local evolution direction up to higher\-order terms\.

GeoIncNO parameterizes this latent increment through a raw increment prediction followed by active\-band geometric correction:

\(28\)δ​zraw\\displaystyle\\delta z\_\{\\mathrm\{raw\}\}=Gψ​\(𝒯θ​\(zt\)\),\\displaystyle=G\_\{\\psi\}\(\\mathcal\{T\}\_\{\\theta\}\(z\_\{t\}\)\),\(29\)Δτpred​\(zt\)\\displaystyle\\Delta\_\{\\tau\}^\{\\mathrm\{pred\}\}\(z\_\{t\}\)=δ​zraw\+η​\(PAB​\(δ​zraw\)−δ​zraw\),\\displaystyle=\\delta z\_\{\\mathrm\{raw\}\}\+\\eta\\left\(P\_\{\\mathrm\{AB\}\}\(\\delta z\_\{\\mathrm\{raw\}\}\)\-\\delta z\_\{\\mathrm\{raw\}\}\\right\),\(30\)zt\+1pred\\displaystyle z\_\{t\+1\}^\{\\mathrm\{pred\}\}=zt\+Δτpred​\(zt\)\.\\displaystyle=z\_\{t\}\+\\Delta\_\{\\tau\}^\{\\mathrm\{pred\}\}\(z\_\{t\}\)\.In this form,Gψ​\(𝒯θ​\(⋅\)\)G\_\{\\psi\}\(\\mathcal\{T\}\_\{\\theta\}\(\\cdot\)\)provides an unconstrained transition direction, whilePABP\_\{\\mathrm\{AB\}\}supplies a structured active\-band correction\. Therefore, the geometric constraint acts on the latent increment that drives state evolution, rather than on the full latent state itself\.

#### C\.1\.3\.Common Assumptions

The following assumptions are used to provide conditional theoretical interpretations of GeoIncNO\. They are not intended to claim that any trained instance of GeoIncNO is guaranteed to outperform existing neural operators\. Instead, they specify the structural conditions under which the proposed increment\-centered and active\-band designs can be related to prediction error control\.

A0\. Function spaces and spectral representation\.We assume that the spatial domainΩ⊂ℝds\\Omega\\subset\\mathbb\{R\}^\{d\_\{s\}\}is either periodic or admits an orthogonal spectral basis compatible with the boundary conditions\. For periodic domains,ℱ𝐱\\mathcal\{F\}\_\{\\mathbf\{x\}\}is thedsd\_\{s\}\-dimensional Fourier transform\. For non\-periodic domains, the Fourier basis can be replaced by a Laplacian eigenbasis, a Chebyshev basis, or another orthogonal basis adapted to the boundary conditions\. All spectral decompositions, orthogonality arguments, and applications of Parseval’s identity below are understood with respect to the chosen basis\.

A1\. Latent reconstruction and decoder regularity\.We assume that the latent representation preserves the information of the ground\-truth solution trajectories up to a bounded reconstruction error\. Specifically, for the output\-side encoderEoutE\_\{\\mathrm\{out\}\}and the unified decoderDD, there existsεlat≥0\\varepsilon\_\{\\mathrm\{lat\}\}\\geq 0such that

\(31\)supuout∈𝒢τ​\(𝒜\)‖D​\(Eout​\(uout\)\)−uout‖𝒳out≤εlat\.\\sup\_\{u\_\{\\mathrm\{out\}\}\\in\\mathcal\{G\}\_\{\\tau\}\(\\mathcal\{A\}\)\}\\left\\\|D\(E\_\{\\mathrm\{out\}\}\(u\_\{\\mathrm\{out\}\}\)\)\-u\_\{\\mathrm\{out\}\}\\right\\\|\_\{\\mathcal\{X\}\_\{\\mathrm\{out\}\}\}\\leq\\varepsilon\_\{\\mathrm\{lat\}\}\.We also assume thatDDis Lipschitz continuous on the relevant latent trajectory set; that is, there existsLD\>0L\_\{D\}\>0such that

\(32\)‖D​\(z1\)−D​\(z2\)‖𝒳out≤LD​‖z1−z2‖𝒵\.\\left\\\|D\(z\_\{1\}\)\-D\(z\_\{2\}\)\\right\\\|\_\{\\mathcal\{X\}\_\{\\mathrm\{out\}\}\}\\leq L\_\{D\}\\left\\\|z\_\{1\}\-z\_\{2\}\\right\\\|\_\{\\mathcal\{Z\}\}\.This assumption allows latent increment errors to be transferred to physical\-space prediction errors\.

A2\. Active\-band calibration\.In GeoIncNO, active bands are calibrated from the spectral energy distribution of the model\-predicted raw incrementδ​zraw\\delta z\_\{\\mathrm\{raw\}\}, whereas the theoretical analysis concerns the spectral structure of the ground\-truth latent incrementΔτ⋆\\Delta\_\{\\tau\}^\{\\star\}\. We assume that the calibrated effective spectral supportΩeff\\Omega\_\{\\mathrm\{eff\}\}covers the dominant spectral support ofΔτ⋆\\Delta\_\{\\tau\}^\{\\star\}up to a calibration errorεcal\\varepsilon\_\{\\mathrm\{cal\}\}:

\(33\)‖ℱ𝐱​\(Δτ⋆\)‖Ωeffc2≤εcal2\.\\left\\\|\\mathcal\{F\}\_\{\\mathbf\{x\}\}\(\\Delta\_\{\\tau\}^\{\\star\}\)\\right\\\|\_\{\\Omega\_\{\\mathrm\{eff\}\}^\{c\}\}^\{2\}\\leq\\varepsilon\_\{\\mathrm\{cal\}\}^\{2\}\.Here,ℱ𝐱​\(Δτ⋆\)\\mathcal\{F\}\_\{\\mathbf\{x\}\}\(\\Delta\_\{\\tau\}^\{\\star\}\)denotes the spectral representation of the ground\-truth latent increment, and∥⋅∥Ωeffc\\\|\\cdot\\\|\_\{\\Omega\_\{\\mathrm\{eff\}\}^\{c\}\}denotes the spectralL2L^\{2\}norm restricted to the complement ofΩeff\\Omega\_\{\\mathrm\{eff\}\}\.

A3\. Approximate low\-rank structure within active bands\.Let\{ℬm\}m=1M\\\{\\mathcal\{B\}\_\{m\}\\\}\_\{m=1\}^\{M\}be the active bands, which are mutually disjoint and coverΩeff\\Omega\_\{\\mathrm\{eff\}\}\. For each bandℬm\\mathcal\{B\}\_\{m\}, we assume that there exists a low\-dimensional subspace𝒰m⊂ℝd\\mathcal\{U\}\_\{m\}\\subset\\mathbb\{R\}^\{d\}with orthogonal projectorQmQ\_\{m\}such that the ground\-truth latent increment has small residual energy outside these subspaces:

\(34\)∑m=1M‖\(I−Qm\)​ℱ𝐱​\(Δτ,m⋆\)‖2≤εband2\.\\sum\_\{m=1\}^\{M\}\\left\\\|\(I\-Q\_\{m\}\)\\mathcal\{F\}\_\{\\mathbf\{x\}\}\(\\Delta\_\{\\tau,m\}^\{\\star\}\)\\right\\\|^\{2\}\\leq\\varepsilon\_\{\\mathrm\{band\}\}^\{2\}\.Here,Δτ,m⋆\\Delta\_\{\\tau,m\}^\{\\star\}denotes the restriction ofΔτ⋆\\Delta\_\{\\tau\}^\{\\star\}to themm\-th active band in the spectral domain\. This assumption specifies the spectral condition under which low\-rank active\-band shaping has a theoretical interpretation\.

A4\. Controlled out\-of\-subspace and out\-of\-band prediction energy\.Approximate low\-rank structure of the ground\-truth increment alone is insufficient to control the model error, since the learned increment may introduce additional energy outside the active low\-rank subspaces or outside the effective spectral support\. We therefore assume that the predicted latent incrementΔτpred\\Delta\_\{\\tau\}^\{\\mathrm\{pred\}\}satisfies

\(35\)∑m=1M‖\(I−Qm\)​ℱ𝐱​\(Δτ,mpred\)‖2\\displaystyle\\sum\_\{m=1\}^\{M\}\\left\\\|\(I\-Q\_\{m\}\)\\mathcal\{F\}\_\{\\mathbf\{x\}\}\(\\Delta\_\{\\tau,m\}^\{\\mathrm\{pred\}\}\)\\right\\\|^\{2\}≤εpred,out2,\\displaystyle\\leq\\varepsilon\_\{\\mathrm\{pred,out\}\}^\{2\},‖ℱ𝐱​\(Δτpred\)‖Ωeffc2\\displaystyle\\left\\\|\\mathcal\{F\}\_\{\\mathbf\{x\}\}\(\\Delta\_\{\\tau\}^\{\\mathrm\{pred\}\}\)\\right\\\|\_\{\\Omega\_\{\\mathrm\{eff\}\}^\{c\}\}^\{2\}≤εoutband2\.\\displaystyle\\leq\\varepsilon\_\{\\mathrm\{outband\}\}^\{2\}\.Here,Δτ,mpred\\Delta\_\{\\tau,m\}^\{\\mathrm\{pred\}\}denotes the restriction of the predicted latent increment toℬm\\mathcal\{B\}\_\{m\}in the spectral domain\. The quantitiesεpred,out\\varepsilon\_\{\\mathrm\{pred,out\}\}andεoutband\\varepsilon\_\{\\mathrm\{outband\}\}control the prediction energy outside the active low\-rank subspaces and outside the effective spectral support, respectively\. This assumption corresponds to the role of the increment geometry regularizer and the projection consistency loss\.

#### C\.1\.4\.Latent Increment Error Controls Physical Prediction Error

We first show that the physical\-space prediction error can be controlled by the latent increment error, provided that the decoder is Lipschitz continuous and the latent representation has bounded output\-side reconstruction error\. This result justifies why GeoIncNO focuses on organizing the latent increment\.

###### Theorem C\.1 \(Latent increment error control\)\.

Under Assumption A1, for anyuin∈𝒜u\_\{\\mathrm\{in\}\}\\in\\mathcal\{A\}, the one\-step physical prediction error of GeoIncNO satisfies

\(36\)‖𝒢^τ​\(uin\)−𝒢τ​\(uin\)‖𝒳out≤εlat\+LD​‖Δτpred​\(zt\)−Δτ⋆​\(uin\)‖𝒵,\\left\\\|\\widehat\{\\mathcal\{G\}\}\_\{\\tau\}\(u\_\{\\mathrm\{in\}\}\)\-\\mathcal\{G\}\_\{\\tau\}\(u\_\{\\mathrm\{in\}\}\)\\right\\\|\_\{\\mathcal\{X\}\_\{\\mathrm\{out\}\}\}\\leq\\varepsilon\_\{\\mathrm\{lat\}\}\+L\_\{D\}\\left\\\|\\Delta\_\{\\tau\}^\{\\mathrm\{pred\}\}\(z\_\{t\}\)\-\\Delta\_\{\\tau\}^\{\\star\}\(u\_\{\\mathrm\{in\}\}\)\\right\\\|\_\{\\mathcal\{Z\}\},wherezt=Ein​\(uin\)z\_\{t\}=E\_\{\\mathrm\{in\}\}\(u\_\{\\mathrm\{in\}\}\),Δτpred​\(zt\)\\Delta\_\{\\tau\}^\{\\mathrm\{pred\}\}\(z\_\{t\}\)is the predicted latent increment, andΔτ⋆​\(uin\)\\Delta\_\{\\tau\}^\{\\star\}\(u\_\{\\mathrm\{in\}\}\)is the ground\-truth latent increment\.

###### Proof\.

Letzt=Ein​\(uin\)z\_\{t\}=E\_\{\\mathrm\{in\}\}\(u\_\{\\mathrm\{in\}\}\)andzt\+1⋆=Eout​\(𝒢τ​\(uin\)\)z\_\{t\+1\}^\{\\star\}=E\_\{\\mathrm\{out\}\}\(\\mathcal\{G\}\_\{\\tau\}\(u\_\{\\mathrm\{in\}\}\)\)\. The predicted and ground\-truth next latent states satisfyzt\+1pred=zt\+Δτpred​\(zt\)z\_\{t\+1\}^\{\\mathrm\{pred\}\}=z\_\{t\}\+\\Delta\_\{\\tau\}^\{\\mathrm\{pred\}\}\(z\_\{t\}\)andzt\+1⋆=zt\+Δτ⋆​\(uin\)z\_\{t\+1\}^\{\\star\}=z\_\{t\}\+\\Delta\_\{\\tau\}^\{\\star\}\(u\_\{\\mathrm\{in\}\}\), respectively\. Using the unified decoderDD, the model prediction is𝒢^τ​\(uin\)=D​\(zt\+1pred\)\\widehat\{\\mathcal\{G\}\}\_\{\\tau\}\(u\_\{\\mathrm\{in\}\}\)=D\(z\_\{t\+1\}^\{\\mathrm\{pred\}\}\)\.

By adding and subtractingD​\(zt\+1⋆\)D\(z\_\{t\+1\}^\{\\star\}\), we obtain

‖𝒢^τ​\(uin\)−𝒢τ​\(uin\)‖𝒳out\\displaystyle\\left\\\|\\widehat\{\\mathcal\{G\}\}\_\{\\tau\}\(u\_\{\\mathrm\{in\}\}\)\-\\mathcal\{G\}\_\{\\tau\}\(u\_\{\\mathrm\{in\}\}\)\\right\\\|\_\{\\mathcal\{X\}\_\{\\mathrm\{out\}\}\}\(37\)≤‖D​\(zt\+1pred\)−D​\(zt\+1⋆\)‖𝒳out\+‖D​\(zt\+1⋆\)−𝒢τ​\(uin\)‖𝒳out\.\\displaystyle\\leq\\left\\\|D\(z\_\{t\+1\}^\{\\mathrm\{pred\}\}\)\-D\(z\_\{t\+1\}^\{\\star\}\)\\right\\\|\_\{\\mathcal\{X\}\_\{\\mathrm\{out\}\}\}\+\\left\\\|D\(z\_\{t\+1\}^\{\\star\}\)\-\\mathcal\{G\}\_\{\\tau\}\(u\_\{\\mathrm\{in\}\}\)\\right\\\|\_\{\\mathcal\{X\}\_\{\\mathrm\{out\}\}\}\.By the Lipschitz continuity ofDD, the first term is bounded byLD​‖zt\+1pred−zt\+1⋆‖𝒵L\_\{D\}\\\|z\_\{t\+1\}^\{\\mathrm\{pred\}\}\-z\_\{t\+1\}^\{\\star\}\\\|\_\{\\mathcal\{Z\}\}\. Sincezt\+1pred−zt\+1⋆=Δτpred​\(zt\)−Δτ⋆​\(uin\)z\_\{t\+1\}^\{\\mathrm\{pred\}\}\-z\_\{t\+1\}^\{\\star\}=\\Delta\_\{\\tau\}^\{\\mathrm\{pred\}\}\(z\_\{t\}\)\-\\Delta\_\{\\tau\}^\{\\star\}\(u\_\{\\mathrm\{in\}\}\), we have

\(38\)‖D​\(zt\+1pred\)−D​\(zt\+1⋆\)‖𝒳out≤LD​‖Δτpred​\(zt\)−Δτ⋆​\(uin\)‖𝒵\.\\left\\\|D\(z\_\{t\+1\}^\{\\mathrm\{pred\}\}\)\-D\(z\_\{t\+1\}^\{\\star\}\)\\right\\\|\_\{\\mathcal\{X\}\_\{\\mathrm\{out\}\}\}\\leq L\_\{D\}\\left\\\|\\Delta\_\{\\tau\}^\{\\mathrm\{pred\}\}\(z\_\{t\}\)\-\\Delta\_\{\\tau\}^\{\\star\}\(u\_\{\\mathrm\{in\}\}\)\\right\\\|\_\{\\mathcal\{Z\}\}\.
The second term is bounded by the output\-side reconstruction error in Assumption A1:

\(39\)‖D​\(zt\+1⋆\)−𝒢τ​\(uin\)‖𝒳out\\displaystyle\\left\\\|D\(z\_\{t\+1\}^\{\\star\}\)\-\\mathcal\{G\}\_\{\\tau\}\(u\_\{\\mathrm\{in\}\}\)\\right\\\|\_\{\\mathcal\{X\}\_\{\\mathrm\{out\}\}\}=‖D​\(Eout​\(𝒢τ​\(uin\)\)\)−𝒢τ​\(uin\)‖𝒳out\\displaystyle=\\left\\\|D\\\!\\left\(E\_\{\\mathrm\{out\}\}\(\\mathcal\{G\}\_\{\\tau\}\(u\_\{\\mathrm\{in\}\}\)\)\\right\)\-\\mathcal\{G\}\_\{\\tau\}\(u\_\{\\mathrm\{in\}\}\)\\right\\\|\_\{\\mathcal\{X\}\_\{\\mathrm\{out\}\}\}≤εlat\.\\displaystyle\\leq\\varepsilon\_\{\\mathrm\{lat\}\}\.Combining the two bounds proves the theorem\. ∎

Algorithm[1](https://arxiv.org/html/2608.11237#alg1)summarizes the complete inference procedure of GeoIncNO\. Given the input physical fields, the model predicts and geometrically refines a latent increment, advances the latent state, and reconstructs the future fields through dual\-branch MFDR with fluctuation\-only phase correction\.

Algorithm 1The inference processing of GeoIncNO1:Input physical fields

uinu\_\{\\mathrm\{in\}\}, ABP operator

PABP\_\{\\mathrm\{AB\}\}, correction strength

η\\eta
2:Predicted future physical fields

u^\\hat\{u\}
3:Step 1: Latent Increment Prediction

4:Encode the input field

zin=Eϕ​\(uin\)z\_\{\\mathrm\{in\}\}=E\_\{\\phi\}\(u\_\{\\mathrm\{in\}\}\)and set

zt=zinz\_\{t\}=z\_\{\\mathrm\{in\}\}
5:Extract transition features

ht=𝒯θ​\(zt\)h\_\{t\}=\\mathcal\{T\}\_\{\\theta\}\(z\_\{t\}\)
6:Predict the raw latent increment

δ​zraw=Gψ​\(ht\)\\delta z\_\{\\mathrm\{raw\}\}=G\_\{\\psi\}\(h\_\{t\}\)
7:Step 2: ABP Correction

8:Apply ABP,

δ​zproj=PAB​\(δ​zraw\)\\delta z\_\{\\mathrm\{proj\}\}=P\_\{\\mathrm\{AB\}\}\(\\delta z\_\{\\mathrm\{raw\}\}\)
9:Compute the refined increment

δ​z=δ​zraw\+η​\(δ​zproj−δ​zraw\)\\delta z=\\delta z\_\{\\mathrm\{raw\}\}\+\\eta\(\\delta z\_\{\\mathrm\{proj\}\}\-\\delta z\_\{\\mathrm\{raw\}\}\)
10:Update the latent state

zt\+1=zt\+δ​zz\_\{t\+1\}=z\_\{t\}\+\\delta z
11:Step 3: Dual\-Branch Reconstruction

12:Compute the input\-anchored prediction

ucorr=α​Aphys​\(uin\)\+DΔ​\(δ​z\)u\_\{\\mathrm\{corr\}\}=\\alpha A\_\{\\mathrm\{phys\}\}\(u\_\{\\mathrm\{in\}\}\)\+D\_\{\\Delta\}\(\\delta z\)
13:Compute the state\-decoded prediction

ustate=Dstate​\(zt\+1\)u\_\{\\mathrm\{state\}\}=D\_\{\\mathrm\{state\}\}\(z\_\{t\+1\}\)
14:Step 4: MFDR Fusion

15:Decompose

ucorru\_\{\\mathrm\{corr\}\}and

ustateu\_\{\\mathrm\{state\}\}into mean and fluctuation components

16:Fuse the mean components with

Γμ\\Gamma\_\{\\mu\}to obtain

μ^\\hat\{\\mu\}
17:Fuse and center the fluctuation components with

Γr\\Gamma\_\{r\}to obtain

r^\\hat\{r\}
18:Step 5: Fluctuation\-Only Phase Correction

19:Apply phase correction to the fused fluctuation and obtain

rphaser\_\{\\mathrm\{phase\}\}
20:Reconstruct the final prediction

u^=repeatt⁡\(μ^\)\+r^\+rphase\\hat\{u\}=\\operatorname\{repeat\}\_\{t\}\(\\hat\{\\mu\}\)\+\\hat\{r\}\+r\_\{\\mathrm\{phase\}\}
21:return

u^\\hat\{u\}

## Appendix DExperiments

### D\.1\.Experiments Setup

#### D\.1\.1\.PDE Benchmarks

Table 3\.Summary of PDE benchmarks\.Table 4\.Prediction accuracy and parameter count across six PDE benchmarks\. Results are reported as mean±\\pmstandard deviation over five independent runs; lower values indicate better accuracy\.We evaluate GeoIncNO on six PDE benchmarks generated following the PDEBench\-style data generation protocol: Burgers, KS, NS, SW equation, 3D compressible Euler \(CE\), and Maxwell\. These benchmarks cover 1D, 2D, and 3D dynamical systems, including convection–diffusion dynamics, spatiotemporal chaos, incompressible vortical flow, free\-surface wave propagation, compressible flow, and electromagnetic wave propagation\.

- •*Burgers*evaluates nonlinear transport and viscous dissipation, and is governed by∂tu\+u​∂xu=ν​∂x​xu,x∈\[0,2​π\)\\partial\_\{t\}u\+u\\partial\_\{x\}u=\\nu\\partial\_\{xx\}u,\\quad x\\in\[0,2\\pi\), whereuuis the scalar velocity field andν\\nuis the viscosity\.
- •*KS*evaluates long\-term prediction under spatiotemporal chaos, and is governed by∂tu\+u​∂xu\+∂x​xu\+∂x​x​x​xu=0,x∈\[0,64\]\\partial\_\{t\}u\+u\\partial\_\{x\}u\+\\partial\_\{xx\}u\+\\partial\_\{xxxx\}u=0,\\quad x\\in\[0,64\], whereuudenotes the scalar state variable\.
- •*NS*evaluates incompressible vortical dynamics, and is governed by the vorticity formulation∂tω\+\(𝐮⋅∇\)​ω=ν​Δ​ω\+f,∇⋅𝐮=0\\partial\_\{t\}\\omega\+\(\\mathbf\{u\}\\cdot\\nabla\)\\omega=\\nu\\Delta\\omega\+f,\\quad\\nabla\\cdot\\mathbf\{u\}=0, whereω\\omega,𝐮\\mathbf\{u\},ν\\nu, andffdenote vorticity, velocity, viscosity, and forcing, respectively\. We use the unforced settingf=0f=0\.
- •*SW*evaluates free\-surface wave propagation and nonlinear transport, and is governed by∂th\+∇⋅\(h​𝐮\)=0,∂t𝐮\+\(𝐮⋅∇\)​𝐮\+g​∇h=0\\partial\_\{t\}h\+\\nabla\\cdot\(h\\mathbf\{u\}\)=0,\\quad\\partial\_\{t\}\\mathbf\{u\}\+\(\\mathbf\{u\}\\cdot\\nabla\)\\mathbf\{u\}\+g\\nabla h=0, wherehhis the fluid height,𝐮=\(u,v\)\\mathbf\{u\}=\(u,v\)is the horizontal velocity field, andggis the gravitational acceleration\.
- •*3D CE*evaluates density, momentum, and energy evolution under compressible dynamics, and is governed by the compressible Euler system \(40\)∂tρ\+∇⋅\(ρ​𝐯\)\\displaystyle\\partial\_\{t\}\\rho\+\\nabla\\cdot\(\\rho\\mathbf\{v\}\)=0,\\displaystyle=0,∂t\(ρ​𝐯\)\+∇⋅\(ρ​𝐯⊗𝐯\+p​I\)\\displaystyle\\partial\_\{t\}\(\\rho\\mathbf\{v\}\)\+\\nabla\\cdot\(\\rho\\mathbf\{v\}\\otimes\\mathbf\{v\}\+pI\)=0,\\displaystyle=0,∂tE\+∇⋅\(\(E\+p\)​𝐯\)\\displaystyle\\partial\_\{t\}E\+\\nabla\\cdot\(\(E\+p\)\\mathbf\{v\}\)=0,\\displaystyle=0,withp=\(γ−1\)​\(E−12​ρ​\|𝐯\|2\)p=\(\\gamma\-1\)\(E\-\\frac\{1\}\{2\}\\rho\|\\mathbf\{v\}\|^\{2\}\), whereρ\\rho,𝐯\\mathbf\{v\},pp, andEEdenote density, velocity, pressure, and total energy, respectively\.
- •*Maxwell*evaluates electromagnetic wave propagation, and is governed in a source\-free vacuum by∂t𝐄=c2​∇×𝐁,∂t𝐁=−∇×𝐄\\partial\_\{t\}\\mathbf\{E\}=c^\{2\}\\nabla\\times\\mathbf\{B\},\\quad\\partial\_\{t\}\\mathbf\{B\}=\-\\nabla\\times\\mathbf\{E\}, with divergence\-free constraints∇⋅𝐄=0\\nabla\\cdot\\mathbf\{E\}=0and∇⋅𝐁=0\\nabla\\cdot\\mathbf\{B\}=0, where𝐄\\mathbf\{E\},𝐁\\mathbf\{B\}, andccdenote the electric field, magnetic field, and wave speed, respectively\.

#### D\.1\.2\.Evaluation Metrics

We report all metrics averaged over the test set to evaluate complementary aspects of prediction quality\. Field\-level accuracy is measured by the relativeL2L\_\{2\}errorRel​\-​L2=‖u−u^‖2/‖u‖2\\mathrm\{Rel\}\\text\{\-\}L\_\{2\}=\\\|u\-\\hat\{u\}\\\|\_\{2\}/\\\|u\\\|\_\{2\}and MSE, where Rel\-L2L\_\{2\}provides scale\-normalized accuracy and MSE reflects absolute pointwise deviations\. To assess whether the predicted solution also preserves local variation, we report the relative Sobolev errorRel​\-​H1=‖u−u^‖H1/‖u‖H1\\mathrm\{Rel\}\\text\{\-\}H^\{1\}=\\\|u\-\\hat\{u\}\\\|\_\{H^\{1\}\}/\\\|u\\\|\_\{H^\{1\}\}, where theH1H^\{1\}norm includes both field and gradient errors\. This metric is particularly useful for derivative\-sensitive dynamics, such as the Kuramoto–Sivashinsky \(KS\) equation\. To evaluate spectral fidelity, we use the weighted log ratio \(WLR\)WLR=∑k∈𝒦wk​\|log⁡Epred​\(k\)\+ϵEref​\(k\)\+ϵ\|\\mathrm\{WLR\}=\\sum\_\{k\\in\\mathcal\{K\}\}w\_\{k\}\\left\|\\log\\frac\{E\_\{\\mathrm\{pred\}\}\(k\)\+\\epsilon\}\{E\_\{\\mathrm\{ref\}\}\(k\)\+\\epsilon\}\\right\|\.Epred​\(k\)=\|ℱ​\(u^\)​\(k\)\|2E\_\{\\mathrm\{pred\}\}\(k\)=\|\\mathcal\{F\}\(\\hat\{u\}\)\(k\)\|^\{2\}andEref​\(k\)=\|ℱ​\(u\)​\(k\)\|2E\_\{\\mathrm\{ref\}\}\(k\)=\|\\mathcal\{F\}\(u\)\(k\)\|^\{2\}denote the spectral energies of the predicted and reference fields at modekk, respectively\. The weightwkw\_\{k\}is computed from the normalized reference spectral energy, andϵ\\epsilonis used for numerical stability\. For time\-dependent rollouts, we additionally report two RMSE\-based stability metrics\. The final\-step RMSE \(F\-RMSE\) measures the prediction error at the last rollout step, whereas the cumulative RMSE \(C\-RMSE\) measures the accumulated trajectory error over all rollout steps\. When boundary accuracy is explicitly evaluated, boundary RMSE \(B\-RMSE\) is further used to measure the rollout error restricted to boundary regions\. To quantify long\-term mean\-field preservation, we further report the relative temporal\-mean drift \(Mean Drift\),‖μ​\(u\)−μ​\(u^\)‖2/\(‖μ​\(u\)‖2\+ϵ\)\\\|\\mu\(u\)\-\\mu\(\\hat\{u\}\)\\\|\_\{2\}/\(\\\|\\mu\(u\)\\\|\_\{2\}\+\\epsilon\), whereμ​\(⋅\)=1T​∑τ=1T\(⋅\)τ\\mu\(\\cdot\)=\\tfrac\{1\}\{T\}\\sum\_\{\\tau=1\}^\{T\}\(\\cdot\)\_\{\\tau\}is the temporal mean over the rollout; it measures the slowly accumulated bias in the long\-term average field\.

In addition to RMSE, Rel\-L2L\_\{2\}, and F\-RMSE, we also report MAE, R2, and FE to provide complementary evaluations of prediction accuracy and rollout stability\. MAE measures the average absolute deviation between the predicted and reference fields,MAE=1N​∑i=1N\|yi−y^i\|\\mathrm\{MAE\}=\\frac\{1\}\{N\}\\sum\_\{i=1\}^\{N\}\|y\_\{i\}\-\\hat\{y\}\_\{i\}\|, where it provides an intuitive measure of pointwise prediction errors\. The coefficient of determination R2evaluates the agreement between the predicted and reference solutions,R2=1−∑i\(yi−y^i\)2∑i\(yi−y¯\)2R^\{2\}=1\-\\frac\{\\sum\_\{i\}\(y\_\{i\}\-\\hat\{y\}\_\{i\}\)^\{2\}\}\{\\sum\_\{i\}\(y\_\{i\}\-\\bar\{y\}\)^\{2\}\}, wherey¯\\bar\{y\}denotes the mean value of the reference solution\. A higher R2indicates that the prediction explains a larger proportion of the variance in the ground\-truth field\. For autoregressive evaluation, FE measures the accumulated forecasting error during multi\-step rollout and reflects the error propagation behavior over the prediction horizon\. Lower values indicate better performance for RMSE, MAE, Rel\-L2L\_\{2\}, F\-RMSE, and FE, while higher values indicate better performance for R2\.

#### D\.1\.3\.Implementation Details

All experiments are implemented in PyTorch and run on a single NVIDIA RTX 4090 GPU with 24 GB memory\. For each benchmark, we use 1,000 trajectories for training and 200 independently generated trajectories for testing\. To keep the optimization setting of baseline neural operators consistent, all non\-physics\-informed baseline models are trained with the MSE prediction loss only\. PINO is treated separately because it is a physics\-informed baseline, and is trained with the MSE prediction loss together with the PDE residual loss following its original formulation\. GeoIncNO is optimized with the objective defined in Section[4\.6](https://arxiv.org/html/2608.11237#S4.SS6), which augments the prediction loss with active\-band increment geometry regularization and projection consistency\. We use the Adam optimizer\(Kingma and Ba,[2015](https://arxiv.org/html/2608.11237#bib.bib19)\)with an initial learning rate of1×10−41\\times 10^\{\-4\}\. A cosine annealing learning\-rate scheduler is adopted to decay the learning rate during training gradually\. Each model is trained for 100 epochs\. To improve long\-horizon stability, all models are trained using the pushforward strategy with a rollout length ofT=5T=5, in which predictions are recursively fed back into the model during training\. For GeoIncNO, we useM=4M=4active frequency bands, set the low\-rank projector dimension tor=8r=8, and set the ABP correction strength toη=0\.05\\eta=0\.05unless otherwise specified\. Across benchmarks, only the dimension\-dependent operators differ: 1D/2D/3D tasks use Conv1d/2d/3d with matching FFTs, while temporal\-window processing is factored out through temporal projection, so no 4D convolution is needed and the ABP transform acts only along spatial axes\.

### D\.2\.Long\-Horizon Rollout Stability

Due to space constraints in the main text, we place the full long\-horizon rollout comparison here\. Figure[7](https://arxiv.org/html/2608.11237#A4.F7)reports the MSE, Rel\-L2L\_\{2\}, and Rel\-H1H^\{1\}of GeoIncNO and all baselines at five rollout checkpoints \(0\.2​T0\.2T–1\.0​T1\.0T\) across the six PDE benchmarks, complementing the analysis in Section[5\.3](https://arxiv.org/html/2608.11237#S5.SS3)\. Across every benchmark and metric, GeoIncNO attains the lowest error at each checkpoint and shows the slowest error growth over the horizon, confirming that explicitly structuring the latent transition increment suppresses long\-horizon error accumulation\.

![Refer to caption](https://arxiv.org/html/2608.11237v1/Figure/97c148eca814cbb1577de1057eb4d466.png)Figure 7\.Long\-horizon rollout errors across the six PDE benchmarks\. Each column is a benchmark and each row a metric \(MSE, Rel\-H1H^\{1\}, Rel\-L2L\_\{2\}\); markers report the error at rollout checkpoints0\.2​T0\.2T–1\.0​T1\.0T\. GeoIncNO \(red\) attains the lowest error at every checkpoint and accumulates error most slowly over the horizon\.
### D\.3\.Cross\-Parameter Generalization

To evaluate cross\-parameter generalization, we conduct experiments on the 2D NS benchmark at resolution256×256256\\times 256\. The training set is constructed by mixing trajectories generated at two Reynolds numbers,Re=1000\\mathrm\{Re\}=1000andRe=10000\\mathrm\{Re\}=10000, while testing is performed on the unseen intermediate settingRe=5000\\mathrm\{Re\}=5000\. This setting evaluates whether the model can generalize to an unseen physical\-parameter regime rather than only fit a fixed Reynolds\-number distribution\. For each training Reynolds number, we use 500 trajectories, and the test set contains 100 independently generated trajectories atRe=5000\\mathrm\{Re\}=5000\. As shown in Table[5](https://arxiv.org/html/2608.11237#A4.T5), GeoIncNO achieves the strongest overall generalization to the unseenRe=5000\\mathrm\{Re\}=5000setting, obtaining the best result on every reported metric\. Compared with UNO, the strongest competing baseline on the main rollout metrics, GeoIncNO reduces F\-RMSE from7\.17​E−037\.17\\mathrm\{E\}\{\-03\}to3\.94​E−033\.94\\mathrm\{E\}\{\-03\}, C\-RMSE from5\.99​E−025\.99\\mathrm\{E\}\{\-02\}to3\.30​E−023\.30\\mathrm\{E\}\{\-02\}, and WLR from8\.46​E−028\.46\\mathrm\{E\}\{\-02\}to4\.65​E−024\.65\\mathrm\{E\}\{\-02\}\. It also obtains the lowest Rel\-L2L\_\{2\}, Rel\-H1H^\{1\}, MSE, and B\-RMSE, indicating better field\-level, derivative\-level, and boundary\-region accuracy under Reynolds\-number interpolation\. These results suggest that GeoIncNO improves cross\-parameter robustness under Reynolds\-number interpolation\.

Table 5\.Cross\-parameter generalization on the 2D NS benchmark under Reynolds\-number interpolation \(trained onRe=1000/10000\\mathrm\{Re\}=1000/10000, tested on the unseenRe=5000\\mathrm\{Re\}=5000\)\. Lower is better for all metrics\.
### D\.4\.Long\-Horizon Rollout Stability under Cross\-Parameter Interpolation

To further evaluate rollout stability under an unseen physical\-parameter regime, we report prediction errors at different rollout checkpoints on the cross\-parameter 2D NS setting\. Following the cross\-parameter protocol in Appendix[D\.3](https://arxiv.org/html/2608.11237#A4.SS3), all models are trained on a mixed dataset at resolution256×256256\\times 256constructed from two Reynolds numbers,Re=1000\\mathrm\{Re\}=1000andRe=10000\\mathrm\{Re\}=10000\(500 trajectories each\), and tested on the unseen intermediate settingRe=5000\\mathrm\{Re\}=5000with 100 independently generated trajectories\. The checkpoints correspond to0\.2​T0\.2T,0\.4​T0\.4T,0\.6​T0\.6T,0\.8​T0\.8T, and1\.0​T1\.0T\(rollout steps44,77,1010,1313, and1616\), which measure how errors accumulate as the autoregressive rollout proceeds under Reynolds\-number interpolation\. As shown in Table[6](https://arxiv.org/html/2608.11237#A4.T6), GeoIncNO consistently achieves the lowest errors across all rollout checkpoints\. At the final checkpoint1\.0​T1\.0T, GeoIncNO reduces MSE from UNO’s5\.14​E−055\.14\\mathrm\{E\}\{\-05\}to1\.56​E−051\.56\\mathrm\{E\}\{\-05\}, Rel\-L2L\_\{2\}from1\.01​E−011\.01\\mathrm\{E\}\{\-01\}to5\.54​E−025\.54\\mathrm\{E\}\{\-02\}, and Rel\-H1H^\{1\}from4\.67​E−014\.67\\mathrm\{E\}\{\-01\}to2\.57​E−012\.57\\mathrm\{E\}\{\-01\}\. The improvement becomes more pronounced as the rollout horizon increases, indicating that GeoIncNO better suppresses error accumulation and preserves both field\-level and derivative\-level consistency under unseen Reynolds\-number dynamics\.

Table 6\.Long\-horizon rollout errors at different prediction checkpoints under Reynolds\-number interpolation \(trained onRe=1000/10000\\mathrm\{Re\}=1000/10000, tested on the unseenRe=5000\\mathrm\{Re\}=5000\)\.Table 7\.Cross\-resolution generalization results on the 2D NS benchmark\. Models are trained at32×3232\\times 32and tested at the training resolution and finer resolutions\. Lower values indicate better performance\.
### D\.5\.Cross\-Resolution Generalization

We further evaluate cross\-resolution generalization on the 2D NS benchmark\. All models are trained at a low spatial resolution of32×3232\\times 32and tested at the training resolution32×3232\\times 32as well as finer resolutions64×6464\\times 64,128×128128\\times 128, and256×256256\\times 256\. This setting examines whether a model trained on coarse\-grid data can transfer to finer spatial discretizations while maintaining stable autoregressive rollout\. We evaluate both one\-step prediction and long\-horizon rollout performance\. For one\-step prediction, we report MSE@1, Rel\-L2L\_\{2\}@1, and Rel\-H1H^\{1\}@1, which measure the immediate prediction error after one rollout step from absolute, relative, and gradient\-aware perspectives\. For long\-horizon prediction, we report MSE@16, nRMSE@16, Rel\-L2L\_\{2\}@16, and Rel\-H1H^\{1\}@16, which measure the prediction quality after 16 autoregressive rollout steps\. Thus, the @1 metrics reflect short\-term resolution transfer, while the @16 metrics reflect long\-horizon stability under cross\-resolution evaluation\. As shown in Table[7](https://arxiv.org/html/2608.11237#A4.T7), GeoIncNO achieves the best performance at the training resolution32×3232\\times 32across both one\-step and 16\-step rollout metrics, with MSE@1 of6\.03​E−086\.03\\mathrm\{E\}\{\-08\}, Rel\-L2L\_\{2\}@1 of2\.67​E−032\.67\\mathrm\{E\}\{\-03\}, MSE@16 of3\.15​E−063\.15\\mathrm\{E\}\{\-06\}, and Rel\-L2L\_\{2\}@16 of2\.20​E−022\.20\\mathrm\{E\}\{\-02\}\. When transferring to finer resolutions, GeoIncNO continues to attain the lowest error on every reported metric, covering both the one\-step \(@1\) and 16\-step \(@16\) settings, indicating that its advantage is preserved under cross\-resolution transfer rather than being limited to the training resolution\. At the highest tested resolution256×256256\\times 256, GeoIncNO obtains MSE@16 of4\.53​E−064\.53\\mathrm\{E\}\{\-06\}, nRMSE@16 of5\.76​E−035\.76\\mathrm\{E\}\{\-03\}, and Rel\-L2L\_\{2\}@16 of2\.62​E−022\.62\\mathrm\{E\}\{\-02\}, and it also achieves the lowest Rel\-H1H^\{1\}@16 of1\.67​E−011\.67\\mathrm\{E\}\{\-01\}, maintaining both field\-level and derivative\-level accuracy after repeated rollout\. These results suggest that the proposed increment\-centered modeling is beneficial for suppressing long\-horizon error accumulation under cross\-resolution transfer\.

Table 8\.Performance comparison on the FSI task from RealPDEBench under different training paradigms\. Results of baseline models are taken from the original RealPDEBench paper\(Hu et al\.,[2026](https://arxiv.org/html/2608.11237#bib.bib14)\)\.Table 9\.Autoregressive evaluation on the FSI task from RealPDEBench\. Results of baseline models are taken from the original RealPDEBench paper\(Hu et al\.,[2026](https://arxiv.org/html/2608.11237#bib.bib14)\)\.
### D\.6\.RealPDEBench FSI Evaluation

We conduct experiments on the fluid\-structure interaction \(FSI\) task from RealPDEBench\(Hu et al\.,[2026](https://arxiv.org/html/2608.11237#bib.bib14)\), which involves coupled fluid and structural dynamics\. Following the RealPDEBench protocol, we compare models under three training paradigms: simulated training, real\-world training, and real\-world finetuning\. We report RMSE, MAE, Rel\-L2L\_\{2\}, R2, and F\-RMSE, where lower values indicate better performance for error\-based metrics and higher values indicate better performance for R2\. The detailed definitions of all evaluation metrics are provided in Appendix[D\.1\.2](https://arxiv.org/html/2608.11237#A4.SS1.SSS2)\. DMD is included as the original non\-training baseline and is reported only under the corresponding inference setting\.

As shown in Table[8](https://arxiv.org/html/2608.11237#A4.T8), GeoIncNO achieves strong performance with only 9\.8M parameters\. Under simulated training, GeoIncNO obtains the lowest RMSE and Rel\-L2L\_\{2\}\. Under real\-world training, GeoIncNO achieves the lowest Rel\-L2L\_\{2\}and F\-RMSE while matching the best RMSE\. Under real\-world finetuning, GeoIncNO obtains the lowest RMSE and Rel\-L2L\_\{2\}and matches the best F\-RMSE\. These results show that GeoIncNO provides consistent gains across different RealPDEBench training settings\.

Table[9](https://arxiv.org/html/2608.11237#A4.T9)further reports autoregressive evaluation results\. GeoIncNO achieves the best RMSE, Rel\-L2L\_\{2\}, R2, and F\-RMSE among all compared models\. Although U\-Net obtains a slightly lower MAE and Transolver achieves the lowest FE, GeoIncNO provides the strongest overall rollout performance across the main prediction\-accuracy metrics\. Compared with DPOT\-L\-FT, GeoIncNO reduces RMSE from1\.79​E−021\.79\\mathrm\{E\}\{\-02\}to1\.64​E−021\.64\\mathrm\{E\}\{\-02\}and Rel\-L2L\_\{2\}from1\.20​E−011\.20\\mathrm\{E\}\{\-01\}to1\.10​E−011\.10\\mathrm\{E\}\{\-01\}while using far fewer parameters, 9\.8M versus 673\.5M\. These results indicate that GeoIncNO maintains stable multi\-step prediction behavior on the real\-world FSI task\.

### D\.7\.Qualitative Visualization

We further provide qualitative visualizations on competitive benchmarks, including 1D Burgers, 2D NS, 2D SW, and 3D Maxwell\. These examples cover different types of PDE dynamics, such as nonlinear transport, vortical flow, free\-surface wave propagation, and electromagnetic wave propagation\. The visual comparisons are used to examine structural preservation, phase consistency, and error accumulation during prediction\.

##### 1D Burgers

Figure[8](https://arxiv.org/html/2608.11237#A4.F8)compares the predicted solution fields and corresponding error heatmaps on the 1D Burgers benchmark\. GeoIncNO closely follows the smooth spatiotemporal structure of the ground\-truth solution, while several baselines exhibit visible distortions along the temporal direction\. The error maps further show that GeoIncNO produces weaker and more uniformly distributed errors, whereas the baselines contain larger localized error regions\. These results indicate that GeoIncNO better captures nonlinear transport dynamics in the Burgers equation\.

##### 2D NS and SW

Figure[9](https://arxiv.org/html/2608.11237#A4.F9)presents qualitative rollout results on the 2D NS and SW benchmarks att=0,5,15,20t=0,5,15,20\. For NS, GeoIncNO better preserves vortical structures and avoids the structural distortion or misplaced high\-vorticity regions observed in several baselines\. For SW, GeoIncNO maintains more accurate fluid\-height evolution and reduces the phase mismatch, amplitude distortion, and enlarged spatial errors produced by other methods\. These results suggest that GeoIncNO provides more stable rollout predictions for two\-dimensional dynamical systems\.

##### 3D Maxwell

Figure[10](https://arxiv.org/html/2608.11237#A4.F10)visualizes the rollout predictions of the magnetic\-field componentBxB\_\{x\}and the electric\-field componentExE\_\{x\}on the 3D Maxwell benchmark\. For each component, XY, XZ, and YZ slices are shown att=0,5,10,15,20t=0,5,10,15,20, together with the ground truth, prediction, and absolute error maps\. GeoIncNO preserves the main wave structures across different slice planes and maintains small errors throughout the rollout\. By contrast, several baselines exhibit over\-smoothed fields, distorted wave patterns, or enlarged errors at later time steps\. These results show that GeoIncNO better preserves 3D electromagnetic dynamics and phase consistency under long\-horizon prediction\.

![Refer to caption](https://arxiv.org/html/2608.11237v1/Figure/Figure8.png)Figure 8\.Qualitative comparison on the 1D Burgers benchmark, including the ground\-truth solution, predicted solution fields, and absolute error maps\.![Refer to caption](https://arxiv.org/html/2608.11237v1/Figure/Figurens2D.png)\(a\)2D NS
![Refer to caption](https://arxiv.org/html/2608.11237v1/Figure/FigurenSW.png)\(b\)2D SW

Figure 9\.Qualitative rollout visualization on 2D dynamical benchmarks att=0,5,15,20t=0,5,15,20\. The top panel shows vorticity\-field predictions on the 2D NS benchmark, and the bottom panel shows fluid\-height predictions on the 2D SW benchmark\.![Refer to caption](https://arxiv.org/html/2608.11237v1/Figure/Figure19.png)\(a\)BxB\_\{x\}component
![Refer to caption](https://arxiv.org/html/2608.11237v1/Figure/Figure20.png)\(b\)ExE\_\{x\}component

Figure 10\.Qualitative comparison on the 3D Maxwell benchmark\. The left panel shows slice\-wise predictions and absolute error maps for the magnetic field componentBxB\_\{x\}, and the right panel shows those for the electric field componentExE\_\{x\}\.