Merging Memory and Space: A State Space Neural Operator
Authors: Nodens Koren
Samuel Lanthaler
Venue: Accepted by TMLR ()

Abstract: We propose the *State Space Neural Operator* (SS-NO), a compact architecture for learning solution operators of time-dependent partial differential equations (PDEs). Our formulation extends structured state space models (SSMs) to joint spatiotemporal modeling, introducing two key mechanisms: *adaptive damping*, which stabilizes learning by localizing receptive fields, and *learnable frequency modulation*, which enables data-driven spectral selection. These components provide a unified framework for capturing long-range dependencies with parameter efficiency. Theoretically, we establish connections between SSMs and neural operators, proving a universality theorem for convolutional architectures with a full field of view. Empirically, SS-NO achieves strong performance across diverse PDE benchmarks—including 1D Burgers' and Kuramoto–Sivashinsky equations, and 2D Navier–Stokes and compressible Euler flows—while using significantly fewer parameters than competing approaches. Our results demonstrate that state space modeling provides an effective foundation for efficient and accurate neural operator learning.

Published in Transactions on Machine Learning Research (03/2026)
Merging Memory and Space: A State Space Neural Operator
Nodens F. Koren
nodersskorenn@nn.cliia..
Deportment of Computer Science ETH Zürich
Samuel Lanthaler
samuel wwnivie. ac.at
Faculty Mathematics
University of Vienna
Reviewed on OpenReview: http:////eeeerrieae///////o//oor
Abstract
We propose the State Space Neural Operator (SS-NO), compact architecture for learning solution operators of time-dependent partial differential equations (PDEs). Our formulation extends structured state space models (SSMs) to joint spatiotemporal modeling, introducing two key mechanisms: adaplive damping, which stabilizes learning by localizing receptive fields, and learnable frequency modulation, which enables data-driven spectral selection. These components provide a unified framework for capturing long-range dependencies with parameter efficiency. Theoretically, we establish connections between SSMs and neural oper- ators, proving universality thcorem for convolutional architectures with a full field of view. Empirically, SS-NO achieves strong performance across diverse PDE benchmarks including Sivashinsky equations, and 2D Navier-Stoke and compressible Euler while using significantly fewer parameters than competing approaches. Our results demonstrate that state space modeling provides an effective foundation for efficient and accurate neural operalor learning.
1D Burgers' and Kuramoto-Sivashinsk
1 Introduction
Many problems in scientific computing require approximating nonlinear operators that map input functions to output functions, often governed by partial differential equations (PDEss). Neural operators (NOs) (Kovachki ct al. 2023 provide a mosh-indopendeen framework for this task, operating directly on functions and generalizing across discretiiatiions.
The Fourier Neural Operator (FNO) (L et al. 2021 is widely used NO, its global convolution kernels effectively capture long-range spatial correlations. However, this fully global design comes with substantial computational and memory costs, scaling poorly in higher dimensions. Factorized variants (Tran et al. 2023
Lehmann et al. 2023 2024 address this by decomposing operator learning into lower-dimensions subspaces, achieving linear scaling in dimension a crucial advantage for higl-dimennionall spatial applications. Yet, these approaches offer limited flexibility in adjusting receptive fields or frequency modes to different spatial regions or temporal dynamics.
In complementary line of work, state space models (SSMs) Gu al. 2022b), such as the diagonalized S4D architecture (Gu et al. 2022a). have achieved remarkable success in modeling long-range temporal dependencics. Like factorized FNOs, SSMs scale linearly in time and memory, but they provide an additional advantage: data-dependent parameterization of convolutional filters that can learn to emphasize different frequency components and temporal scales. Despite these parallels, the connection between FNOs (primarily spatial) and SSMs (primarily temporal) has remained underexplored.
This paper makes that connection explicit. We introduce the State Space Neural Operator (SS-NO), a unified operator learning framework that applies SSMs directly joint spatiotemporal domains. SS-NO

---

Published in Transactions on Machine Learning Research (03/2026)
can be interpreted in two complementary ways: extending SSMs from purely temporal modeling to spatiotemporal operator learning, generalizing FNOs into more expressive architectures with adaptive receptive fields, frequency content, and temporal causality.
Theoretically, we prove that any continuous operator can be approximated arbitrarily well by convolutional NOs with full receptive fields, a condition satisfied by both FNO and SS-NO. Importantly, the spatial-only variant of SS-NO subsumes factorized FNO as special case while adding two key capabilities: learned damping coefficients for stability control and data-driven frequency modulation. Empirically, we validate SS-NO on challenging 1D and 2D benchmarks, including chaotic Kuramoto Sivashinsky dynamics, variants of Kolmogorov flow, Richtmyer-Messhe SS-NO surpasses existing methods using fewer parameters. These results highlight SS-NO as a principled, scalable, and practical architecture for data-driven modeling in engineering and the physical sciences.
instability, and gravitational Rayleigh-Taylor turbulence. Across tasks,
2
Related Work
2.1
Neural Operators
Data-driven neural operators have emerged a powerful framework for learning PDE solution operators
(Chen & Chen
1995 Bhattacharya et al., 2021 Lu et al. 2021 Kovachki et al. 2023 Theoretical work
established that several neural operators can universally approximate nonlinear operators (Chen & Chen
1995
2024: Lanthaler et al. 2024), including DeepONet. (Lu et al. 2021 Lanthaler et al. 2022) and the Fourier Neural Operator (FNO), which implements global convolution via the Fourier
Kovachki et al.
2021
transform (Li et al.) 2021). Extensions such as the Factorized Fourier Neural Operator (FFNO) (Tran et al.
2023) enhance efficiency and expressivity in practice, but their theoretical expressivity remains unknown. Alternative inductive biascs have also been explored, including U-Net-based convolutional operators (Raonic
et al. 2023 Gupta & Brandstetter] 2023 Rahman et al. 2023 Takamoto el al. 2022) and aotention-basey
operators like the Galerkin Transformer (GKT) [Cao
2021 and the FactFormer (Li et al. 2023a
While some studies suggest that purcly Markovian models can be performant for learning timo-cvolution PDEs
(Tra et al. 2023 Lippe et al. 2023), recent work has shown that incorporating past states improves accuracy, particularly under low-resolution or noisy conditions. This has led the integration of memory mechanisms into neural operators. Buitrago al. (2025) introduce the Memory Neural Operator (MemNO) framework, which embeds temporal S4-based recurrence within general neural operator architectures, motivated by the Mori-Zwanziz formalism (Mori 1965 Zwanzig 1961). Our work extends this direction by examining scenarios with missing contextual information and introducing a spatiotemporal factorization that distributes memory across both time and space.
2.2 Structured State Space Models (SSMs)
Structured state space models (SSMs) have recently achieved strong results on long-range sequence tasks in natural language processing and vision. The Structured State Space sequence model (S4) (Gu et al.) 2022b () uses a continuous-time linear SSM layer to capture very long-range dependencies efficiently, and the Mamba model (Gu Dao 2024) extends this idca by making the SSM parameters input-dependent, further improving performance on large-scale language modeling tasks. These SSM-based architectures have been shown to match exceed Transformer performance on variety of benchmark tasks.
SSMs have also begun to be used in operator learning. For example, Hu et al. (202)) incorporate Mamba-style SSM layers lo predict dynamical systems efficiently. In the context of neural PDE solvers, recent works like Cheng et al. (2024) and Zheng et al. (2024) embed Mamba SSM modules to capture spatial correlations in solution fields. These methods primarily focus on spatial modeling, applying SSM layers across spatial dimensions or treating time steps sequentially. By contrast, our work focuses on applying SSMs jointly in both space and time, employing unified spatiotemporal SSM architecture for neural operator learning.
2

---

Published in Transactions on Machine Learning Research (03/2026)
3 Methodology
3.1 Problem Formulation
Let Ω C Rd be a bounded spatial domain, and consider the solution u E C( T|: (2); of a time- dependent partial differential equation (PDE), where C' denotes the space of continuous functions and L denotes the space of square-integrable functions. We assume access to a dataset of N solution trajectories (1, generated from varying initial or parametric conditions.
In practice, we discretize the spatial domain Ω using an equispaced grid S of resolution f. and the temporal domain [0, using an equispaced grid T with N₂ + points. Let Tin < T denote the fixed input horizon.
Given the trajectory segment u(4, (((,,))<e(oTTm,) on S, our goal is to predict the future evolution u(t, (),,))eee (eET(T..TTT on the same grid. Rather than forecasting single step, we aim to learn the full system dynamics conditioned the initial segment. Formally, the NO model defines parametric map Vo : C(Tiu (9Q; R)) C(Tou; (22 )),, where Ti = [0, mm) and Tout [m T], mapping the input segment to the predicted solution.
3.2 Structured State Space Models (S4 and S4D)
We begin with the contimuous-time linear state space model (SSM) defined as:
d
dt vt) Av(t) + Bu(t),
(1)
y(t) Cv(t) Du(t).
where v(t) E R is the hidden state, H is the number of states, e(t)) E R is the input, and y(t) E R is the output. The matrices RR B E Hx1 Ce RRIH and
define the system dynamics.
The Structured State Space (S4) model (Gu et al. 2022b introduces parameterization where A is structured to allow efficient computation of the convolution kernel k(t) corresponding to the system's impulse response. Specifically, S4 utilizes diagonal plus low-rank (DPLR) structure, A A PQ where A E " is diagonal, and P, Q E #xr with r < H. This structure enables computation of k(t) via the Cauchy kernel and allows for efficient implementation with the Fast Fourier Transform (FFT).
To further simplify the model, S4D (Gu et al. 2022a) restricts A to be diagonal, i.e., A = diag(A,, AN), and initializes the cigenvalues {A.} approximate the behavior of the original $4 model. This diagonalization reduces the compulational complexity and memory footprint while relaining the ability to model long- long-range dependencies.
In both S4 and S4D, after computing the convolution with the kernel s<(), a pointwise nonlincarity a(-) (c.g..
GELU) is applied, followed by residual connection.
3.3 Markovian Spatial State Space Model
To model spatial dependencies in PDEs, we extend the S4D framework to spatial dimensions. A key requirement. for universality of factorized convolutional neural operators is a full field of view, as formalized in Theorem 4.1 A unidirectional SSM does not satisfy this condition because each spatial location only aggregates information from one side, limiting expressivity. Empirically, we confirm Appendix C.1 that. unidirectional scans underperform, motivating the use of bidirectional processing.
For 1D spatial domain discretized into X points, define two SSM models: Myud and Mbwd, where the fund and aad labels serve as identifiers rather than indicating scan direction. The processing is computed as
Dfod M Mum Youd flip ((pp)(())
with the final output y = ffd + Yowd- For multi-dimensional grids, such as X Y, bidirectional SSMs are applied sequentially along each axis (first x, then y), ensuring every spatial location has full field of view and satisfying Theorem 4.1
3

---

Published in Transactions on Machine Learning Research (03/2026)
Comparison with Existing Methods. Unlike Vision Mamba (Zhu al. 2024 which processes 2D data in vigzag manner using bidirectional SSMs, our approach applies SSMs separately along each spetial dimension. While Factorized Fourier Neural Operators (F-FNOO (Tran et al. 2023) process each spatial dimension in parallel and sum the results, our method applics SSMs sequentially, allowing more expressive modeling of spatial interactions. Due to memory constraints, we do not employ Mamba-based SSMs directly but explore alternative 2D and factorized spatial SSM configurations using S4D Appendix D
Forward SSM
Input
Reshape
Output
Shape = (B*T*Y, (, X)
Backward SSM
Flip
Flip
Figure 1: Detailed illustration the spatial bidirectional SSM module. B: batch size, T: temporal length, X: spatial dimensions, C: imput channels, pointwise nonlinearity, and +: element- wise addition. The input is processed through both forward spatial SSM and flipped backward spatial SSM. Each path includes residual connection and nonlinear activation, and their outputs are aggregated to form the final output.
Spatial SSM
Spatial SSM
Temporal SSM
Spatial SSM
Spatial SSM
Input
Projection
Projection
Output
Figure 2: Architecture combining Markovian ID spatial SSM modules with single temporal SSM following the MemNO framework. The spatial SSMs are applied sequentially across spatial dimensions, while the temporal SSM's position within the stack is tunable hyperparameter.
3.4 Full SS-NO Model with Temporal Memory
The SS-NO architocture combines Markovian spatial SSM layers (Soction 3.3 with a single non-Markovian temporal memory module following the MemNO framework to capture spatiotemporal dependencies in PDEs. Spatial layers process local interactions sequentially across dimensions, while the temporal module aggregates information across previous timesteps. During inference, the temporal module maintains and updates a hidden state at each step, which is combined with spatial processing for next-step predictions.
The number spatial layers and the placement of the temporal layer within the network stack are tunable hyperparameters optimizod for specific PDE characteristics. Figures[anm illustrate SSM modules perform bidirectional processing with residual connections, while the full model integrates these with the temporal module. Complete implementation details provided in Appendix G
iinllusrtt the architecture: spatial
4

---

Published in Transactions on Machine Learning Research (03/2026)
4 Theory
We view the proposed SS-NO architecture as an instance of popular neural operator paradigm Kovachki et al. (2023), combining two types of layers: (a) pointwise composition with nonlinear activations, and (b) application of nonlocal convolution operators, and with prototypical hidden layers of the form,
L: v(z) + (We(x) + JD k(z ()(z=-))(()) v) dy).
(3)
(a) nonlinear
(b) nonlocal
where v: D - TR is a hidden state, a is a standard activation function (e.g., GELU), W E RHH is weight matrix, b R a bias vector and k D + k is a learnable integral-kernel.
To give a unified discussion, we will say that V is a convolutional NO, if it is of the form vu = Q 1 o L1 R(u), with hidden layers La as in equation 3 and a choice of parametrized kernel ke(x)) ke(s; 0). In addition, have lifting layer (x))(x) = R(u(z), and projection layer Q(v(z)) Q(v(z)) by composition with shallow neural networks R, We next discuss three architectures exemplifying this approach, highlighting their commonality and giving a unified description. This results in a sharp criterion for universality for any such convolutional NO architecture.
Fourier neural operator (FNO). The FNO parametriies the convolution kernel by a truncated Fourier series, x(x) ((r)=_np - &ge with cut-off parameter K, and [k max |kyl The Fourier coefficients & E " are complex-valuee matrices. The parametrization of the integral kernel x(x) of FNO in d dimensions entails a considerable memory footprint, requiring O(LI(2 parameters. This can render FNO prohibitive in high-dimensional applications.
Factorized FNO (F-FNO). To lessen the computational demands of FNO, so-called "factorized" architec- tures have been introduced. Here, convolutions are taken along one dimension a time. Mathematically, this corresponds to choosing kernels of the form k(x) = (x,)) Ik 6(xx), where K, (2y) is 1d FNO kernel and 8(xx) the Dirac delta function. Thus, integration is effectively only performed with respect to zz. The resulting F-FNO only requires O(LII² Kd) parameters Tran et al. (2023).
Spatial SSM. We next consider the spatial convolution layers of the proposed SS-NO architecture. In this case, solution of the ODE system equation] [ leads to an explicit formula for the corresponding kernel. For a 1d spatial domain, this results in a kernel k(x)) = K4 (x) + (x),,, where ky (z) correspond the backward and forward scans, respectively, and Ka (I) is of the form
Ky (x) -Ip) (x)
Σ Gellll twne C& rk- Real(As), wh Imag(As), BR, Cx E
(4)
Here 1R2 (x) is the indicator function of R2 {+r 0}. The main difference with FNO lies in the parametrization of the convolutional kernel K, which now has tunable frequency parameters w3 E R, and tunable damping parameters < O. The parameter count of the resulting factorized spatial SSM architocture in d-dimensions requires at most O(L + K)d) parameters.
A comparison of the model parameter count entailed by the above choices is summarized below:
Model
FNO
F-FNO spatial-SSM
# Parameter (scaling)
kd
KH
L(H² + HK)d
4.1 A Sharp Criterion for Universality of Convolutional NOs.
As highlighted above, FNO, factorized FNO and SSMs all share common structure, distinguished by the chosen kernel parametrization. What can be said about the expressicity of such architectures? Our first goal is to derive a sharp, general condition for the universality of such convolutional NO architectures.
A minimal requirement for the universality of neural network architectures that the value of each output pixel must be informed by the values of all input pixels, i.c., the model needs have "full field of view". For convolutional NO architectures, this leads to the following definition (ef. discussion in App. B.2
5

---

Published in Transactions on Machine Learning Research (03/2026)
Definition. A convolutional NO architecture with layer kernels R1....;KL KL has full field of view, if the ilerated kernel R == Kr, NI.-1 K1 is non-wanishing, i.e. R(x - y) + 0 for all I; y E D.
The next result shows that this "minimal" condition is actually already sufficient for the universality of a convolutional architecture, requiring no additional assumptions on the convolution operators:
Theorem 4.1. A (faclorized) concolutional NO archilecture is universal if it has full field of view.
We refer Appendix B for the precise statement and proof. To the best of our knowledge, the criterion identified above is both the simplest and most widely applicable result for convolutional NO architectures; it implics universality of (vanilla) FNO, factorized FNOs and combinations of (F-)FNO and SSMs, SS-NO and even gives a sharp criterion for the universality of localized convolutional NOs (similar to the CNO Raonic et al. (2023) ) for which Kę has localized support. A theoretical basis for both factorized and localized architoctures had been missing from the literature. The above result closes this gap. Appendix contrasts empirical results with a unidirectional SSM that violates the criterion.
4.2 Adaptivity and Enhanced Model Expressivity of SS-NO
Based on the formulation in equation the factorized spatial SSM structurally subsumes the FNO (in ID) and F-FNO architectures. We provide a rigorous proof of this capability in Appendix []] Fundamentally, the SSM parameterizes continuous ID convolutional kernel of the form
K Σ
<<x)
cee c -ke--Pllu p
Pk 202 Wk R.
(5)
This formulation generalizes the standard Fourier basis, allowing the model to learn adaptive frequencies and damping rales beyond the fixed harmonics of the DFT. However, since the damping and frequency are tunable parameters within SS-NO, this adds further adaptivity: (1) choice of pk 0 allows effective kernel localization; the model can optimize the support of its convolutional kernel, interpolating between the global kernels of FNO (px R 0) and very localized CNN-like kernels (pk. >> 1). (2) additional adaptivity comes from adaptive mode-fittering; the model can optimize the frequencies Wi...., WK: most relevant for nonlocal processing. Thus, in theory, this added adaptivity implies enhanced model expressivity of SS-NO over F-FNO.
5 Experiment Setup
5.1 Setup and Dataset Description
We evaluate our models on a suite of one-dimensional (ID) and two-dimensional (2D) partial differential equation (PDF) benchmarks commonly used operator learning. These datasets span a variety of dynamical regimes from chaotic behavior to turbulent flows and are designed to assess the models' ability to capture long-range temporal dependencics. For 1D problems, we use datasets based on the Burgers' equation (nonlinear convection ddifuusioon and the Kuramoto Siiasshisskyy equation (chaotic dynamics), following the same data sources as Buitrago et al. 2025) and generating low-resolution versions through uniform spatial downsampling.
widely used benchmarks, including the TorusLi, TorusVis, and TorusVisForce datascts (Li et al., 2021 Tran et al., 2023), as well as compressible Euler benchmarks such the Richtmeyer Meshkov insta bility (CE-RM) and the Rayleigh Taylor instability with gravitational forcing (GCE-RT) (Herde et al. 2024h.
For 2D problems, we evaluate several
The TorusVis and TorusVisForce datasets incorporate variable, time-dependent forcing under randomly sampled viscosities in the range v E (10-5,10-4). Together, these benchmarks test the models' ability to capture turbulence, long-range interactions, shocks, chaotic dynamics, and interface instabilities.
For evaluation, we compare SS-NO against suite of popular bascline models including U-Net (Gupta Brandstetter 2023), GKT (Cao 2021), Factformer (Li et al.) (2023a), and FFNO (Tran et al. 2023 Buitrago
et al. 2025 for 1D benchmarks, and additionally FNO with a two-dimensional kernel (2D FNO) [Li et. al.
6

---

Published in Transactions on Machine Learning Research (03/2026)
2021 for 2D problems. Following the MemNO framework (Buitrago et al. 2025 , augment all baselines and SS-NO with temporal S4 module using memory window of K 4 for fair comparison. The only exception is GKT, for which we use a multi-input variant with K = 4 where the temporal dimension is mixed with features, as we found the model struggled to converge with an additional temporal S4 module.
Data Preprocessing and Training Sctup. We follow standard practices from prior work (Tran et al., 2023), normalizing all data to the range [0, 1] and adding Gaussian noise with variance 10-3 during training for stability. Models are trained to minimize the step-wise normalized relative &g loss.
Unless otherwise specified, all models use four blocks with consistent hidden dimensions. Full details on baseline models, data generation, preprocessing, and training are provided in Appendices and
5.2 Training Objective.
We train the model by minimizing the empirical risk over the dataset of trajectories. Given a loss function f: L (2) x (2; + R, we solve:
Σ (t,) Vos
0- arg min
(6)
Nout
(ETous
where (-) denotes the prediction at time & E Tout As our loss L, we employ the relative L2-error, discussed below. The formulation in equation enables the model to learn the entire future evolution of the system dynamics from finite observed window, rather than stepwise or autoregressive forecasting.
5.3 Evaluation Metric.
Performance reported using the relative lg error:
Relative Ez(u(t, 2), i(t, [|u(t, a(t,,))l2 ||u(t,
(7)
where ||- ll2 denotes the squared norm over all spatial locations and time steps in the output horizon.
All results reported in the main text represent the mean performance, measured by the best validation loss, across five independent runs with different random seeds. To contextualize performance relative to model complexity, we also report the number of paramcters for each model, defined as the total count of learnable parameters. We note that for all 1D experiments, the standard deviation is less than +5% of the reported mean at resolutions 128 and 64, and less than +7% at resolution 32, while for all 2D experiments, the standard deviation below +0.3% in relative & error. The exceptions to this are GKT and Factformer (2D), which exhibit slightly higher variability.
6 Results
6.1
1D Burgers' Equation
We evaluate model performance on the canonical 1D Burgers' equation with v = 0.001 across temporal resolutions of 128, 64, and 32. As shown Figure 3 SS-NO achieves competitive accuracy at all resolutions while demonstrating exceptional parameter efficiency.
Notably, SS-NO delivers superior performance across all tested regimes. Unlike Fourier-based models whose complexity often scales with input resolution and the number of spectral modes, SS-NO maintains consistent architectural footprint independent of the input grid size. This structural advantage is particularly valuable for applications requiring seamless meptoyment across multiple spatial or temporal scales.
The Burgers' equation exhibits smooth dynamics dominated by diffusion and mild nonlinearity, and SS-NO's consistent performance across resolutions highlights its robustness in handling such well-behaved systems. The model's ability to maintain accuracy while drastically reducing parameter count represents a significant advancement. in efficient operator learning.

---

Published in Transactions on Machine Learning Research (03/2026)
Ks (v=0.075)
KS (v=0.1)
10
10-1
Relative la Error
10"
32
128
256
512
32
64
128
256
512
KS (v=0.125)
Burgers (v=0.001)
Relative Ea Error
10-2
10-²-
10a²
256
Resolution (N)
Resolution (M)
U-Net
GKT
FactFormer
FFNO
SS-NO (Ours)
Figure 3: Resolution Dependence Analysis. Relative & error of SS-NO and baseline models on ID benchmarks across varying spatial resolutions (N , 512)). Left three panels: Kuramoto- Sivashinsky (KS) equation with increasing viscosity coefficients {0.075, 0.1, 0.125). Right pancl: Burgers' equation (v = 0.001), limited to N = 128 due to dataset constraints. While U-Net performance stagnates at higher resolutions (N 2 128), operator learning methods generally improve. Notably, SS-NO (Ours) achieves superior accuracy early at 128 and maintains the lowest error floor at 512, demonstrating robust resolution efficiency.
6.2 1D Muramoto-Sivasiinsyy Equation
We next evaluate models on the 1D Kuramoto Sivashinsky (KS) equation across three viscosities (v 0.075, 0.1, 0.125). The KS equation presents challenging benchmark characterized by increasingly chaotic dynamics and multiscale spatiotemporal features as viscosity decrcases.
Performance Analysis. As shown in Figure 3 SS-NO demonstrates consistent. performance improvements, achieving the best results among evaluated methods across all viscosity regimes. the standard resolution of N 128, SS-NO outperforms the second-best model (FFNO) by 22% v 0.075, 42% at 0.1, and 58% 0.125. These improvements are particularly pronounced the low-viscosity regime (v 0.075), where chaotic behavior dominates and traditional methods often struggle to track long-term dependencies.
8

---

Published in Transactions on Machine Learning Research (03/2026)
Resolution Dependence and Spectral Efficiency. To analyze the asymptotic behavior the architectures,
extended our evaluation to higher spatial resolutions of N = 256 and N = 512. In general, we observe diminishing returns in performance improvement for most architectures beyond N = 128, suggesting that the macroscopic dynamics of the system are largely resolved at this scalc. However, a distinct trend emerges when comparing spectral and state-space approaches:
1. Explicit Spectral Scaling (FFNO): The FFNO demonstrates a slight but continuous reduction in error at higher resolutions. We attribute this to the explicit nature the Fast Fourier Transform; as the grid resolution increases, the Nyquist frequency rises, allowing the FFNO capture high-frequency details that were previously truncated.
Implicit Adaptive Learning (SS-NO): In contrast, SS-NO achieves its optimal performance floor earlier, effectively converging at N = 128, and maintains this superior accuracy at N = 512 even with fixed state dimension. This stability supports our hypothesis regarding the adaptive frequency learning capabilities of State-Space Models. Unlike FFNO, which requires finer spatial discretization to resolve high-frequeet features via FFT, SS-NO's recurrent mechanism appears capable of extracting and encoding these critical dynamics efficiently within its latent state space.
This efficiency is further evidenced at coarse resolutions. At N 32, SS-NO maintains 21% improvement over FFNO at 0.075 and 30% at 0.1. These results collectively demonstrate that SS-NO achieves optimal balance of accuracy and conustess--ccoolling in chaotic regimes while maintaining parameter efficiency across varying resolutions.
6.3 Ablation Study
We conduct a comprehensive ablation study to understand the individual contributions of model capacity, learnable frequencies, and explicit damping mechanisms. Four distinct configurations are evaluated: All, where both frequencies and damping are learnable; Damping only, where frequencies remain fixed while damping is learnable; Frequency only, where frequencies are learnable but damping is explicitly to zero (absent); and Fixed, where neither frequencics nor damping are learnable.
Experiments are conducted on three variants of the KS equation with viscosities E {0.075, 0.1, 0.125), all at. resolution 128. Performance is quantified using relative l2 error. Results are presented in Table [ to isolate the effects architectural components across different capacity regimes.
Table 1: Relative l error on KS benchmarks varying state sizes and training settings.
State Size
Setting
Parameters
Relative la Error
KS
V 0.075
V 0.1
D 0.125
All
203,713
0.0086
0.0026
0.0013
Damping only
187,329
0.0086
(-0.0000)
0.0027
(+0.0001)
0.0013
(+0.0000)
64
Frequency only
187,329
0.0110
(-0.0024)
0.0030
(+0.0003)
0.0014
(+0.0001)
Fixed
170.945
0.0116
(+0.0030)
0.0033
(+0.0006)
0.0016
(+0.0003)
All
154,561
0.0090
0.0027
0.0012
Damping only Frequency only
146,369
0.0094
(+0.0004)
0.0028
(+0.0001)
0.0014
(+0.0002)
32
146,369
0.0137
(-0.0047)
0.0038
(+0.0011)
0.0018
(+0.0006)
Fixed
138, 177
0.0158
(-0.0068)
0.0042
(+0.0015)
0.0020
(+0.0008)
All
129,985
0.0115
0.0038
0.0017
Damping only Frequency only
125,889
0.0131
(-0.0016)
0.0044
(+0.0006)
0.0019
(+0.0002)
16
125,889
0.0198
-0.0083)
0.0069
(+0.0033)
0.0031
(+0.0014)
Fixed
121,793
0.0227
(-0.0112)
0.0068
(+0.0030)
0.0031
(40.0014)

---

Published in Transactions on Machine Learning Research (03/2026)
6.3.1 Observations
Explicit damping enables remarkable parameter efficiency. The most striking finding is that proper damping mechanisms allow models to maintain strong performance even under extreme capacity constraints. The 16-state damping-only model achieves performance competitive with much larger models, demonstrating that damping serves as powerful regularization and stabilization mechanism that dramatically improves parameter efficiency. This efficiency is particularly evident in challenging regimes (u = 0.075), where the 16-state damping-only model approaches the performance of fixed-configuration models with 2-4x more states, despite its siguificantly reduced capacity.
Model capacity dramatically affects damping requirements. Higher-capacies models can implicitly compensate for architectural limitations through additional spectral modes, while lower-capacity models rely critically explicit damping for stability. The performance degradation in 16-state models without damping versus 64-state models confirms that damping becomes increasingly essential as model capacity decreases. This progressive dependency relationship underscores damping's 's role as a crucial stabilization mechanism in capacity-connttrained settings.
Damping significance escalates with problem difficulty. The performance gap between models with and without damping grows substantially with increasing chaos (decreasing viscosity). For v 0.075, proper damping provides significant absolute improvement, while for == 0.125, the benefit reduces substantially. This differential effect confirms that chaotic dynamics benefit substantially from explicit damping mechanisms, while smoother regimes can tolerate their absence.
Learnable frequencies provide consistent spectral adaptation benefits. Across all configurations, the frequency-only setting consistently outperforms the fixed configuration, demonstrating the value of data-driven spectral basis adaptation. This benefit is particularly pronounced in constrained settings and complex regimes like KS 0.075, where learnable frequencies allow models to dynamically allocate their limited spectral resources to the most relevant temporal modes. Rather than being constrained to fixed Fourier basis, models can adapt their frequency representations to focus on the specific oscillatory patterns most critical for the target dynamics, leading to more efficient temporal representation learning. While these gains are generally more modest than those provided by explicit damping particularly in challenging regimes- confirm that learnable frequencies provide meaningful improvements in spectral efficiency.
Results are not driven by overparameterization. The consistent performance pattern across capacity levels confirms that our findings are not artifacts of excessive model size. Even the smallest models (16 states) achieve competitive performance when equipped with appropriate inductive biases, demonstrating that architectural design rather than parameter count drives performance improvements.
6.3.2 Damping Analysis: Connecting Ablation Results to Learned Behavior
The ablation results naturally raise the question: how do these performance differences manifest in the actual learned damping behavior? To answer this, we conduct a detailed analysis of the learned damping distributions across different configurations, aggregating results from five independently trained models with different random seeds for each case to ensure statistical robustness.
Architectural constraints drive stronger damping. Quantitative analysis reveals that damping-only models learn 24% stronger mean damping than full models (0.726 vs. 0.586, p= 5.2 x 10-118), representing a 163% larger relative incrcase from initialization (45.3% vs. 17.2% increase from 0.5), as shown in Figure 4 This compensatory strengthening explains the ablation results: when frequency adaptation is disabled, models intensify their damping mechanisms to maintain stability, directly corresponding to the observed patterns in Table 1
Problem difficulty modulates damping strength. Models trained on more chaotic dynamics (u = 0.075)
learn 14% stronger damping than those on smoother regimes (v = 0.125, 0.586 vs. 0.514, p = 5.8 x 10 10-71), with a 537% larger relative increase from initialization (17.2% vs. 2.7% increase from 0.5), as visualized in Figure 5 This dramatic differential learning directly mirrors the ablation findings where damping provides exponentially greater benefits in challenging regimes.
10

---

Published in Transactions on Machine Learning Research (03/2026)
State Size = 64, KS 0.075
Damping Distribution Overview Layer Groups)
Cumulative Distribution Overview Layer Groups)
3.0
0.8
2.5
Damping only
Full: 0.586
Mcan Domping only. 0.726
Initialized: 0.5
1.5
0.4
memory_o
memory_1 1
1.0
memory_2 2
0.2
memory
0.5
0.0
Actua Damping
Actaal Demping
Figure 4: Damping coefficient distributions for full vs. damping-only models 0.075 with 64 states.
State Size = 64 (Fixed)
Damping Distribution Overview Layer Groups)
Cumulative Distribution Overview
Layer Groups)
1.0
3.5
0.8
v-0.12
2.5
v-0.075
2.0
Mean v=0.075: 0.586
Initialized: 0.5
0.4
memory_ 0
memory_1
1.0
momory 2
0.2
memory_3
0.5
0.0
0.00
0.25
0.50
0.75 1.00
1.25
1.50
1.75
2.00
0.00
0,25 0.50
0.75 1.00
1.25 Actual Damping
1.50 1.75 2.00
Actua Damping
Figure 5: Damping cocfficient distributions for 0.125
0.075 models with 64 states.
KS == 0.075
Damping Distribution
Cumulative Distribution Overview (All Layer Groups)
(All Layer Groups)
1.0
3.5
0.8
N-64
Cumulattv Probabiity
2.5
0.6
Mean N=64; 0.586
2.0
Mcan N-16: 0.652
Initialized: 0.5
0.4
memory 0
memory.1
1.0
memory 2
0.2
memory.3
0.5
0.0
0.0
0.00
0.25
0.50
0.75
1.00
1.25
1.50
1.75
2.00
0.00
0.25
0.50
0.75
1.00
1.25
1.50 1.75 2.00
Actual
Figure 6: Damping coefficient distributions at / = 0.075 for different state sizes.
Capacity constraints amplify damping requirements. Lower-capacity models (16 states) learn 11%
stronger damping than higher-capacity counterparts (64 states, 0.652 0.586, p 2.2 x 10-19), with
11

---

Published in Transactions on Machine Learning Research (03/2026)
77% larger relative increase from initialization (30.4% vs. 17.2% increase from 0.5), demonstrated in Figure This intensified damping in capacity-constrained models explains their ability to maintain stability despite reduced parameter counts, directly connecting to the efficiency results in our ablation study.
The consistent pattern across all three analyses supported by overwhelming statistical significance (p < 10-16 in all cases)-ddmonstratess demonstrates that damping serves as a crucial stabilization mechanism that scales intelligently with architectural and problem constraints. Models automatically learn strengthen damping in response to constraints, providing mechanistic explanation for the performance patterns observed in our ablation study. This adaptive damping behavior represents key innovation that enables both parameter efficiency and robustness across diverse dynamical regimes.
Table 2: Relative & error on 2D Navier-Stokes datasets at 64 x 64 resolution.
Architecture
# Parameters
Relative & Error
TorusLi
TorusVis
TorusVisForce
CE-RM
GCE-RT
U-Net
7, 783, 777
0.0611 1
0.0474
0.0558
0.0644
0.0300
Factformer (2D)
1, 006, 433
0.0628
0.0337
0.0395
0.0642
0.0315
GKT
8, 418,, 049
0.0619
0.0404
0.0742
0.0691
0.0172
FNO2D
67, 197, 700
0.0452
0.0466
0.0444
0.0717
0.0155
FFNO
2, 192, 897
0.0409
0.0231
0.0326
0.0688
0.0196
SS-NO (ours)
369, 665
0.0345
0.0218
0.0263
0.0583
0.0138
6.4 2D Navier-Stokes Equations
Navier- Stokes with Fixed Viscosity. We evaluate model on the TorusLi dataset, which contains vorticity fields generated by solving the 2D incompressible Navier Stokes equation on the unit torus with fixed forcing term and low viscosity y = 10-5, This regime is particularly challenging due to reduced diffusion, leading to highly nonlinear and chaotic dynamics.
Table reports the relative & errors various baselines. Among them, FNO2D achieves an error 0.0452. In contrast, SS-NO significantly outperforms this baseline, reducing the error to 0.0345 representing a 23.7%
improvement in predictive accuracy.
This performance gap is notable given the differences in scaling. FNO2D leverages full 2D Fourier transform to model frequency interactions jointly across both spatial dimensions, with cost scaling as O(ND), where D is the number of dimensions and N is the resolution per dimension. In contrast, SS-NO applies separate 1D operators along each axis in a factorized manner, achicving more favorable O(N D) complexity while retaining or improving accuracy. These results demonstrate the advantage of temporal-stale modeling even in time-invariant contexts, highlighting the role of state-space modeling for capturing fluid dynamics low-viscosity regimes. Among several factorized connection variants tested, our configuration delivered the best performance (see Appendix Dj.
2D Navier-Stokes: Problems with Varying Viscosities and Forces. We also evaluate our model on the TorusVis and TorusvisForce datasets from the FFNO benchmark (Tran et al.) 2023), where both viscosity and external forcing vary across trajectories. In TorusVis, viscosity is sampled uniformly between 10-5 and 10-4 4 and the forcing function timo-indopendent. TorusVisForce introduces a time-varying forcing function with phase shift 8 0.2.
Table 2 summarizes the results. SS-NO achieves the lowest errors across both datasets, consistently outper- forming all basclines. This demonstrates that SS-NO effectively captures spatiotemporal dynamics across mixed-viscosily regimes, where larger viscosily values in the data are significantly easier to predict, while still handling the challenging low-viscosity cases.
Compressible Euler Benchmarks (CE-RM and GCE-RT). We also evaluate SS-NO on compressible Fuler benchmarks: the Meettneyyr-MMehhkon instability (CF-RM) and the Rayleigh-Taylee instability with
12

---

Published in Transactions on Machine Learning Research (03/2026)
gravitational forcing (GCE-RT). These problems involve complex interface dynamics and shocks, making them particularly challenging for predictive modeling.
Table 2 shows that CE-RM has relatively high errors across all models, reflecting its intrinsic difficulty. Interestingly, on GCE-RT, GKT and FNO2D perform better compared to other factorized baselines, suggesting that explicitly modeling full 2D spatial interactions can be beneficial for this problem.
Despite being factorized, SS-NO achieves the lowest relative &2 errors on both CE-RM (0.0583) and GCE-RT (0.0138), demonstrating that it can efficiently capture the dominant spatiotemporal dynamics even in problems with shocks and complex interface interactions. These results highlight that SS-NO combines the efficiency of factorized operations with strong modeling capacity, making it competitive with full 2D approaches in capturing challenging fluid behavior.
7 Conclusion and Future Work
We have presented the State Space Neural Operator (SS-NO), a compact neural operator designed to efficiently learn solution operators for time-dependent partial differential equations. Grounded theoretical universality rosult for convolutional operator learning with full field of view, SS-NO generalizes and improves upon prior laclorizations such FNO by integrating adaptive mechanisms including damping for receptive field localization and learnable frequency modulation for data-driven mode selection-wiin maintaining temporal causality. Our spatial-only variant exactly recovers F-FNO, yet SS-NO dramatically reduces parameter count by exploiting a linear O(N. D) scaling compared to the O(ND) complexity of high-dimennsonal convolutions.
Empirically, SS-NO achieves competitive performance across a suite of ID and 2D benchmarks, including Pavier-Stokon and Euler oquations, demonstrating strong accuracy
the Burgers', Sirramoto-SSvashisskee with significantly fewer parameters. In addition, we proposed a dimensionally factorized variant of SS-NO, which scales favorably in higher dimensions and achieves competitive results on challenging 2D fluid dynamics datasets. These findings highlight that state space modeling is not only compatible with operator learning but also offers unique advantages in long-term memory, adaptivity, and efficiency.
While our formulation provides clear benefits for structured spatiotemporal domains, one limitation is the lack of a systematic connection between damping- receptive field localization and irregular geometries. Developing principled extensions of SS-NO to unstruct ured meshes and complex domains remains an important. direction for future work. Beyond this, we aim apply SS-NO to higher-dimensional and multi- -uuti-physics systems, including astrophysical, cosmological, and space plasma dynamics, where efficient and generalizable PDE operators are critical.
References
Kaushik Bhattacharya, Bandad Hosseini, Nikola B Kovachki, and Andrew M Stuart. Model reduction and neural networks for parametric PDEs. The SMAI journal of computational mathematics, 7:121-157, 2021.
Florent Bonnet, Jocelyn Mazari, Paola Cinnella, and Patrick Gallinari. Airfrans: High fidelity computational fluid dynamics dataset for approximating Reynolds-averagel davier-Stoker solutions. Tn Advances in Neural Information Processing Systems (NeurIPS), 2022.
Johannes Brandstetter, Max Welling, and Daniel E Worrall. Lie point symmetry data augmentation for neural PDE solvers. In The International Conference on Machine Learning (ICML), 2022.
Ricardo Buitrago, Tanya Marwah, Albert Gu, and Andrej Risteski. On the benefits of memory for modeling time-dependent PDEs. In The International Conference on Learning Representations (ICLR), 2025.
Shuhao Cao. Choose transformer: Fourier or Galerkin. In Advances in neural information processing
systems (NeurIPS), 2021.
Tianping Chen and Hong Chen. Universal approximation to nonlincar operators by neural networks with arbitrary activation functions and its application to dynamical systems. IEEE transactions on neural networks, 6(4):911-99) 1995.
13

---

Published in Transactions on Machine Learning Research (03/2026)
Chun-Wun Cheng, Jiahao Huang, Yi Zhang, Guang Yang, Carola-Bibiane Schönlieb, and Angelica I Aviles- Rivero. Mamba neural operator: Who wins? Transformers vs. state-space models for PDEs. Xiv preprint 21..0:2410.02113, 2024. URL
http://rrxivv..or/bb///1111 22111/0/111111113133
Albert Gu and Tri Dao. Mamba: Linear-time sequence modeling with selective state spaces, In The Conference on Language Modeling (COLM), 2024.
Albert. Gu; Karan Goel, Ankit Gupla, and Christopher Ré. On the parameterization and initialization of diagonal state space models. In Advances in Neural Information Processing Systems (NeurIPS), 2022a.
Albert Gu, Karan Goel, and Christopher Re. Efficiently modeling long sequences with structured state spaces. In The International Conference on Learning Representations (ICLR), 2022b.
Jayesh K Gupta and Johannes Brandstetter. Towards multi-sptiioeemporalssaaall generalized PDE modeling. Transactions on Machine Leorning Research, 2023.
Maximilian Herde, Bogdan Raonic, Tobias Rohner, Roger Käppeli, Roberto Molinaro, Emmanuel de Bézenac, and Siddhartha Mishra. Poseidon: Efficient foundation models for PDEs. Advances in Newral Information Processing Systems (NeurIPS), 2024.
Zheyuan Hu, Nazanin Ahmadi Daryakenari, Qianli Shen, Kenji Kawaguchi, and George Em Karniadakis. State-space models are accurate and efficient neural operators for dynamical systems. arXiv preprint SF5::2409.0299/ 2024. URL https::/arxxv..rgg/ /bwpgg//bb//29999 2222/23222930111111
Sukjum Hwang, Aakash Lahoti, Ratish Puduppully, Tri Dao, and Albert Gu. Hydra: Ridirectional state space models through generalized matrix mixers. Advances in Neural Information Processing Systems (NeurIPS), 2024.
Nikola Kovachki, Samuel Lanthaler, and Siddhartha Mishra. On universal approximation and error bounds for Fourier Neural Operators. Journal of Machine Learning Research, 22( 76, 2021. URL http: html
Nikola Kovachki, Zongyi Li, Burigede Liu, Kamyar Azizzadenesheli, Kaushik Bhattacharya, Andrew Stuart, and Anima Anandkumar. Neural operator: Learning maps between function spaces with applications to PDEs. Journal of Machine Learning Research, 24(89):1-977 2023.
Nikola B. Kovachki, Samuel Lanthaler, and Andrew M. Stuart. Operator learning: Algorithms and analysis. In Handbook of Numerical Analysis. Elsevier, 2024. doi: 1111///....a..02//.......... 2024..002..0. URL https: / /www sciencedirect. /coo://cienee//rtiicllpp////////
Samuel Lanthaler, Siddhartha Mishra, and George E Karniadakis. Error estimates for DeepONets: A deep
learning framework in infinite dimensions. Transactions of Mathematics and Its Applications, 6(1)ttnac001, 2022.
Samuel Lanthaler, Zongyi Li, and Andrew M. Stuart. Nonlocality and nonlinearity implies universality in operator learning, 2024. URL httpp::://aaxxvoorr/////aa
13221
Fanny Lehmann, Filippo Gatti, Michaël Bertin, and Didier Clouteau. Seismic hazard analysis with a factorized Fourier Neural Operator (F-FNO) surrogate model enhanced by transfer learning. In NeurIPS 2023 Al for Science Workshop, 2023. URL https:// wtwpr/rrrwiieeeeeeoouuuriiieeeee www.neee/formmiixxxxinyrrrBaa
Fanny Lehmann, Filippo Galti, Michaël Bertin, and Didier Clouteau. 3d elastic wave propagation with a Factorized Fourier Neural Operator (F-FNO). Computer Methods in Applied Mechanics and Engineering. 420:116718, 2024. ISSN 0045-7825 doi:
200..mm112...116111. URL hwtp::tts///// //www
com oeg.reeeeiiccccccoooocoeaaattaee
Zijie Li, Dule Shu, and Amir Barati Farimani. Scalable transformer for PDE surrogate modeling. In Advances in Neural Information Processing Systems (NewrIPS), 2023a.
14

---

Published in Transactions on Machine Learning Research (03/2026)
Zongyi Li, Nikola Borislavov Kovachki, Kamyar Azizzadenesheli, Kaushik Bhattacharya, Andrew Stuart, Anima Anandkumar, et al. Fourier Neural Operator for parametric partial differential equations. In The International Conference on Learning Representations (ICLR), 2021.
Zongyi Li, Daniel Zhengyu Huang, Burigede Liu, and Anima Anandkumar. Fourier neural operator with learned deformations for PDFs on general geometries. Journal of Machine Learning Research, 24(388): 1 1-26, 2023b.
Phillip Lippe, Bas Veeling, Paris Perdikaris, Richard Turner, and Johannes Brandstetter. PDE-refiner: Achieving accurate long rollouts with neural PDE solvers. In Advances in Newral Information Processing Systems (NewrIPS), 2023.
Lu Lu, Pengzhan Jin, Guofei Pang, Zhongqiang Zhang, and George Em Karniadakis. Learning nonlinear operators via Deeponet based the universal approximation theorem of operators. Noture machine intelligence, 3(3):218-22) 2021.
Hazime Mori. Transport, collective motion, and Brownian motion. Progress of theoretical physics, 33(3): 1965.
Vinod Nair and Geoffrey Hinton. Rectified linear units improve restricted Boltzmann machines. In The International Conference on Machine Learning (ICML), 2010.
Korbinian Poppel, Richard Freinschlag, Thomas Schmied, Wei Lin, and Sepp Hochreiter. PLSTM: parallelizable Linear Source Transition Mark networks. In ICML 2025 Workshop on Long- Conteet Foundation Models, 2025.
Md Ashiqur Rahman, Zachary Ross, and Kamyar Avizwaeneshelli. U-NO: U-shaped neural operators. Transactions on Machine Learning Research, 2023.
Bogdan Raonic, Roberto Molinaro, Tim De Ryck, Tobias Rohner, Francesca Bartolueci, Rima Alaifari, Siddhartha Mishra, and Emmanuel de Bézenac. Convolutional neural operators for robust and accurate lcarning of PDEs. In Advances in Newral Information Processing Systems (NeurIPS), 2023.
Makolo Takamoto, Timothy Praditia, Raphael Leiteritz, Daniel MacKinlay, Francesco Alesiani, Dirk Pflüger, and Mathias Niepert. PDEBench: An extensive benchmark for scientific machine learning. In Advances én Neural Information Processing Systems (NeurIPS), 2022.
Alasdair Tran, Alexander Mathews, Lexing Xie, and Cheng Soon Ong. Factorized Fourier Neural Operators. In The International Conference on Learning Representations (ICLR), 2023.
Chloe Wang, Oleksii Tsepa, Jun Ma, and Bo Wang. Graph-Mamba: Towards long-range graph sequence modeling with selective state spaces. arXiv preprint arXiw:2/02.00788, 2024.
Haixu Wu, Huakun Luo, Haowen Wang, Jianmin Wang, and Mingsheng Long- Transolver: A fast transformer solver for PDEs on general geometries. In The International Conference on Machine Learning (ICML), 2024.
Jianwei Zheng, Wei Li, Ni Xu, Junwei Zhu, and Xiaoqin Zhang. Alias-free Mamba Neural Operator. In Advances in Newral Information Processing Systems (NeurIPS), 2024.
Lianghui Zhu, Bencheng Liao, Qian Zhang, Xinlong Wang, Wenyu Liu, and Xinggang Wang. Vision Mamba: Efficient visual representation learning with bidirectional state space model. The International Conference on Machine Learning (ICML), 2024.
Robert. Zwanzig. Memory effects in irreversible thermodynamics. Physical Review, 124(4):983, 1961.
15

---

Published in Transactions on Machine Learning Research (03/2026)
A
Appendix
B
Universality of Convolutional NOs
We here provide a rigorous definition of convolutional NOs having "full field of view", and provide the precise statement of our sharp universality condition.
B.1 Full Field of View
We describe the required full field of view property in greater detail, which we informally stated as definition before Theorem 4.1 in the main text.
The following assumptions were implicit in that definition: 'he convolutional NO architecture is of the form
v(u) Q Cr ... o R(u),
where the lifting layer R(u)(x)= == R(((x)).)) x)) and projection layer 2(v(x)) = Q(v(x)) are given by composition with shallow neural networks R, Q of width H (both applied pointwise to the respective input functions). The hidden layers Le are given by
(e(v)(x)= (Weo(x) + be + / ke(z y: (:0))u()) dy)
with channel width H. We will assume that H can be chosen as large as required. For the integral kernel we assume that, for any choice of H and matrix A E H acting on the channel dimension, there exists a setting the parameters 0f, such that ke(x;07) H a scalar multiple of A, ie. the kernel ke(:;0)) is the form
(<( 07) E E HxH ge(x) Ac, ge: D + R some scalar function.
(8)
To our knowledge, all proposed architectures, including FNO, factorized FNO, SSM, as well (UNet-style) convolutional NOs with localized kernels have this property.
The precise definition of the full field of view property is then the following:
Definition (rigorous). An architecture is said to have full field of view, if for any channel width H and A, E HHH with product Ay ... A₁ + O, we can find porameters such that equation holds, and the iterated kernel
R(x) (KL * xx)) = (gL ....... g1)(2)) Ar A
(9)
has full support, in the sense that the suppor of the function (x, y) + k(x - y) is D D. Remark In particular, R satisfies the full field of view property, if is nowhere-vanishing, i.c for any y, D, we have 6(I y) + O.
B.2 Intuitive Explanation of the Full Field View Property
Before stating our universality theorem, based on the "full field of view" property, we discuss its intuitive meaning in this section.
A general (non-linear) operator G u + G(u) takes input function u(x) and maps it to another function v(x)= = (x)) In general, such operators are non-locol, meaning that the value of the output function v(x) at point x depends on the values of the input function u(y), for y E D, across the entire domain D.
A simple example of an operator with non-local dependency is an operator with a general integral kernel:
((w)(x)= == J. k(I,, v)u(y) dy.
We recall that the support general function f(z) is the closure of the set {=1[(=) #0).
16

---

Published in Transactions on Machine Learning Research (03/2026)
To determine the value of the output function at for such G, it is clear that we generally require knowledge of u(y) at all points y E D. While this illustrative example defines linear operator 9, more general forms of non-locality are present in most operators of interest (such as solution operators of non-linear PDEs), which generally are both non-local and non-linear.
Thus, it intuitively clear that a universal operator-learning architecture must provide mechanism that enables information to flow from the input function u(y) at any point y D to the output function v(z) any D.
For the convolutional architectures discussed in this work, the only mechanism to communicate information non-locally during the forward pass is via the convolutional integral kernels v(x) + I Ke(2 - y)v(v) dy- The iterated kernel equation is a mathematical way to capture the flow of information across all layers: Indeed, if x(œ y) # O, this means that (at least some) information about the imput function at y can flow to the output function x E D. However, a priori, is unclear whether this specific mechanism to transfer information non-locally is sufficiently informative. Indeed, one might assume that the x-dependency of the kernels <²(x)= ge Az, i.e. the functions ge(z)) themselves, need to be very specifically tuned to ensure all relevant information on the input side to correctly reach the output side.
The main theoretical insight of the present work is to show that no specific tuning of ge(x) is required for universality. In fact, any functions 91..... 9L will do, as long as they ensure a full field of view and as long as we can tune the matrices A,...., Ar E IR xx (and allow arbitrarily large hidden channel dimension H). Given this discussion, we view the "full field of view" property as a "minimal" property to ensure the needed flow of information through our architecture. In Theorem B.1 below, will rigorously prove that this "minimal" property is alrcady sufficient.
B.3 A Sufficient Condition for Universality
We now state and prove our universality result for general convolutional NOs.
Theorem B.1. Let G C(D; B + IDD: be continuous operator, with either D C Rd compact or D = T the periodic torus. Let K. c C(D;R4) be a compact set of input functions. Let (u)) be a convolutional NO architecture with a full field of view. Then for any > 0, there erists a channel width H, and a setting of the weights of v, such that
sup sup (G(x(()) - (()(x))s Se. =ED
B.3.1 Proof of sufficiency of "full field of view" property (Theorem B.1
Proof. For simplicity, we will assume that d, = do The proof easily extends to the more general case. The space of continuous operators G K C(D), u(x) + (z)(,,, is isometrically isomorphic to the space C(K's x D), by identifying ((u,z) == (z)) By assumption, both K. and D are compact sets.
Given the choice of convolutional NO architecture, we now fix KI,..., such that R is nowhere vanishing. Let us now introduce the set consisting of all operators that can be represented by a choice of channel width Hl and a choice of tunable parameters:
= = C .... C1 o R, v is a convolutional NO with kernels KI...., KL and channel width H
A A=
Our goal is to show that A C(K x D) is dense. We denote by A the closure of A in C(K D) (set of limit points). To prove that A c C(K D) is dense, will use the following result, which follows from the Stone- Weierstrass theorem:
Lemma B.2. subset CCC x D) is dense, if
The constant function 1 A,
separates points: For any (u,, 21), (uz, x2) E D, there exists E such that (u x1) # V X2).
17

---

Published in Transactions on Machine Learning Research (03/2026)
& is an approrimate vector subolgebra:
A is closed under addition and scalor multiplication (i.e. vector subspace),
A approzimately closed under multéplication: For any 11 22 the product V1 12 E h, where
(V1 (2) =
is the pointwise multiplication.
Our goal is to show these properties for A to conclude that A C CCC x D) is dense.
A contains constants. It is very easy to show that the constant function 1 E A, by defining a convolutional NO that disregards the input and has constant output =1
A is an approximate subalgebra. To show the other properties, we first note that is closed under scalar multiplication and under addition, i.e.
II U
+ 22 EAA, E,, 12 E.
+
If H1 and H2 are the channel widths of Vi and 12, respectively, this conclusion follows by simple parallelization of 11 and 22, and employing the last projection Q-layer to sum (and scale) the parallelized results. Thus, is a vector subspace of C(K x D).
show that A is approximately closed under multiplication, iee. 1 2/2 E implies that Vi V22 E A, we arbitrary 11, 2 E A. We will denote by V the NO v for j 1,2, but where the last projection-layer has been removed:
v1 = Q1 11,
1 = Q2 2
Assuming wlog that the channel width Hi -2= H (otherwise, we can pad the channel width by zeros), note that these incomplete NOs v and 62 define continuous mappings K. D - R Since K. x D is compact, it follows that also the images 1 x D) c R are compact (since compact sets mapped to compact sets under continuous maps). In particular, there exists B > 0, such that
@ (u. (2(u, 2) ( B|,, VuEK,
Consider now Q₁ and Q2, ie. the last layers of 1 and V2, respectively. The following product mapping
[|B x -I-B, + R, (21, + (vi) @2(v2),
is continuous. By the universality of shallow neural networks, for any > 0, there thus exists a shallow net Q 2 + R, such that
sup Jes 1
1Q(1 02) - Q1(1r) 22(12) Sc.
Define now new NO as (u, = ( (1, x)(,, i.e. Q applied to the parallelization of i and V2 Then
sup - ( V2 (2(u,x)) K
|Q( (u, x), (2(u, - (1, (2((2(u,=))
sup
KxD
sup
Q(v 22) - (1())) (x(r)
Se.
Since e was arbitrary, this shows that V V2 is a limit point of v e A, thus V, V2 A. We conclude from the above that A approximate vector subalgebra.
18

---

Published in Transactions on Machine Learning Research (03/2026)
A separates points. To conclude our proof, it only remains to show that A separates points. Let (11, + (2) be two distinct elements in K x D. There are two cases: Fither + 22, 22 and # W2.
Case 1: 21 + 22. We want to construct v E A. such that (ui, 21) & V (ug, Ø2). This is easy, since we can always use the lifting layer R to eliminate the dependency of V on the u-variable. Once the weights acting on the u-variable have been set to zero, we then have that ((, = +xxx) is an ordinary multilayer perceptron. Since 21 + 22: (1) (x1) # (x2)) L (u2, x2).
easily choose such that c.g. (b(1)) = 0 and (x22) 1. Thus
Case 2: T1 x2 x- and 21 A U2- This is the more difficult case. In the following argument, we will appeal to the universality of shallow neural net works, and their compositions, multiple times. We will forego the tedious e-Ō estimates, and instead sketch out the general idea; filling in the details is a straightforward exercise that would not provide any additional insight.
The underlying idea is that roughly by linearization of the hidden layers of the underlying non-linear convolutional architecture we can construct v which approximately takes the form,
(())()) = / g(z 11x-vp (( dy
J (gr. -1--1 *............................... 91) - Ri (u(u), v) dy
where R is an arbitrary multilayer perceptron. Given this simplified form, and since g has full support by assumption, turns out that we can always construct suitable R₁ that ensures that the we can separate the input functions un and u2 ultimately ensuring that (u))(()) o V (uz)(x)) We provide the details of the construction below, starting with the "linearization trick".
Linearization trick: We construct l as follows: First, given matrices A, and functions ge as in the (rigorous) definition of "full field of view", we choose the weights of the hidden layers Le, such that
(=())(z)= = -0 (A1 ge( - y)v(g) dy +bb) be).
This is possible by assumption on the convolutional NO architecture, cp. equation g By universality of shallow neural networks, we can choose H sufficiently large, and choose matrices C., A. and biases O.,, S., such that for all relevant input vectors ę = (E1, .)) (effectively one-dimensional), we have
C.a (A.E 3.) + +a. t0. E,
(10)
i.e. the resulting shallow neural network is an approximation of the identity on one-dimensional &, any desired accuracy.
We now momentarily focus on the case L = 2. In this case, we choose A = A., bi = b.. then choose Ag = A.C and bias by a, f - 3) dy + Be, so that
(2)((()((z) o A2 sols dy
where
A2 92(2 - (c)))) Ar / (C + dy b..
We now note that, by construction, we have for any hidden state function of the form v(y) = [ui(y), "..")......."
C.Li( +a, = CCo (A. w)(() + B.) a. R 1 g1 ())=)))
19

---

Published in Transactions on Machine Learning Research (03/2026)
where we have used equation 10 and note that this approximation is to any desired accuracy. It follows that also
A, / g2(2- ---) y) / 91 2)((2) dz dy + b.
C2 (1()(x)) a
= A. / (1)== y)v(y) dy + b.)
to any desired accuracy. This shows that for two hidden layers L1, La there exists a choice of the parameters, such that for hidden state v(x) (2),,,,,.... we obtain arbitrarily good approximation
An /A (92 /// gi)(2 y)u(y) dy + b.).
L2 o (i()()))
Tterating this argument, for general L 2, we start from Lt. Lr.-1 ... Li, and we first choose the weights of LL and CL-1 such that
(alle--)(())) - v)u(y) dy - b.).
An / (9r.
o (
91. 1)(x (9.
This can be done by the argument employed for the case = 2. Next, apply the same argument again to choose the weights of LI-2, such that
L1. Lt.-1 LI. (z)) (Cr, o ot. Lr-
(A. A,
(gL IL-1 * (0-22( g)v(v) dy by
and continue similarly, until
2 A. f g1) -nlo--)v())())) dy b.).
C, o.. 22. (i(v)(x)
any desired accuracy. This argument is based on the assumption that v(x) (x)),,,..... We next add projection layer Q to this composition, of the form Q(v)(x) Q(C.v(()) ta.), where Q is shallow neural network, to obtain
Qo o... o (1())(z) Q (/, g(z v)u(e) dy). g(z) = (gl *.... (x))
Recall that the support of (I, + is D D, which, roughly speaking, might be thought of g(r y) for all x,y E D. The non-vanishing support is by assumption (full field of view).
Pre-composing with lifting layer (x))(z) = R(u(z), which we choose to have only a non-vanishing first component on the oulput-side, i.e. R(u(x),) = (), 0), it follows that we can construct (x)(x) = Q LL ....o Li o R, such that
(w)(r)~ ( (()))))((u())).)
where the approximation error can be made arbitrarily small. Here, we are still free to choose Q and Ri- Our goal choose them such a way that for our given functions in * ug and the given point x² 21 xz E D, we have (x1))())) (wz)(=)) + (uz)(=*) (2z)) We will choose Q as an approximation of the identity on the first component, that
(()((a)) = / g( Ri((())) dy.
20

---

Published in Transactions on Machine Learning Research (03/2026)
Separating (u1,21) from (u2, when U1 F U2, 21 D2 x=: By assumption, g(2" y) +
has full support for y E D. Since Mı * u2, there exists yo E D such that "I(2o) # w2(w)). We may wlog assume that (16) < 71 < < u2(20) for some TI, 72 E R. By contimuity of u1, U2, there exists > 0, such that (y(o) U1 (1) < T1 < 72 < (y(o) u2(Y)). By the universality of the shallow network R1, we can choose R1, such that we approximately have
(1( v) R psly (r
where ps(. yo) is a non-negative function supported inside Bs(yo), and (-) is non-negative function, such that
1,
(n > T2):
h A (n) =
0,
(1) < <))
Since y - 0(x - v) has full support, and since Po(y yo) is an arbitrary function apart from the constraints we imposed above, we can further refine our choice of po: to ensure that
J, g((a - -ypss - vo) dy + O.
It then follows that and
(1(1(y)))) poly y
yo) no (u1(y)) 0,
R1 (2((())))) p5ly yo (u2(v)) ps(y yo)-
In particular, we conclude that to any desired accuracy we can construct v(u), such that
v(un)(=') = f,0(0- ()
and
((aa)(=') /
| gla" RR (uz(v),,))
J g(x y)po yo) dy # 0.
In particular, given the above two approximate identities, upon making the approximation errors sufficiently small, it follows that there exists V E A such that (ur)((z)) = (=))(()) * (u2)((z)) = ((z) This finally shows that A separates points inputs (u,,zi) and (u2,22), under the assumption that #n j u2 and 21 = 22 = x", and concludes our proof for Case 2.
Our proof of the universality of A C C(C x D) now conchudes by application of the Stone- Woine-Weissssaass thcorem (cp. Lemma B.2
C Further Ablation Studies
C.1 Unidirectional Spatial SSM
We conduct an ablation study to evaluate the importance of the bidirectional spatial module in our SS-NO architecture on 1D Burgers' Equation. To ensure a fair comparison that isolates the effect directionality, we construct a unidirectional model with the same parameter count as the full SS-NO. This done by stacking layers of forward-only spatial SSM blocks and inserting a single temporal SSM layer the middle.
As shown in Figure | the unidirectional variant fails to capture the dynamics effectively. Its training curve plateaus early, and the relative & error stagnates around 0.075. In contrast, the full bidirectional model continues to improve and ultimately reaches a much lower error of approximately 0.007-more than ten times better. These results highlight that unidirections ional spatial SSMs are fundamentally limited in their representational capacity for 2D PDEs and act as non universal approximators in this setting - indeed, this unidirectional spatial SSM is lacking a full field of view, violating the criterion for universality in Theorem Bidirectional context is essential to capture the long-range spatial dependencies needed for accurate forecasting.
21

---

Published in Transactions on Machine Learning Research (03/2026)
Comparison of Unidirectional vs. Bidirectional Spatial SSM
0.200
Unidirectional Spatial SSM
Bidirectional Spatial SSM
0.175
0.150
Relative L2 Full Loss
0.125
0.100
Nu Will
0.075
0.050
0.025
0.000
100
200
300
400
500
Epoch
Figure 7: Training curves comparing the unidirectional and bidircctional spatial SSMs on the TorusLi dataset. The unidirectional variant plateaus early, failing to caplure global spatial dependencies.
C.2 Effect of Missing Contextual Information and the Role Temporal Modeling
To better understand the role of temporal modeling in architecture, we evaluate the full SS-NO against ablated variant that removes the temporal state-space module (S4), well as two strong baselines (FFNO and FNO2D). Figure S shows the relative & error across three levels of contextual information: all contert, only forcing, and no context, on both the TorusVis and TorusVisForce datasets.
Three main observations emerge from these results. First, the temporal S4 component is necessary in the most challenging No Contert setting, where both viscosity and forcing inputs are missing. Without temporal modeling, SS-NO suffers from very high errors (e.g., 0.1156 vs. 0.0411 on TorusVis), whereas the full model remains comparatively robust by leveraging temporal dependencies to compensate for missing information. These results exlend the observations made by Buitrago el al. (2025). who showed that temporal memory helps mitigate the effects of resolution loss and noisy data. In our case, we show that memory-based modeling also helps compensate for missing physical parameters highlighting an additional benefit of incorporating temporal memory into spatiotemporal PDE models.
Second, in the Only Force setting, the presence of forcing information appears be sufficient for accurate predictions in our architecture; both the full and ablated SS-NO achieve consistently low errors, while FFNO exhibits slight degradation. This suggests that temporal modeling is not strictly necessary when forcing is available, even the ablated variant performs competitively.
Finally, in the extreme case where all context is missing, FFNO surprisingly outperforms SS-NO (full), despite performing worse in other settings. We hypothesize this could be the result of the difference between how FFNO and FNO/SS-NG make residual connections, since seems like the degree of performance loss is similar between SS-NO (full) and FNO2D, but we do do further investigation as is of the scope of this work.
C.3 Effect of Varying the Memory Window Size
We investigate how the choice of memory window size K affects the performance our SS-NO model on the KS dataset. with 0.075. This hyperparameter controls how many past spatial feature maps accessible
22

---

Published in Transactions on Machine Learning Research (03/2026)
Performance Under Varying Context Availapility
TorusVis
TorusVisForce
0.12
Relative Error
SS-NO (w/o temporal 54)
0.08
S5-NO (full)
FFNO
FNO2D
0.04
0.00
Context
Only
Context
All Context
Context
Figure 8: Relative l error across three levels of contextual information: full context, only forcing, and context. Lower is better.
the model during temporal state updates. All results reported in this section are obtained from models trained without teacher forcing to better reflect deployment conditions.
As shown in Figure increasing K consistently improves performance up around K 8. Beyond this point, the performance gains become increasingly marginal, and the validation loss begins to plateau. Notably, setting K = O, which corresponds to no temporal memory, results in significantly degraded performance. These results highlight the critical importance of temporal memory in modeling complex spatiotemporal dynamics, while also suggesting diminishing returns beyond a moderate window size.
Validation Loss over Epochs for Different Memory Window Sizes
0.200
Memory Window Size
0.175
= 0
K
K
0.150
K = 12
Relative L2 Full Loss
K
0.125
K
0.100
0.075
0.050
0.025
0.000
200
300
Epoch
Figure 9: Validation l2 loss over training epochs for different memory window sizes K on the KS dataset. (y 0.075). Increasing the window improves performance until around K = 8, after which gains diminish.
C.4 Effect of Teacher Forcing on SS-NO Performance
Case Study: 1D KS with 0.075 at Low Resolution. While prior work such as FFNO Tran et al. 2023) suggests that teacher forcing can improve model performance by stabilizing training, we observe a different trend for our SS-NO architecture, particularly in low-resolution regimes. To investigate this,
23

---

Published in Transactions on Machine Learning Research (03/2026)
Validation Loss over Epochs for Different Memory Window Sizes (Teacher Forcing)
0.200
Memory Window Size
0.175
=
K
K 8
0.150
K
Relative Full Loss
K
0.125
K = 20
0.100
0.075
0.050
0.025
0.000
100
200
300
400
500
Epoch
Figure 10: Validation & loss over training epochs the KS dataset (v 0.075) using teacher forcing for various memory window sizes K. Compared to Figure g training is less stable and results in higher loss.
revisit the 1D Kuramoto Sivashinsky (KS) dataset with v 0.075, where reduced spatial resolution poses a significant challenge.
We train SS-NO with teacher forcing across varions memory window sizes (K 0 to 20), and present the results in Figure 10 Compared to the non-teacher-forced variant (Figure g training with teacher forcing leads to noticeably less stable optimization, slower convergence, and consistently higher final relative & errors. Morcover, while performance still improves with longer memory, the gains platcau carlier, and the benefits beyond = 8 diminish more rapidly than in the non-teacher-forceo case.
These findings suggest that in this setting, SS-NO benefits more from fully autoregressive training, likely due to its strong inductive bias for modeling long-term temporal dependencies without relying on ground-truth guidance.
General Trends Across 1D and 2D Benchmarks. To better understand the broader impact of teacher forcing, compare SS-NO models trained with and without teacher forcing across all 1D and 2D benchmarks (Tables and
In ID settings, results are mixed: at higher resolutions (e.g., 128), teacher forcing slightly improves performance in some cases (e.g., KS with / = 0.1 and v = 0.125), but at lower resolutions, non-teacher--orte models generally perform better- especiallyy in the challenging autoregressive training may confer more robustness in low-resolution or harder regimes, where model predictions diverge more easily from ground truth.
0.075 KS setup. These results suggest that
In contrast, across all 2D Navier-Stokes benchmarks, teacher forcing consistently yields better performance. This suggests that for ligher-dimenusoonal systems with complex spatiotemporal interactions, teacher forcing can help stabilize training and guide the model toward more accurate trajectories.
Overall, while the effect of teacher forcing appears to be dataset- and resolution-depndeent, we find that non-teacher-forces training can be highly effective in 1D regimes, particularly under low resolution. In contrast, for complex 2D flows, teacher foreing remains valuable tool for improving predictive accuracy.
24

---

Published in Transactions on Machine Learning Research (03/2026)
Table
Relative & error on ID benchmarks at varying resolutions.
Resolution
Architecture
Relative & Error
Burgers'
0.075
v == 0.1
V == 0.125
= 0.001
SS-NO (Teacher Forcing)
0.0086
0.0027
0.0013
0.0070
128
SS-NO (No Teacher Forcing)
0.0076
0.0034
0.0019
0.0070
SS-NO (Teacher Forcing)
0.0143
0.0048
0.0029
0.0121
64
SS-NO (No Teacher Forcing)
0.0135
0.0042
0.0023
0.0114
SS-NO (Teacher Forcing)
0.0472
0.0162
0.0088
0.0200
32
SS-NO (No Teacher Forcing)
0.0358
0.0135
0.0071
0.0171
Table 4: Relative l error on 2D Navier Navier-Stokes datasets. All models are evaluated at a spatial resolution of
x 64.
Architecture
Relative l Error
TorusLi
TorusVis
TorusvisForce
SS-NO (Teacher Forcing)
0.0345
0.0218 0.0385
0.0263 01.071
SS-NO (No Teacher Forcing)
0.0546
Comparison of Factorized Architectural Variants
In this section, we explore how different architectural factorization strategies affect performance on the TorusLi dataset. Specifically, we examine two popular alternatives to our default SS-NO spatial block.
SS-NO-VM (Vision Mamba style). This variant adapts the Vision Mamba architecture by replacing the Mamba core with our spatial SSM block. The key idea is to preserve Vision Mamba's flattened zigzag scanning structure in 2D, allowing bidirectional processing over a 1D sequence obtained from flattening the spatial grid. After bidirectional processing over the first spatial axis, the result is passed through a mixing linear projection before second pass along the alternate axis, again using bidirectional spatial SSMs and linear mixing. This model maintains Vision Mamba's emphasis alternating directional fusion while leveraging the structure of our SSM.
SS-NO-FF (FFNO style). In this variant, we discard the original spatial block of SS-NO and instead adopt the FFNO-style connection, where spatial context is aggregated through two 1D sweeps: one across the z-axis and one across the y-axis. Each sweep involves forward and backward spatial SSM applied independently along that axis, and the outputs from both axes are summed together. This approach highlights the role of directional decoupling in factorized Fourier models and contrasts with the fused zigzag sweep in SS-NO-VM.
Pseudocode for both SS-NO-FF and SS-NO-VM variants provided in Figure
FNO2D-R (Reduced Fourier Kernel). To verify that the gains of SS-NO are not trivially due to kernel simplification, we construct a reduced FNO2D model, denoted as FNO2D-R. this variant, we remove the output channel mixing in the Fourier kernel by modifying the kernel from:
weifhs nn. ParPaaamtterrsell scale torch. rand( in_channels, out_channels, modes1, modes2, baape=ttrrc.fffllatt..
return torch. einsum("bixy, ioxy->boxy", input, weights)
25

---

Published in Transactions on Machine Learning Research (03/2026)
SS-NO-VM (Vision Mamba Style)
SF-NO-FF (FFNO (FFNO-style)
x_flat flatten_spatial(N
residual
First zigzag sweep
# X-axis sweep
x_ fwd = SSM_flat_ forward (x_flat)
x1 = sSM_k_forward(x1)
x x_bwd = wlwt_backaaarddeecee
x1 x1_fwd = SSM_x_
21 permute_channell (x_ fwd d + x_bwd)
x1_bwd = reverse bew... (x1)))
linearT(z1) x_flat
xl_combined = reshape_backkx1 x1_bwd)
Second zigzag sweep x_fwd2 ssM_flat_forwarde
Y-axis sweep
y1
x_bwd2 fst_flat.. hew._aaaaarrdrreeerrrr
y1_fwd ssM_y_forward
permute_channels (x_ fwd2 x_bwd2)
y1_bwd reverse baw..yy.aacaarrrreeereeee
linear12(z2)
yl_comiined reshape_ y1_bwd)
x_out reshape_spatialee
Combine and project
z x1_combined + y1__cmbines
Listing 1: SS-NO-VM (Vision Mamba style)
z backcast_ff(2)
output 2 residual
Listing 2: SS-NO-FF (FFNO-style connection)
Figure 11: Pseudocode for two spatial architectural variants of SS-NO: Vision Mamba-style (left) and FFNO-style (right).
to:
self.weights nn. aarae * forch.rands in_channels, modes1, modes2, (xype=torch..flott)
return input weights
This isolates the contribution of full channcl-wise Fouricr mixing and allows for a more controlled comparison against SS-NO.
As shown T our original SS-NO achieves the best performance among all tested variants, demonstrating that both our architectural choices and structured state-space modeling contribute meaningfully to the observed improvements.
Table 5: Relative la error on the TorusLi dataset for different factorized architectural variants. All models are evaluated at spatial resolution of 64 x 64.
Architecture # Parameters
Relative l Error
TorusLi
SS-NO-VM
402,945
0.0495
SS-NO-FF
503,297
0.0403
FNO2D-R
1, 1,115, 841
0.0718
SS-NO (ours)
369,665
0.0345
E Data Preprocessing and Training Details
E.1 Data Preprocessing and Training Setup
We normalize all input data to the range [0, 1] and add fixed-variance Gaussian noise during training. Through empirical validation, we found a 0.005 optimal for Burgers' equation stabilize training while preserving signal integrity, while 0.001 works best for all other datasets. All models are trained for 500 epochs using
26

---

Published in Transactions on Machine Learning Research (03/2026)
the AdamW optimizer with an initial learning rate of 10-3, weight decay of 10-*, and cosine annealing learning rate schedule. We employ teacher forcing during training unless otherwise noted.
To simulate low-resolution scenarios, we generate downsampled versions of the original datasets by uniformly subsampling in space. All models are trained to minimize the normalized step-wise relative &g loss, and evaluation reported using the full-sequence relative l error.
E.2 Model Architecture Specifications
All models use four blocks with a consistent hidden dimension of 64 across both 1D and 2D problems. Intermediate and output linear layers maintain a hidden size of 128. With the exception of FFNO and our SS-NO on ID problems, we adopt the hyperparameter settings reported by [uitrago et al. (2025 Unlike their setup which uses hidden dimension 128 for 1D and 64 for 2D problems, we found no statistical difference in performance with a unified hidden dimension of 64, simplifying architecture design while maintaining competitive results.
All models except GKT follow the MemNO framework and incorporate a temporal S4 module with window size K = 4 positioned in the middle of the network stack. For GKT, we use a multi-input variant where the temporal dimension is mixed with features.
We adopt a simple spatial positional encoding scheme shared across all models. In 1D, for a grid with f equispaced points over the interval [0, L), the positional encoding E E Rf is defined as E, = f. for S i << f. In 2D, for a ff x f grid over [0, Lg] x [o, Lyl. the encoding is Eij = ( +.) The input lifting operator Tin is shared linear layer that maps the concatenation of input features and positional encoding to the hidden space, applied pointwise across the spatial grid.
E.3 Baseline Model Architectures
Unless otherwise specified, all modcls use four blocks with channel width of 64. Intermediate and output linear layers have a hidden size of 128 (i.e., twice the channel width). With the exception of FFNO, multi-input FNO, and our SS-NO on ID problems, we adopt the "optimal" hyperparameter settings (i.e., number of layers, hidden dimension, and expansion size of linear layers) reported by Buitrago et al. (2025). The U-Net we use is even larger, containing more parameters than this standard baseline. For memory- memory-augmenteed architectures, the memory module is inserted after the first two blocks and before the last two. For all non-Markovian models, we fix the memory window size to K across experiments; the impact of varying K on SS-NO is explored in Appendix C.3
We adopt simple spatial positional encoding scheme shared across all models. In 1D, for grid with f equispaced points over the interval [0, L), the positional encoding E E R is defined as E, = for 0 S sii < In 2D, for f x f grid over [0, T.]] x [0, Tyy, the encoding is F.jj = ( 2,). The input lifting operator Pin is shared linear layer that maps the concatenation of input features and positional encoding from R2+h, where I is the number of features, to the hidden space R², applied pointwise across the spatial grid. Similarly, a shared decoder Your maps from 1 back to R. For all models containing an FNO or FFNO module, we follow the setup in Buitrago et al. (2025) and retain all available Fourier modes by setting kmax = [6 where f denotes the spatial resolution. A notable exception the TorusLi dataset, where an 8-mode FNO configuration perfoormee higher mode counts (16, 24, and 32). We therefore present the results based on 8 modes for this specific case.
Factorized Fourier Neural Operator (FFNO). The Factorized Fourier Neural Operator (FFNO),
introduced by TTran et al. (2023), extends the original Fourier Neural Operator (FNO) (Li et al. 2021) by modifying how the spectral integral kernel is applied. We use the standard FFNO implementation from Each FFNO layer operates on a spatial grid of size |s| and a hidden dimension h, and is defined as a residual block:
https://githhb. cow/alasdairtran/ wwwraa//foorrrrllooo
((v) + Linear o Linearw-h o K(N),
27

---

Published in Transactions on Machine Learning Research (03/2026)
where denotes the ReLU activation function Nair & Hinton 2010 and h is an intermediate dimensionaliity used within the nonlinear mapping. The integral operator K. transforms v in the Fourier domain and is computed by aggregating across spatial dimensions:
K(w) =
Σ IFFT RRa FFTa(v)),
where FFT, and IFFT,, are the (inverse) discrete Fourier transforms applied along the a-th spatial axis, and Ra E Che xkmas are learned complex-valued projection matrices for each dimension.
Fourier Neural Operator (FNO). For 2D problems, we use a modified version of the original Fourier
Neural Operator (FNO) (Li ct al. 2021 where the output of the integral kernel is passed through a nonlinearity before the residual skip connection, following the implementation al htttts://githe:i
com/ 1wllxx11///// 8/or//rrforrrr retrda/.poraborbbbbo////trrrrrrr 2d_____d__________eeee py] which we found to perform better.
Galerkin Transformer (GKT). We use the Galerkin Transformer (GKT) implementation from https: //github. wttphuuw.ooo///cooooii/o/ cttp.///ccccccaaaa/aaaeiiiiieeee transiiontttaassoomerrrr We employ a multi-input variant where the model is made Mon-Markovian by providing the last K = 4 steps as imput and processing the temporal dimension as normal channels concatenated with other features, rather than as an independent dimension. We use hidden size of 32 with dropout rates of 0.05 in attention layers and 0.025 foodforward layers.
Factformer 2D.
As part of our evaluation, we include the original Factformer model from Li et al. (023a), which is designed for spatiotemporal modeling over 2D spatial domains. While Buitrago el al. 2025) only evaluated the 1D variant, we adopt the same architectural configuration for the 2D case to ensure fair comparison: four attention layers, hidden dimension of 64, and attention heads with total projection dimension of 512. Factformer 2D applies linear attention sequentially along each spatial axis. Given a hidden tensor w E RRS.xS,xH, ; two separate MLPs (MLP, and MT.Py) are applied to compute keys and queries across S₂ and S,, respectively. After averaging over the complementary axis, we obtain gz: ka E SR and qy, ki E S. Values e E SS, are computed through a shared linear projection, and attention is applied first along x and then y. Our implementation follows the original GitHub repository https://github. cww wwt.//aaraataa/aaaaoronaaaa
with minimal changes limited to data format compatibílity.
We also evaluate the 1D version of Factformer, as introduced in Buitrago al. 2025 using the same configuration of four attention layers, hidden dimension 64, and 4 attention heads (total projection dimension 512). Unlike the 2D case, the 1D variant processes inputs over single spatial dimension. As such, only one MLP (MLP,) is used to compute queries and keys, and no spatial averaging is applied. A single linear attention operation is performed per layer along the 1D spatial axis. Values are projected from the input using a linear layer. This model implemented within the same codebase as the 2D version to maintain consistency, with minor modifications to handle 1D inputs.
Factformer 1D.
U-Net Neural Operator (U-Net). Our U-Net implementation follows the typical encoder-decoder architecture with skip connections (Gupta Brandstetter, 2023). It consists of four downsampling con- volutional blocks, bottleneck convolutional block, and four upsampling blocks with residual connections linking corresponding encoder and decoder layers. We use first hidden dimension of 32, and apply chan- nel multiplicrs of [1,2, 4, 8) in the encoder path. No time embeddings are used. This setup is adapted to process spatiotemporal data by flattening the spatial dimensions and independently applying the net- work to each time step. Our implementation is based on the publicly available codebase from PDEBench bittps://gtthe: com/ wwmccoom/ieeeeeeeeeee wwwpeench/pEBeeteeeeeeeeeeeeeeee variant used prior work by Buitrago al. (2025) employs a similar architecture but with channel multipliers of 2,
(Takamoto et al. 2022). For reference, the U-Net Neural Operator
28

---

Published in Transactions on Machine Learning Research (03/2026)
F Data Generation
F.1 1D Burgers' Equation
We consider the one-dimensional viscous Burgers' equation, expressed as:
0,u + u u, v Unaul (1, [0, T [],
where is the viscosity coefficient. For our experiments, utilize the publicly available 1D Burgers' dataset from PDEBench (Takamoto et al. 2022), with viscosity 0.001. The original dataset contains 10,000 spatiotemporal samples, each defined on a spatial grid of 1, 024 points and 201 time steps over the interval [0, 2.01] seconds. We restrict our usage to the first 140 time steps and uniformly downsample them to obtain 20 steps spanning up to 1.4 seconds. This truncation is motivated by the observation that, after this point, the dissipative effect of the diffusion term O,zu suppresses high-Irequency dynamics, causing the solution evolve slowly (Buitrago al. 2025). For training, we use the first 2,048 samples from the dataset and reserve the last 1,000 samples for testing.
F.2 1D Kuramoto Sivassii-- Equation
The -Svamott-Sivashissky (KS) equation given by:
0,u + u ud2u + Oxxl + Uezzz 0, (t, E [0, 7] x [),
with periodic boundary conditions. We follow the exact same setting as provided by Buitrago et al. (2025) using their public repository at httpss///iituu/ http///thuub.ooommbuiitee cowm/rrbuittago/LPSPSS
which builds upon the implementation of Brandstetter al. (2022). The spatial domain [0, 64] discrotized into 512 points, and the temporal domain 00 2.5] is divided into 26 equispaced lime steps with fixed A 0.1.
Initial conditions are generated random superpositions of sine waves:
22rk;
[o(2)
Σ A sin
L
==
where for each trajectory, the amplitudes A, are sampled from continuous uniform distribution [-0.5, 0.5], the frequencies & are drawn from discrete uniform distribution over (1. 8), and the phases ø; are sampled uniformly from [0, 2m).
We consider three different viscosities: v 0.075, 0.1, and 0.125. For each viscosity, we generate 2, 048 training samples and 256 validation samples.
F.3
2D Navier Stokes Equations
equations on the unit torus T = [0, 1]2 in vorticity form. The evolution of the scalar vorticity field w(x y, t) is governed by:
The TorusLi Dataset. We consider the two- dro-dmmensionll incompressible
0, u Ve = vAw
where u u= (u, v) is the velocity field satislying the incompressiiility condition V 0, > 0 is the kinematic viscosity, and f is fixed external forcing function.
We directly reuse the dataset released by Li et al. 2021 referred to as TorusLi, which was originally developed to benchmark the FNO. The simulations were generated using a pseudospectral Crank Nicolson second-order time-stepping scheme on rishutesolution 256 256 grid, and subsequently downsampled to
64 64. All trajectorics use a constant viscosity of v 10-5 (corresponding Reynolds number Re 2000), and share the same external forcing:
f(x, v) = 0.1 (2 v)) cos + v))]
29

---

Published in Transactions on Machine Learning Research (03/2026)
Initial vorticity fields wo are sampled from Gaussian random field:
wo ~ N
with periodic boundary conditions. The numerical solver computes the velocity field by solving a Poisson equation in Fourier space and evaluates nonlinear terms in physical space with dealiasing applied. The nonlinear term is handled explicitly in the Crank Nicolson scheme.
The dataset consists of solutions recorded every == 1 time unit, with total of 20 time steps per trajectory, corresponding to final time horizon of T = 20. This makes the task relatively long-range compared to other PDE benchmarks. The spatial resolution is fixed at 64 x 64 for all experiments in this paper.
The TorusVis and TorusVisForce Datasets. We utilize two additional datasets, TorusVis and TorusVisForce, introduced by Tran et al. (2023), to evaluate model generalization under varying phys- ical regimes. Both datasets are generated using the same TorusLi, maintaining consistency in numerical methodology.
pseudospectral solver used in
These dalasets extend the Navier Stokes setting by incorporating variability in viscosity and external lorcing. Specifically, each trajectory uses a randomly sampled viscosity v between 10 and 10-4. The external forcing function is defined as:
2 E.
f(t,z,y) = 0.1 p=1 ==
[ape sin (2mplix + + + 6t) + Spis cos (2mp(ir + jy) + 8t]],
where Spij ~ 24 1] are sampled independently for cach trajectory. The parameter 8 controls the temporal variation in the forcing: for TorusVis, 8 = 0, resulting in time-invarianl forcing; for TorusVisForce, 0.2, producing a time-varying force.
As with TorusLi, the spatial resolution is fixed at 64 x 64, and trajectorics consist of 20 time steps sampled every / time unit, yielding total time horizon of T 20.
F.4 2D Richtmeyer-Meshkov (CE-RM) Problem
We consider the 2D Richtmeyer-Mestt Meshkov (CE-RM) benchmark for the compressible Euler equations (7 = L.4) on the unit square [0, 1]2 with periodic boundary conditions. The initial state contains a high- pressure circular region and heavy fluid on one side of perturbed interface (Herde et al. 2024). Specifically, the initial pressure and densily are given by
20,
V2²+z² < 0.1,
Po(r, y) =
(11)
1,
otherwise,
2,
|a|| < y,,
po(x,y) =
(12)
1, otherwise,
with initial velocities vs = 0g = 0. The interface T (x, w) perturbed by a random Fourier series:
10
I 1(x, w. = 0.25 + Σ (w) sin(2n((2, v) (1, + (,)))),
(13)
where aj and by are independent uniform random amplitudes and phases (with Σ, aj = 1). We integrate the Euler equations up to time T = 2, saving 21 equally spaced snapshots in time.
The publicly released CE-RM dataset (Herde et al.) 2024) contains 1260 trajectories on a 128 x 128 grid, which we spatially downsample by factor of two to 64 64. Each snapshot includes five fields (density p. velocity ve, Uy, pressure p, and passive tracer). In our setup, the passive tracer is assumed available as an input at every time step, while the model predicts only the conserved variables [p, Uz, By p). We follow the
30

---

Published in Transactions on Machine Learning Research (03/2026)
dataset split of Herde et al. (2024), using the first 1000 trajectories for training and the last 200 for validation. We further apply temporal downsampling by factor of 2: we take the solution at times indexed / = 0, 2, 4, (the first four snapshots) as imput and predict the next 7 steps at indices 8, 10,..... 20. In other words, given the fields at time steps tu, t2, ta, to: the model prodicts the fields at .8,
F.5
2D Gravitational Rayleigh-Taylor (GCE-RT) Problem
The GCE-RT problem adds gravitational forcing to the compressible Euler equations on [0, 1]2 with periodic boundaries. We use the two-dimensional Rayleigh Sayleigh-Taylor Taylor setup from astrophysics: a 7 2 polytropic equilibrium on a model neutron star perturbed at a random interface (Horde et al. 2024) In cylindrical
symmetry (r = vx² + y²), the umperturbed pressure and gravitational potential are
Po pooii(()))) or Po (arr))
p(r) = Ko
(14)
(() = -2Ko
(15)
ar
with Ko = po/p. and a = (2Ko) (with G 1). The initial density profile is
Ko Po sonin(ar)
p(r) = V (r)
(16)
where (r) Ky for reT and ((r) (1 A/(1 + K) for TRT. The Rayleigh- Taylor interface radius is given by
(RT v) 0.25 (1 cos cos 2) + b)).
(17)
with random amplitude e [-1,1] and phase b [- n]. We also perturb the central density Po. pressure Po: and Atwood number A via
Ps 1 0.2c,
(18)
Po = 1 + 0.2d,
(19)
= 0.05 + 0.2e),
(20)
with e, d, e ~ U[11, Initial velocity is set to zero. We evolve this setup to T = 5 and save 11 snapshots
(every A/= = 0.5).
The GCE-RT dataset (Herde al. 2024) likewise contains 1260 trajectories on a 128 128 grid, which we spatially downsample
64. Each snapshot includes six fields (p, Uz, Uy, P. tracer, and the gravitational potential Φ). In our setup, the passive tracer and gravitational potential are assumed available as inputs every time step, while the model predicts only [p, Uxy Uy p). We follow the datasct split of Herde et al. 2024 using the first 1000 trajectories for training and the last 200 for validation. Since the raw data has Lime frames, we take indices 0-3 as input and predict indices 4-10 (i.e. given the first snapshots we predict the next 7).
G Pseudocode of 2D SS-NO
class SS-NO (nn. Module):
Notation:
B: batch size
T: temporal length
Y: spatial dimensions
input channels
hidden dimension
31

---

Published in Transactions on Machine Learning Research (03/2026)
def __init__ (self,
lifting_layer: nn.Module,
projection_layer: nn. Module,
nenory-pre_forward.xx: Module,
nemory_pre_bbckward.. nn. Module,
nemory_pre_forward..y nn. Module.
nemory_pre_backward..y nn. Module,
memory_t: mo
nemory_post_formerd..x nn. Module,
memory_post__ackwarr_.x nn on.Medl
nemory-post_fooward. nn. Module,
semory_.oost_bakwarddd nn :odulee
(uper()...iiit...3
self =iftiing_lly
=
projection_layer
self.remery.pre_.ooreee.l
remory_pre_forward.
self.renrry.pre..ccewrr. nemory-pre.backward..
saf..remcry_rreefoorrrddy nemory_pre_forward..
sel..eeoory.ree.accarrry nemory_preebackwarrd.y
self.nemory_
eelf.remorryppsttforrarr. nemory.post_forward..x
celf.nemmrryypsssbaawwrr.r... menoryy_post.bacwerd..
self sea.f.meoorypoosfforwrrrr memory_post_forwardy
self meeorrypoott_acckeeeeee nenoryypost_bacckard_.
def forlard(self, x: Tensor) >> Tensor:
Args:
x: Input sequence of states (B, C, X, Y, T)
Returns:
Predicted next state (B, X, Y, 1)
if telfitraining
==x torch. ranen_like(x) noise
rearrange c x)yc) self.p(x)
x rearrange 'bb ()
hbt))xy, t=T)
Pre-menory spatial SSM x rearrange '(b t) h x -> (b x', t=T)
x_ fwd wrw.fmemory_peeeeee. m (xx
x_ bwd torch. flip wew.ry..re..ccwwrr
flip(x, dims=(-1))) dims=(-1])
x = x_fwd x_bwd
x = rearrange (x, ((b y) h (b
t=T)
x = rearrange (x, '(b x) >> (b h y', (=T)
f_f self.memory memory_oree y_bwd torch. flip menory.ore.baccwardyyooo
dims=(-1])) dims=[-1])
= y_fwd y_bwd rearrange(x, tbb > hb t=T)
Temporal SSM
rearrange (x, (bb t) x -bx t', t=T)
x rearrange (x, *b y h xb t²)
x self. memory_t(x)(e))
x rearrange (x, '(b y) h b t', x=X, y=y) t=T)
x rearrange (x, "b y h t >> (b t)
# Post-memory spatial SSM
x = rearrange '(b t) h -> (b h x', t=T)
x_fwd = self. melf.memory_post. www-posttforward_x(177
x x_bwd = torch. flip(self. menory..post.backwrdd (torch. fxip(x, diss= dims=[-1])
x = x_fwd x_d
x = rearrange(x. (bb y) h > <bbbbbb (b()x (=T)
rearrange (x, '(b t) x > (b t=T)
32

---

Published in Transactions on Machine Learning Research (03/2026)
y_fwd self. peemory.oostfforwrr.yyyoto
y_bwd torch. flip .emorry.pstt (torch. fxip dins=[-1)))(0,, dims=[-11)
y_fwd y_bwd
rearrange(x, ((b x) (b tb t=T)
x self.q(x)
x rearrange '(b t) yxyy---- bbxytc' t=T)
return x
Listing 3: SS-NO pseudocode
Computational Efficiency and Cost-Accuracy Analysis
To rigorously evaluate the practical deployment viability of the proposed method, we conducted a comprehen- sive analysis of the trade-offs between predictive accuracy, theoretical complexity, and real-world inference latency.
H.1 Benchmarking Methodology
All benchmarks were performed on a single NVIDIA A6000 GPU with batch size of to simulate real-time inference conditions. We report two key metrics across spatial resolutions N E (32, 64, 128, 256, 512}:
Theoretical Complexity (GFLOPs): The total number floating-pon point operations (in billions) required for a single forward pass, calculated using the thop profiler with custom handlers for spectral operations.
Inference Latency (Time): The wall-clock time in milliseconds, averaged over 100 iterations following a GPU warmup phase.
H.2 Results and Discussion
The results are summarized in Table 6 and visualized as Pareto frontiers in Figure 12
State-Space Neural Operator (SS-NO). Our model demonstrates the most favorable cost-accuracy profile. achieves the lowest inference latency across all resolutions, consistently clocking under 0.85 ms. Notably, its runtime is effectively constant with respect to resolution in this regime, as the computation is dominated by fixed kernel launch overheads and the highly parallelizable associative scan of the S4 backbone.
Galerkin Transformer (GKT). We acknowledge the theoretical strength of the Galerkin Transformer, which exhibits the lowest GFLOP count of all models (0.03 GFLOPs at N=512). This efficiency stems from its linear attention mechanism, which avoids the quadratic complexity of standard Transformers. However, theoretical efficiency does not translate directly to wall-clock speed; GKT's inference latency (=2.26 ms) is nearly 3x higher than SS-NO due the memory-bound nature of attention layers and the overhead of frequent reshaping operations.
FFNO and U-Net. The FFNO remains strong competitor, offering low theoretical complexity (0.12 GFLOPs at N=512) and respectable runtime (<1.39 ms). However, SS-NO outperforms it in both speed (x1.7 faster) and accuracy. The U-Net and FactFormer architectures are significantly more expensive, with the U-Net requiring an order of magnitude more FLOPs (1.41 GFLOPs) for comparable or worse accuracy, highlighting the inefficiency of dense convolutional encoders for this class of problems.
In summary, while models like GKT and FFNO offer specific theoretical advantages, SS-NO achieves the optimal practical balance, delivering the highest accuracy with the lowest real-world latency.
33

---

Published in Transactions on Machine Learning Research (03/2026)
Complexity vs.
Latency
Relative Error (Log Scale)
Relative I, Error (Lag Scale)
10-2
10-3
Computational Complexity (GFLOPs)
Latency
SS-NO
FANO
FactFormer
U-Net
GKT
Figure 12: Cost- Accuracy Trade-off Analysis (ID KS, v 0.075). We compare predictive error against. computational cost across resolutions N 512)...,11). Left: Theoretical Complexity (GFLOPs). SS-NO (Red) dominates the Pareto frontier, achieving the lowest error with significantly fewer operations than FactFormer and U-Net. Right: Empirical Inference Latency (ms). On NVIDIA A6000 GPU, SS-NO is the fastest model (< 0.85 ms) and maintains constant runtime scaling. Note the discrepancy for GKT (Purple): despite low theoretical FLOPs, it suffers higher real-world latency due to attention overheads.
Table 6: Computational Cost Summary. Inference latency (ms) and theoretical complexity (GFL.Ps) for all models across the full rango of spatial resolutions (N = to 512). SS-NO consistently demonstrates the lowest inference latency (+s 0.81 ms) and competitive scaling. While GKT shows very low theoretical FLOPs, its real-world latency is significantly higher due to attention overheads.
N=32
N=64
N=128
N=256
N=512
Model
Time
GFLOPs
Time
GFLOPs
Time
GFLOPs
Time
GFLOPs
Time (ms)
GFLOPs
(ms)
(ms)
(ms)
(ms)
U-Net
1.29
0.088
1.24
0.176
1.22
0.352
1.25
0.704
1.24
1.407
GKT
2.30
0.002
2.26
0.004
2.26
0.007
2.25
0.015
2.26
0.030
FactFormer
3.36
0.165
3.39
0.330
3.38
0.659
3.43
1.319
3.37
2,638
FFNO
1.37
0.016
1.35
0.033
1.35
0.044
1.35
0.067
1.39
0.120
SS-NO (Ours)
0.81
0.026
0.82
0.052
0.81
0.104
0.209
Formal Proof of FNO Recovery by S5-NO
In this section, we rigorously prove that the State Space Neural Operator (SS-NO) kernel formulation structurally subsumes the Fourier Neural Operator (FNO) kernel a special case. Specifically, we show that by fixing the SS-NO parameters to specific values (zero damping, harmonic frequencies), the effective spatial convolution kernel of SS-NO becomes mathematically identical to the implicit spatial kernel of the FNO defined via the Discrete Fourier Transform (DFT).
1.1 The FNO Kernel Formulation
The Fourier Neural Operator (FNO) operates by computing convolution in the spectral domain. For 1D input function u(z) discretized on grid of N points {In N J the FNO computes the output v as:
N-1
-= (R (R-
where R is a learnable weight matrix in the frequency domain. By the Convolution Theorem, this is equivalent to circular spatial convolution KFNO, where the spatial kernel RENO the IDFT of the spectral
34

---

Published in Transactions on Machine Learning Research (03/2026)
weights R. Let Rg denote the spectral weight matrix for the k-th frequency mode. The spatial kernel RENO evaluated at a grid point is:
N-1 Rg tg
KFNO
(Equation A)
k=0
1.2 The SS-NO Kernel Formulation
The SS-NO applies convolution with continuous kernel derived from the state equation. As detailed in Equation of the main text, the effective spatial kernel is:
RSSM (x)
Σ mm
Cm Bue pm gumee
Evaluating this on the discrete grid In
(assuming forward direction):
N
M
KSSM(n] = Σ
(C Dm)e Pm dum
(Equation B)
m=1
1.3 Proof of Equivalence
We demonstrate that Equation B (SS-NO) can exactly recover Equation A (FNO) under the following
restrictions: Zero Damping (p = 0): Set Pm = for all stales, eliminating the decay term (e = 1). Harmonic Frequencies (w 2rk): Set the continuous frequency parameter Wm 2nk. On the grid, the exponent becomes (2(2m) & e Rukn recovering the DFT basis. Spectral Weights (SVD Decomposition): Unlike FNO, which uses a full-rank matrix R2 for channel mixing, single SSM state provides a rank-1 mixing term Cm However, any matrix Rş can be decomposed into a sum of J rank-1 terms (e.g., via Singular Value Decomposition): Rg = EE= Bx To recover FNO, we set the number of SSM states to M N J. We assign states to each frequency mode k indexed by j. By initializing the parameters such that Σÿ-1 B. Rs, Equation becomes:
N-I
N-1
RSSM(N
Σ A==
Ckj
d cei
Rgd
2
k=
Conclusion. This mathematically identical to Equation A. Thus, the SS-NO formulation strictly subsumes the FNO kernel parametrization. SS-NO is a generalization because it allows for p > 0 (localized kernels), w # 2nk (adaptive frequencies), and low-rank approximations where fewer than states are used per frequency.
J Extensions to Irregular Geometries and Complex Unknown Systems
A primary motivation for Neural Operators their ability to model systems where the governing PDEs are cither unknown, partially known, or computationally intractable solve directly. In such regimes, the operator must learn the underlying physics purely from observational data. In this section, demonstrate SS-NO's capability in these challenging scenarios by extending it to irregular domains via coordinate transformation and slicing.
J.1 Symmetry, Equivariance, and Geometric Inductive Bias
We first explicitly acknowledge a limitation raised by reviewers: the factorized scanning mechanism employed in main text (processing dimensions sequentially, e.g., & v) is not inherently equivariant to rotation
35

---

Published in Transactions on Machine Learning Research (03/2026)
reflection. Unlike the isotropic kernels of standard CNNs the global spectral bases of FNOs, the SSM scan introduces a directional inductive bias. For example, rotating fluid field by 90 changes the order in which the SSM encounters features, potentially altering the output. Similarly, for periodic boundary conditions, the standard causal scan (t ++ 1) breaks translation symmetry unless specific circular padding schemee are used.
To reconcile this, we explore combining SS-NO with existing coordinate transformation frameworks, specifically utilizing the diffcomorphic mappings introduced in Geo-FNO (Li et al. 2023b). By decoupling the scanning grid from the physical domain, we can maintain the computational efficiency of the factorized scan while respecting the geometric properties of the data.
J.2 Approach 1: Recovering Symmetry via Transformations (Shallow Water Equations)
To validate SS-NO on system with complex topology and no explicit PDE supervision, we test on the spherical Shallow Water Equations (SWE).
To address the symmetry breaking discussed above, we adopt the coordinate transformation approach from Goo-FNO (Li et al. 2023b). This method learns a diffcomorphism that maps the physical manifold (the sphere) to uniform computational latent grid. Crucially, because this mapping is diffeomorphic, it can preserve spatial symmetries such as rotation and reflection- in the physical space, even if the internal computational scan remains sequential. The transformation layer effectively "unwraps" the domain into a canonical representation where the directional bias of the SSM is mitigated.
Experimental results against the Spherical FNO (SFNO) suggest that SS-NO's adaptive frequency learning is uniquely suited for this transformed space. In a mapped domain, physical waves often manifest as warped, non-harmonic oscillations that fixed Fourier bases struggle to represent sparsely. SS-NO automatically adapts its continuous frequency parameters wm to these distorted geometries, effectively learning physics-compliant basis from data.
To rigorously test whether this combination has learned the underlying operator, we utilized a strict evaluation protocol: models were trained on a short horizon (10 steps) but evaluated on stability during long rollouts up to t = 100. As shown in Table 7 SS-NO demonstrates superior stability, achieving ~ 2.2x lower error than the baseline at l == 100, confirming that the learned diffeomorphism successfully preserved the necessary global dynamics.
J.3 Approach 2: Complex Aerodynamics via Slicing (AirFRANS)
We further evaluate SS-NO on the AirfRANS dataset (Bonnet et al., 2022), representing a "partially known" or complex system where analytic solutions unavailable. This task involves predicting velocity, pressure, and skin friction fields over NACA airfoils. Here, the "true" PDE (Navier-Stokes) is computationally prohibitive for design optimization, and the model must act as surrogate learned entirely from RANS simulation data.
For this irregular mesh problem, we adopt the slasing-based approach from Transolver (Wu et al. 2024). We integrated the SS-NO block into the Transolver architecture, replacing its core attention mechanism. As shown in Table 8 SS-NO achieves competitive regression performance, significantly lowering Volume Error (0.0017 vs 0.0037) and Surface Error (0.0046 vs 0.0142) compared to the baseline.
However, the model lags slightly in derived physical cocfficients (Lift Cz). We hypothesize this is due to the Linear Time-Invarian (LTT) nature of the standard S4 backbone, which applies uniform dynamics globally. Complex aerodynamics require spatially varying dynamics (e.g., distinguishing laminar leading edges from turbulent wakes). We anticipate that replacing the S4 backbone with input-dependent State Space Models (e.g., Manba (Gu & Dao 2024)), which can dynamically modulate parameters based on the local flow state, would bridge this gap.
36

---

Published in Transactions on Machine Learning Research (03/2026)
Table 7: Spherical Shallow Water Equations (SWE) Benchmark. Comparison of relative & errors on spherical SWE forecasting over long time horizons. Models were trained on a short horizon (10 steps) and evaluated up to 100. SS-NO demonstrates superior stability and generalization, achieving 2.2x lower error than the Spherical FNO (SFNO) at t 100, indicating it has effectively learned the global atmospheric dynamics from data.
Evaluation Step
SFNO (Baseline)
SS-NO (Ours)
Improvement
= 20
6.56 x 10-2
2.63 x 10-2
2.5x
= 50
5.99 x 10-1
2.54 x 10
2.4x
100
2.82 10
1.28 10
2.2x
Table 8: AirRANS Benchmark (Data-Driven Aerodynamics). Comparison of SS-NO against Transolver on the AirfRANS dataset. This task represents "partially known" system where the model must learn complex acrodynamic mappings purcly from data. SS-NO excels at field reconstruction error) but requires input-dependent gating (e.g., Mamba) to better capture global coefficienss like Lift.
Surracc
Model
Error Vol. (4)
Error Surf. (4)
Lift Coef. CL (1)
Spearman PL (t)
Transolver
0.0037
0,0142
0.1030
0.9978
SS-NO (Ours)
0.0017
0.0046
0.1098
0.9973
K
Background: The MemNO Framework
In this work, we adopt the training and inference strategy proposed the Memory Augmented Neural Operator (McmNO) framework Buitrago al. (2025). As this framework serves as the temporal backbone for our experiments, we provide sunmary of its theoretical motivation and architectural design here.
K.1 Theoretical Motivation: Mori-Zwanziin and Coarse-Graining
The primary motivation for MemNO stems from the Mori-Zwanzig formalism in statistical mechanics. This theory addresses the problem of modeling a dynamical system when the full state is not observable Anservable-seccificlllll when the system is projected onto a lower-dimensional or coarser representation (e.g., discretizing a PDE onlo a coarse spatial grid).
Even if the underlying PDE (the microscopic dynamics) is Markovian the future depends only the current exact state-the course-grained dynamics generally non-Markovian. The resolved (observable) scales interact with the unresolved (sub-grid) scales, and the influence of these unresol scales manifests as a "memory" of the system's history. Mathematically, the evolution of the coarse variables depends on an integral over the past states (the memory kernel).
Standard autoregressive Neural Operators often ignore this, treating the coarse mapping as Markovian (de+1 = (u.)). NO argues that to accurately model coarse-grained PDE trajectories, the architecture must explicitly account for this memory term to compensate for the loss of information due to spatial discretization.
K.2 Architectural Implementation: S4 Temporal Layers
MemNO employs $4 layers along the temporal dimension. The temporal S4 mechanism allows the model efficiently compress the entire history of the input function into compact latent state. This enables the model to effectively learn the non-Markovian memory kernel required close the system dynamics, essentially using temporal information to recover spatial details lost during coarse-graining.
37

---

Published in Transactions on Machine Learning Research (03/2026)
K.3 Adoption in This Work
We adopted the MemNO framework for our temporal benchmarks for two key reasons:
1. Handling Coarse Resolutions: Many of our benchmarks (e.g., Navier-Stokes at 64 x 64) operate in regimes where the spatial grid does not fully resolve the turbulence. The McmNO framework allows all models to utilize temporal history to compensate for this aliasing, ensuring a physically rigorous evaluation.
Standardized Temporal Backbone: By fixing the temporal architecture to use MemNO's S4- based time mixing for all baselines (SS-NO, FNO, FactFrrmer, ctc.), we ensure a fair comparison. This isolates the contribution of the spatial mixing architectures, confirming that the performance gains reported for SS-NO arise from its superior spatial processing rather than advantageous temporal integrator.
L
Limitations
Despite the significant reduction in model size and strong performance achieved by SS-NO, certain limitations remain inherent to the current formulation.
First, the accumulation of gradients during autoregressive training poses a scalability bottleneck. While SS-NO is parameter-efficient,, increasing the state dimension D or the sequence length L linearly increases the memory required for backpropagation through time. This restricts the fcasibility of adopting extremely large latent states without specialized techniques. However, unlike language modeling where infnite-horizon dependencies are crucial, PDE dynamics are often effectively Markovian, suggesting that techniques such as truncated backpropagation or gradiont checkpointing could provide effective mitigation.
Second, the sequential scanning mechanism introduces a directional inductive bias that not inherently equivariant. Unlike isotropic convolutions or global spectral layers, the state-space scan processes data sequentially (e.g., + y), which breaks rotational and reflection symmetries. While we have shown in Appendix
g that this can be managed via symmetric embedding layers coordinate transformations, the core mixing block itself does not guarantee equivariance by construction.
M Possible Extensions and Future Work
We see several promising directions to extend the SS-NO framework, particularly in enhancing its expressivity for complex dynamics and further generalizing its geometric capabilities.
Input-Dependent Backbones (Mamba and Hydra). A distinct advantage of the SS-NO framework is its modularity; the core scanning backbone can be swapped for more advanced architectures. While the current implementation utilizes efficient Linear Time-Invariant (LTT) systems, a natural extension incorporate input-dependent backbones such Mamba (Gu & Dao 2024) or Hydra (Hwang al. 2024). These models allow state transition matrices to be functions of the input field, enabling the operator to dynamically modulate its effective receptive field and frequency response based on local flow features (e.g., shock waves or boundary layers). This would be particularly beneficial for multi-regime systems like the AirfRANS benchmark discussed in Appendix J
Native Graph and Manifold Learning. While successfully demonstrated extending SS-NO to irregular domains via slicing and coordinate transformation, these methods essentially project data onto regular latent structures. A more fundamental approach is to generalize the scanning operator to work natively on graphs or meshes. Recent works like Graph Mamba (Wang et al. 2024) and pLSTM (Pöppel ct al. 2025) define propagation orders over graph topologies. Adopting these backbones would allow SS-NO to process unstructured meshes without the need for interpolation or slicing, offering a principled way handle complex engineering geometries.
38

---

Published in Transactions on Machine Learning Research (03/2026)
True Geometric Recurrence. A more theoretical direction is to redefine the recurrence relation itself using geometric differential operators. This would entail replacing the ordinary time derivative in the SSM with a covariant derivative,
Va h = Ah + Bu,
incorporating notions of parallel transport evolve the hidden state along geodesics. Such an extension would enable genuinely Riemannian state evolution, preserving vector field symmetries on curved manifolds by construction, though it introduces significant mathematical complexity.
Memory-Efficient Training Strategies. Finally, improving the memory efficiency of autoregressive training remains a priority. Given that turbulent flows often exhibit chaotic but bounded attractors, exact gradient retention over long horizons may yield diminishing returns. We plan to investigate truncated backpropagation and implicit differentiation strategics specifically optimized for Neural Operators, aiming to balance the fidelity of long-tern gradients with the computational constraints of high-resolution 3D simulations.
39


Reviews:
[{"review_id": "JdkcfBdlHQ", "reviewer": "Reviewer_V25i", "summary": "", "questions": "", "limitations": "", "rating": null, "confidence": null, "soundness": null, "presentation": null, "contribution": null, "strengths": "", "weaknesses": "", "quality": null, "clarity": null, "significance": null, "originality": null, "strengths_and_weaknesses": "", "was_revised": false, "final_justification": "", "rebuttal": "", "claims_and_evidence": "Yes", "methods_and_evaluation_criteria": "", "theoretical_claims": "", "experimental_designs_or_analyses": "", "relation_to_broader_scientific_literature": "", "essential_references_not_discussed": "", "other_strengths_and_weaknesses": "", "other_comments_or_suggestions": "", "questions_for_authors": "", "summary_of_paper": "", "summary_of_recommendation": "", "technical_quality": "", "clarity_of_presentation": "", "potential_impact": "", "robotics_focus": "", "quality_of_limitations": "", "questions_for_rebuttal": "", "reasons_to_accept": "", "reasons_to_reject": "", "questions_to_authors": "", "paper_topic_and_main_contributions": "", "excitement": null, "reproducibility": null, "ethical_concerns": "", "main_review": "", "summary_of_the_paper": "", "summary_of_the_review": "", "strength_and_weaknesses": "", "limitations_and_societal_impact": "", "review_text": "", "summary_of_contributions": "The main contributions of this work are:\n1. Connects state space models to neural operators, and in particular, FNOs. It is shown that in 1d a SSM can exactly recover a FNO that can also modulate the frequency components, giving it the ability to more easily localize.\n2. Proposes a factorized SSM by applying a 1d SSM in each direction similarly to a factorized transform like the FFT.\n3. Proves a general universal approximation theorem for integral operators of the convolution type, showing that all that is needed is for the composition of all the kernels in each layer to have global support.\n4. Extensive numerical experiments on Burger's equation, the Kuramoto–Sivashinsky equation, and the Navier-Stokes equation show the method is promising and outperforms 4 to 5 baselines.", "claims_and_evidence_yn": "Yes", "claims_explanation": "The theory is well-formulated and the poofs given in the appendix are fully rigorous. The numerical experiments are conducted over 5 random initialization of the model weights and show statistically significant improvements.", "audience_yn": "Yes", "audience_explanation": "The paper is significant for the field of data-driven surrogates models for PDEs, connecting promising methods in NLP to operator learning.  This is a significant audience for TMLR.", "requested_changes": "1. In 1d, the resolution dependence of all the models can be tested more rigorously. The maximum resolution of 128 is quite small and doesn't demonstrate the resolution dependence of the U-net likely because the architecture tested still has a fixed field of view at that resolution. It would be helpful to put these results on a plot instead of in tables as it makes comparisons much easier.\n2. The run-time cost of any of the models is not consider. Parameter count is not really good proxy for this. Cost-accuracy trade-off curves would be very useful, especially as the resolution grows (in 2d and in 1d). This can be done w.r.t. a theoretical calculation for the number of flops for each method or the run-time on the same hardware (or both). \n3. The authors claim in the text that the kernel having a \"full field of view\" is a necessary and sufficient condition for universality but the theorem is stated as the condition only being sufficient. Could the authors clarify and provide extra discussion if condition is indeed not necessary.", "broader_impact_concerns": "There are no ethical implications.", "additional_comments": "", "extra_scores": {}, "extra_text": {}}, {"review_id": "nT1BUfmL4w", "reviewer": "Reviewer_Dod9", "summary": "", "questions": "", "limitations": "", "rating": null, "confidence": null, "soundness": null, "presentation": null, "contribution": null, "strengths": "", "weaknesses": "", "quality": null, "clarity": null, "significance": null, "originality": null, "strengths_and_weaknesses": "", "was_revised": false, "final_justification": "", "rebuttal": "", "claims_and_evidence": "No", "methods_and_evaluation_criteria": "", "theoretical_claims": "", "experimental_designs_or_analyses": "", "relation_to_broader_scientific_literature": "", "essential_references_not_discussed": "", "other_strengths_and_weaknesses": "", "other_comments_or_suggestions": "", "questions_for_authors": "", "summary_of_paper": "", "summary_of_recommendation": "", "technical_quality": "", "clarity_of_presentation": "", "potential_impact": "", "robotics_focus": "", "quality_of_limitations": "", "questions_for_rebuttal": "", "reasons_to_accept": "", "reasons_to_reject": "", "questions_to_authors": "", "paper_topic_and_main_contributions": "", "excitement": null, "reproducibility": null, "ethical_concerns": "", "main_review": "", "summary_of_the_paper": "", "summary_of_the_review": "", "strength_and_weaknesses": "", "limitations_and_societal_impact": "", "review_text": "", "summary_of_contributions": "This paper introduces the State Space Neural Operator (SS-NO), which extends SSMs to spatiotemporal operator learning for time-dependent PDEs. The main contributions are: 1. a bidirectional spatial SSM architecture that applies S4D-like convolutions along each spatial dimension sequentially, 2. a learnable damping coefficients and frequency modulation that provide adaptive receptive fields, and 3. a universality theorem showing that convolutional neural operators with \"full field of view\" are universal approximators.\n\n**strengths**\n1. The parameter efficiency is well notable.\n2. The ablation study on damping mechanisms is thorough and provides good insight into when/why damping helps\n3. The theoretical connection between SSMs and FNOs through the kernel formulation (Eq. 4-5) is interesting\n\n**weaknesses**\n1. The claim that SS-NO \"subsumes\" F-FNO is not clearly demonstrated\n2. The universality proof, while technically sound, has unclear practical implications\n3. Breaking of spatial symmetries (rotation, reflection, translation) is not discussed\n4. Benchmarks are relatively easy and do not test settings where neural operators are most useful (PDEs not known)", "claims_and_evidence_yn": "No", "claims_explanation": "The experimental results are in general well-presented and the ablation studies are convincing for the claims about damping. However, several theoretical claims are not adequately supported:\n\n1. The paper claims SS-NO \"subsumes\" F-FNO as a special case, but the actual mechanism for recovering F-FNO kernels is hand-wavy. Setting $\\rho_k = 0$ and $\\omega_k = k$ in Eq. 5 does not obviously recover the F-FNO parameterization since the kernel structure differs (exponential decay vs truncated Fourier series).\n\n2. The universality theorem (Theorem 4.1) is proven rigorously in the appendix, but its purpose is unclear. Both FNO and SS-NO trivially satisfy the full field-of-view condition, so what does this result add? The paper states it \"closes a gap\" for factorized architectures, but doesn't explain why this gap mattered or what practical insight we gain.\n\n3. The claim about \"enhanced model expressivity\" over F-FNO (Section 4.2) is theoretical but the empirical comparison doesn't isolate this effect. SS-NO also has different architecture choices beyond the kernel parameterization.", "audience_yn": "Yes", "audience_explanation": "The neural operator community may find this work relevant, particularly those interested in efficient architectures for PDE learning. The connection between SSMs and neural operators is timely given the recent success of SSMs in sequence modeling. The practical parameter efficiency gains are meaningful for deployment scenarios. The damping analysis also provides useful insights that could inform future architecture design.", "requested_changes": "**critical**\n1. Clarify or rephrase the claim that SS-NO subsumes F-FNO. Currently it reads as if one can exactly recover F-FNO, but the kernel parameterizations are clearly different. Either provide a rigorous proof of this subsumption or reframe as \"inspired by\" or \"related to.\"\n2. Add discussion on spatial symmetry breaking. The sequential application of 1D SSMs along x then y (and forward and backward) breaks rotational and reflection symmetries (and also translation symmetries for periodic boundaries). For many physical systems these are important inductive biases. The authors should discuss: (a) how this affects results on rotationally/reflectively/translationally symmetric problems, (b) whether data augmentation helps, (c) comparison to architectures that preserve symmetry.\n3. Clarify the purpose of the universality theorem. As stated, it's not clear what new insight this provides beyond confirming that both FNO and SS-NO are universal. If the main point is about factorized architectures, please explain why this was previously uncertain.\n4. Include experiments on problems with unknown or partially known PDEs. This is arguably the most important use case for neural operators. Showing that SS-NO works when we only have data (no explicit PDE) would strengthen the paper considerably.\n5. Test on more challenging 2D/3D problems. The 1D benchmarks dominate the paper but are relatively simple. Higher-dimensional turbulence or multi-physics problems would be more convincing.\n6. Fix grammatical issues, e.g., \"the formulation equation 6\" should be \"the formulation in Equation 6\" or similar. Also check consistency in notation.\n\n**strengthen**\n1. The MemNO framework is referenced repeatedly but the connection could be explained more clearly for readers unfamiliar with that work.\n2. Some baseline comparisons may be unfair. FNO2D with 67M parameters is clearly overparameterized for these tasks. A capacity-matched comparison would be informative.", "broader_impact_concerns": "None", "additional_comments": "", "extra_scores": {}, "extra_text": {}}, {"review_id": "uufKVkNWhE", "reviewer": "Reviewer_JDKn", "summary": "", "questions": "", "limitations": "", "rating": null, "confidence": null, "soundness": null, "presentation": null, "contribution": null, "strengths": "", "weaknesses": "", "quality": null, "clarity": null, "significance": null, "originality": null, "strengths_and_weaknesses": "", "was_revised": false, "final_justification": "", "rebuttal": "", "claims_and_evidence": "Yes", "methods_and_evaluation_criteria": "", "theoretical_claims": "", "experimental_designs_or_analyses": "", "relation_to_broader_scientific_literature": "", "essential_references_not_discussed": "", "other_strengths_and_weaknesses": "", "other_comments_or_suggestions": "", "questions_for_authors": "", "summary_of_paper": "", "summary_of_recommendation": "", "technical_quality": "", "clarity_of_presentation": "", "potential_impact": "", "robotics_focus": "", "quality_of_limitations": "", "questions_for_rebuttal": "", "reasons_to_accept": "", "reasons_to_reject": "", "questions_to_authors": "", "paper_topic_and_main_contributions": "", "excitement": null, "reproducibility": null, "ethical_concerns": "", "main_review": "", "summary_of_the_paper": "", "summary_of_the_review": "", "strength_and_weaknesses": "", "limitations_and_societal_impact": "", "review_text": "", "summary_of_contributions": "The paper introduces the State Space Neural Operator (SS-NO), a compact and efficient neural architecture that unifies spatiotemporal operator learning for time-dependent partial differential equations (PDEs) using structured state space models. The SS-NO model merges insights from neural operators and structured state space models (SSMs) to learn mappings between function spaces defined by PDEs. Traditional models like the Fourier Neural Operator (FNO) capture global dependencies but suffer from high memory costs, while state space models like S4 efficiently model long-range temporal patterns. SS-NO generalizes FNO by incorporating adaptive damping (which localizes the receptive field for stability) and learnable frequency modulation (which selects spectral modes dynamically), enabling efficient and expressive spatiotemporal modeling. Architecturally, SS-NO employs bidirectional 1D spatial SSMs along each axis and integrates a temporal memory module following the MemNO framework. A theoretical universality theorem is proven, showing that convolutional neural operators with full field-of-view—including SS-NO—can approximate any continuous operator. The model is evaluated on several 1D and 2D PDE benchmarks, including Burgers’, Kuramoto–Sivashinsky, Navier–Stokes, and compressible Euler equations. SS-NO outperforms existing methods in accuracy while using significantly fewer parameters. Extensive ablation studies analyze the role of damping, frequency modulation, temporal memory, and architectural variants, revealing that SS-NO’s adaptive mechanisms are critical for performance, especially in chaotic or low-resolution regimes.", "claims_and_evidence_yn": "Yes", "claims_explanation": "SS-NO is reportedin the paper to consistently achieve state-of-the-art accuracy across 1D and 2D PDE benchmarks while reducing parameter count by up to 180× compared to baselines. It excels particularly in chaotic or low-resolution scenarios, where adaptive damping and spectral modulation enhance stability and generalization, validating both its theoretical expressivity and empirical robustness. The paper provides a sharp universality criterion for convolutional neural operators with full field-of-view, filling a theoretical gap in the operator learning literature. The paper also carefully isolates the impact of each architectural component (damping, frequency learning, temporal memory) across multiple regimes and capacities. Finally, the benchmarks include diverse dynamical systems with varying dimensionalities, physics regimes, and complexities, making the conclusions broadly relevant.\nThe paper also has some weaknesses. The method assumes regular spatiotemporal grids, with no extensions to unstructured meshes or irregular geometries. Also, training remains memory-intensive due to backpropagation through long sequences, especially when increasing state dimensionality. Finally, while adaptivity to missing context is shown, the role of SS-NO under real-world noisy or multimodal inputs (e.g., from sensors or images) could be further clarified.", "audience_yn": "Yes", "audience_explanation": "Many researchers in the TMLR audience would find modeling state space models for time-dependent partial differential equations (PDEs) useful.", "requested_changes": "I would suggest the authors to adress some of the comments related to how the paper could be improved above. One specific suggestion concerns extending SS-NO to irregular domains or unstructured meshes. Is it possible to reformulate SS-NO on graph domains using graph-based state-space models (e.g., graph SSMs or message-passing variants of S4). Alternatively, the analysis could be applied to extend the universality theorem by defining “full field-of-view” over geodesic or graph distances on manifolds or meshes, enabling theoretical guarantees beyond Euclidean grids.. This would allow SS-NO to operate over unstructured spatial inputs. If so, that could strengthen the model to irregular domains or unstructured scenarios. If a theoretical approach is not possible an empirical meta-analysis could help. For instance, evaluate SS-NO variants on PDEs defined over triangulated meshes (e.g., FEM or CFD simulations on irregular domains), or combine SS-NO with graph neural operators (e.g., MGNO, Diffusion GNNs) and test interpolation/generalization on mesh-refinement or boundary deformation tasks.", "broader_impact_concerns": "NA", "additional_comments": "The paper is well written, so I have no additional comments, except for a minor typo:\nThe paper used both 'field-of-view' and 'field of view'. Choose one form consistently.", "extra_scores": {}, "extra_text": {}}]

Author rebuttal:

Decision: Accept (poster)