Hierarchical Implicit Neural Emulators
Authors: Ruoxi Jiang
Xiao Zhang
Karan Jakhar
Peter Y. Lu
Pedram Hassanzadeh
Michael Maire
Rebecca Willett
Venue: NeurIPS 2025 poster ()
Keywords: Neural emulators; Long-term stability
Primary area: machine_learning_for_sciences

Abstract: Neural PDE solvers offer a powerful tool for modeling complex dynamical systems, but often struggle with error accumulation over long time horizons and maintaining stability and physical consistency. We introduce a multiscale implicit neural emulator that enhances long-term prediction accuracy by conditioning on a hierarchy of lower-dimensional future state representations. Drawing inspiration from the stability properties of numerical implicit time-stepping methods, our approach leverages predictions several steps ahead in time at increasing compression rates for next-timestep refinements. By actively adjusting the temporal downsampling ratios, our design enables the model to capture dynamics across multiple granularities and enforce long-range temporal coherence. Experiments on turbulent fluid dynamics show that our method achieves high short-term accuracy and produces long-term stable forecasts, significantly outperforming autoregressive baselines while adding minimal computational overhead.

Hierarchical Implicit Neural Emulators
Ruoxi Jiang'.6 Xiao Zhang' Karan Jakhar Peter Y. Lu
Pedram Hassanzadeh Michael Maire' Rebecca Willett 1.2
! DDepartment of Computer Science, The University of Chicago
Department of Statistics, The University of Chicago
Department of Geophysical Sciences, The University of Chicago
'Department of Physics, The University of Chicago
Artificial Intelligence Innovation and Incubation Institute, Fudan University Department of Mechanical Engineering, Rice University Department of Electrical and Computer Engineering, Tufts University
Abstract
Neural PDE solvers offer a powerful tool for modeling complex dynamical systems, but often struggle with error accumulation over long time horizons and maintaining stability and physical consistency. We introduce multiscale implicit neural emula- tor that enhances long-term prediction accuracy by conditioning on a hierarchy lower-dimensional future state representations. Inspired by the stability properties of numerical implicit time-stepping methods, we developed an approach that lever- ages predictions several steps ahead in time at increasing compression rates for next-timestep refinements. By actively adjusting the temporal downsampling ratios, our design enables the model to capture dynamics across multiple granularities and enforce long-range temporal coherence. Experiments on turbulent fluid dynamics show that method achieves high short-term accuracy and produces long-tenn long-term stable forecasts, significantly outperforming non non-hierarchical autoregressive base- lines while adding minimal computational overhead. The codebase is available this link!
Introduction
Machine learning (ML)-based surrogate modeling of dynamical systems has spurred great interest in recent years due to transformative applications in climate modeling [1-5], molecular dynamics [6-8], cosmology [9-11], and beyond. These surrogate models or emulators are often orders of magnitude faster than classical numerical solvers and can be trained to directly emulate real-worke data. Although existing developments in neural surrogate models have made notable strides in short-term prediction, often benchmarked using one-step mean square еитог (MSE), their ability to deliver stable long-term predictions remains critical challenge [12-14]. For scientific practitioners, achieving a long-term stable forecast valuable and expands the potential for research use, especially where classical numerical solvers are too costly to run. For example, climate scientists can use accurate long-term emulators to predict the probability of extreme events [15]. However, the error accumulation in unrolling the predictions for complex, nonlinear, and multiscale dynamical systems often causes existing surrogate models to degrade catastrophically during long-term forecasting, resulting in unphysical collapsed or exploding solutions [14, 16).
Recent efforts address long-term stability explore two primary avenues. From the temporal evolution perspective, frameworks centered on multi-step ahead prediction aim to shift focus from single-step accuracy extended horizons by training models on multi-step targets [17-20]. Yet,
"Correspondence roxie_jiang deciiin@@udan..co www.aa.uu...nnnn...mmmmmmm whhnnn@@chiccgoocoommm edu, and ciilett@uchicago.co 'https://github httrp:///gitabb..omo/o///e)/ www/roxie62 wtw.eeeeeeeeeeeeriiiiaaaaaaa Implicit CeriiccctNNeurrl. Emultorrs
39th Conference on Neural Information Processing Systems (NcurIPS 2025).

---

Autoregressive Rollout and Residuals
Extreme Long-term Rollout
Baseline
Ours
25
Difference
50
Difference
75
Difference
5000
10K VL 50K
100K
Figure I: In a chaotic, turbulent system, our emulator achieves accurate short-to-mid-term rollout predictions (left). Even over extremely long sequences (up to 109 emulation steps, right), it captures the physical jet patterns, while baseline autoregressive methods quickly drift and break down.
such methods can underperform even one-step ahead autoregressive baselines for short-term complex dynamics [21]. From the spatial domain perspective, architectures that use frequency decompositions long-term statistics training targets advance the performance of surrogate modeling to better handle the multiscale structure of the data [22-26]. However, they often face trade-off between short-term accuracy and long-term robustness, or require additional physical knowledge.
In this work, we draw inspiration from numerical analysis and reframe neural surrogate modeling through the lens of time-stepping methods for solving differential equations. Viewing the conventional autoregressive models from this perspective, we can interpret their failure patterns as analogous to the instability of explicit schemes (e.g., forward Euler), especially for larger timesteps and in stiff systems. To improve stability (Fig. 1), we propose to enhance the autoregressive system based on an analogy with the implicit time-stepping methods, which are known for their unconditional stability in challenging integration scenarios [27].
The comerstone of our approach is the use of future frame information to guide next-step predictions, fundamentally departing from existing autoregressive frameworks constrained to the use of past temporal states. While inspired by numerical implicit schemes that iterative solvers to refine pre- dictions, our framework avoids the corresponding computational overhead by effectively performing the iterative refinement steps in parallel. This is achieved by simultaneously forecasting a hierarchy of future states during autoregressive rollout (Fig. 2), resulting in nearly the same number of forward passes compared the one-step ahead baseline model. Furthermore, unlike prior multi-step ahead approaches [17-19], our framework actively compresses future states into a hierarchical representa- tion, balancing accuracy and computational tractability while enabling feedback from future time frames. On 2D turbulence with multiscale and chaotic structure, our method demonstrates greatly enhanced long-term stability compared to the baseline models (Fig. right), while introducing minimal amount of additional computational overhead.
We summarize contributions follows:
Neural emulation with an implicit schema and a hierarchical structure. We propose a novel emulator that uses augmented state built from hierarchy of coarse-grained representations of future states. This approach captures multiscale interactions and - by autoregressively applying the emulator to the augmented state efficiently mimics the iterative refinement mechanism of implicit solvers, achieving better accuracy and stability.
Accurate simulation of chaotic systems. Our approach is well-suited for turbulent flows, which chaotic and have complex multiscale structures. On 2D Navier-Stokes system with Reynolds number of 10 our model achieves substantial improvements - reducing mean squared error by over 50% for predictions 25 50 steps, compared to the standard 1-step ahead baseline.
Stability for extremely long rollouts that are 10x the training length. Our emulator exhibits strong stability over extremely long sequences, For the 2D Navier-Stokes system at 10* Reynolds number, it rolls out up to 2 x 10 steps (one time the length of the training set) without empirical failure. Even at 2 x 10a steps (ten times the length of the training set), it maintains a low failure rate of just 7.0 7.0%, far outperforming existing baselines, which typically fail much earlier and for all initial conditions.

---

timestep
timestep
Y1
71 n+
n +
72 1
?
++2
n +
U-Net
(2)
fo
((2
8n42
*n43
2(1)
R442
-10)
Un Zn
un
Figure 2: Diagram of our hierarchical implicit neural emulator. Our model approximates the dynamics using an implicit formulation with hierarchical latent representations. For each trajectory Un. we augment its state with multiscale latent features zn obtained by downsampling un. The model predicts the future state ú and long-horizon latent representations 2 2 enabling it to capture long-term dynamics. Following the implicit schema, the model conditions on both the past trajectory u, and future latent variables =n- 2n+2 during training, while using predictions of future states computed in the previous step of an autoregressive rollout during inference, providing richer context to effectively mitigate enror accumulation for long-term predictions.
(2)
2 Related Work
Hierarchical networks. Many physical systems exhibit patterns across spatial and temporal scales. To effectively emulate such dynamical systems, the model architecture should reflect their inherent hierarchical structure. A U-Net [28], for instance, introduces skip connections between different spatial resolutions, enabling effective reasoning across both high-level semantics and fine-grained details. Many subsequent approaches similarly integrate hierarchical structure into a system's 's data representations or computational elements. Feature Pyramid Networks [29] fuse multiscale features to improve prediction accuracy. Multiscale Vision Transformers [30] construct hierarchical representations tailored for attention mechanisms. Inspired by multigrid solvers [31], multigrid CNNs [32, 33] maintain multiscale representations at each layer and introduce cross- scale convolution operations to enable efficient information routing via bidirectional connections between adjacent scales. Hierarchical VAEs (HVAEs) [34] extend latent variables across multiple scales, significantly improving generative model quality. However, HVAES often face challenges such as high training variance and instability. Nested Diffusion Models [35] address similar challenges by pretraining and freezing the hicrarchical latents, then training a collection conditional diffusion models to generate each level of the latent hierarchy conditioned on those above, concluding with image generation.
Neural PDE solvers. Neural PDE solvers offer dataen alternative to traditional numerical methods, with promising applications in real-world applications, such as weather forecasting (e.g.. FourCastNet [1], GenCast [36)). These methods fall into two main classes: (1) neural operators, which map between function spaces (e.g., FNO [22], DeepONet [37]), and (2) grid-based architectures, which operate on discretized grid or mesh (e.g.. U-Nets [28]). While neural operators promise mesh- agnostic inference, grid-based methods remain prevalent due to their simplicity and popularity in computer vision. Although our current implementation assumes grid-structured inputs, the framework can be naturally extended to mesh-based representations, as demonstrated in prior work [38-40].
major challenge for both long-term stability, where small errors amplify over time. Stabilization efforts include: (1) dltt--riven methods like pushforward training [18] and PDE-Refiner [41], which improve rollouts but suffer from instability and high computational cost; (2) physics-informed strategies such as multi-time stepping [42], which embed physical priors but limit data-driven flexibility; and (3) architectural advances like FNOs [22] and Wavelet Neural Operators [23], which improve spatial modeling but struggle with temporal coherence and optimization. Despite progress, current approaches fall short in handling multiscale chaotic systems, such as turbulent flows, where spatiotemporal complexity demands both efficient learning and robust long-horizon predictions.
3

---

Problem Formulation
Dynamical systems. Consider a nonlinear dynamical system evolving over time as:
Ou - f(, Ui V..... uaa, 0) uo:
Öt
where / encodes the unknown governing physics. Let u(x, 1) E IR denote the continuous trajectory of dynamics in function space Ut and u, R represent discretized snapshots (i.e., sample from a finite spatial grid with m points) at time t. To numerically approximate the dynamics, we adopt
discrete methods (c.g., finite differences, spectral methods), transforming the PDE into system of ODEs:
du: = f(ue, t), Ug E R
(2)
dt
where corresponds to the temporal dynamies discretized states. Given N + consecutive observations {uo, Ui..... uNJ. our goal is to learn neural emulator fe that approximates the underlying dynamics and predicts future states. The emulator parameterized by
Autoregressive models. Autoregressive methods learn a transition operator fo R - R to iteratively predict next state ú,+1 from Un
in+1 = fe(un);
(3)
and chain predictions via recursive rollouts:
ün+* = folfolfo... (u.))) for ls steps.
(4)
While efficient and accurate for short-term predictions, autoregressive methods often suffer from error accumulation over long horizons, where small discrepancies compound over time [18, 14]. Our approach adheres to this autoregressive framework, but introduces mechanisms mitigate these long-term errors, as described in Section 4.
Method
In this section, we show how to improve short-term forecast and preserve long-term coherence in an autoregressive model by building implicit multiscale framework.
4.1 Implicit Neural Emulator
Explicit implicit methods. For the semi-discretized system in Eqn. explicit time-stepping methods [43] compute future states directly through:
Un+1 Ri F(u.);
where the mapping F-) represents an explicit update rule depending only on the current state Up.- For instance, in the forward Euler method [44], it takes the form:
F(u.)) u. Atf(u.)
(6)
to compute an estimate of un++- This approach can be connected to the autoregressive learning seen in Eqn. 3, where predictions rely solely on past states. And similar autoregressive methods, is known to suffer from the instability using a larger time step Al or facing systems with stiffness [45], where adaptive and higher-order variants usually required ensure stability of the solver, which substantially increases the computational cost.
Implicit methods are designed solve these issues: they leverage both current states u. and future states +1 to solve the following implicit equation:
U
F(un; Un+1).
Here F(u,, uh++1) u, + At Atf in the backward Euler method. Solving u441 in Eqn. requires extra computations to find the root. A classic example Newton's method, which iteratively updates the estimate via:
F F(uhp: w
8
(uu u 2

---

where i indicates the iteration index. Compared to explicit methods, implicit approaches can accommodate much larger timesteps while maintaining stability, offering computational efficiency for long-time integration.
Implicit neural emulator with two-step prediction. Drawing inspiration from the robustness of implicit solvers, we propose neural network architecture that combines the stability of implicit methods with the simplicity of explicit methods. Our goal is address the compounding error typically seen in autoregressive models by introducing a structure that implicitly reasons about future states. Rather than solving the implicit equation in Eqn. directly - which would entail iterative root- solving at each step we adopt relaxed version. We introduce latent variable 2x+1 = T(un+1), which represents an abstract encoding of the future state un+1- This latent variable can be directly computed during training, but must be predicted at test time, leading to the following formulation:
un+11 2n = So (Un, 1
where the network is trained to simultaneously predict the next physical state ún+1 and the latent representation 2, that encodes information about future state.
This architecture mirrors the iterative refinement process used in implicit solvers. Our model given an initial approximation of state n + 1 represented by En+1 and then refines this approximation by applying fun,, zn+1) to produce the predicted physical state un+1- By also providing 2x+2 as auxiliary output, the model is able to perform this refinement process as part of the standard autoregressive rollout across timesteps, effectively treating (u,,, En+1) as an augmented state variable.
In practice, we choose the transformation T to be downsampling or coarse-graining operation. Predicting a coarse- coarse-grained ε rather than the full physical state ün+2 encourages the emulator to first capture large-scale spatial features and, as shown in our experiments, leads to better prediction accuracy. We will elaborate on the design of the transformation T in a later subsection.
4.2 Hierarchical Multi-step Prediction
The architecture above naturally extends to a hierarchical multi-step modeling framework in which predictions are conditioned multiple latent representations of anticipated future states. We denote these representations as ( ( (u,)). where indexes increasing levels of abstraction and we
(0) The model is then trained to predict across I, hierarchical levels as follows:
assume um
(1) Un 2222 -447 = So Un, 2nt1............................ ...........
(1)
(L-1) -.
(10)
As shown in Figure 2, our model effectively operates in an augmented state space formed by the concatenation of the current state with future latent states of increasing levels of abstraction. In practice, this simple strategy proves scalable and robust, and easily extensible to deeper hierarchies.
This generalized framework enhances the analogy to the multi-step iterative refinement used in implicit solvers while also enabling the model to capture broader temporal dynamics, which helps mitigate error accumulation by informing each prediction with richer view of future latent dynamics.
Leveraging future context in a hierarchical fashion. Using future context offers two complementary benefits. From encoding standpoint, it allows the model to process input structured sequence where immediate states like u, retain fine-grained detail, while distant latents such as z(() coarse-scale insights like aggregate dynamics long-term trends. This mirrors the philosophy of multi-step implicit methods like
#++ provide
[46], which integrate information from multiple past states to produce more accurate future estimates, though with slight distinction that conditioning on the future frames, while being able to extend to include the history states as well.
From a decoding standpoint, learning to predict distant states becomes progressively harder due to increased uncertainty and error amplification. By supervising the model across multiple levels of abstraction, we encourage balanced learning process where both local precision and global structure are prioritized. Viewed through the lens of implicit solvers, this can be interpreted as executing steps of iterative refinement, distributed across temporal frames and abstraction levels.
4.3 Training Details
Choosing abstract transformations. A key component of training our model is the choice of the transformations T()), which map physical states un to latent representations zn. There
5

---

Confi. Name
Model
Config. Name
Model
Ours
a as 2
(t Snos)
I.
Mniii #0421 4n43 Solu, 4 x
Spatial Hierarchy
Ma+t1 2 [44) faluen 1/ 1
Baseline
2-Sup Ahead
2 (Mbr)
History Hierarchy
fo(Mn Fy Fu.22)
2-Step History [19)
dutt1 to (Me (by))
Step History [19]
Mp+1 fo (M 2: Ma- st.))
Table 1: Model and baseline configurations. We conduct experiments with our model designs
(L = 2, L = 3) and several baseline configurations for comprehensive analysis.
2-step
History-tieraachy
Spatialkieerarhy
2-step History
3istep History
1.0
0.0
Stabiitt
0.4
0.2
0.0
20
40
81
100
100
Rollouts time steps
Rolkut time stoss
(a) Performance
estimation measured by mean squared error (MSE)
(b) Stability for long-term rollout, up to 2 x steps
100-step autoregressive rollout. Solid lines indicate the average performance, and shaded regions represent one standard deviation. Right: We evaluate the stability rate over 10 times the training sequence length across 100 trials with various initial conditions. Stability is defined by maintaining physical statistics (conserved kinetic energy). Our method with L 3 achieves 93% stability rate, significantly outperforming all baselines, with the next best, 2-Step History, at just 38%, indicating superior short-term accuracy and long-term robustness.
Figure 3: Short-term accuracy long-term robustness. Left: MSE trend
several viable options for T such as a learnable encoder or a predefined informatiionrreduction function (e.g., spatial pooling). While learned encoder offers expressiveness, it introduce instability during training due to the constantly shifting latent representations. Instead, we focus on the multiscale structure of spatiotemporal systems and choose (()) to be a set of fixed spatial downsampling transformations:
reduction
(1) (1) (un) = DownSample(un: r"),
(11)
Here p denotes the predefined downsampling rate for the 1-th level with of < ,h, providing sequence increasingly coarse-grained states.
Training objective. Finally, our training loss is designed to supervise both the predicted physical state and the associated abstract latents:
{(O) d(ún+1 un+1) + Σ d(>)) 2
(12)
where d(.) denotes a distance metric such as [1 or (2 loss for simplicity, though this could be replaced with domain domain-specific alternatives. This objective ensures that the model's predictions remain aligned with both immediate and long-term dynamics, across all abstraction levels.
Experiments
In this section, we benchmark our implicit modeling framework for turbulence prediction. This canonical flow has been extensively used for testing novel ML-based schemes [47-49, 41. 50].

---

Baseline: L-1
Ours: L-3
Histor- Hierarchy
2-step History --- Reference truth
L=2
2-stop Ahead
Spatial-Hicrarchy
History
Error Energy: Spectruum
Energy Spectrum
26--
14-*
yoer woodt and that
Wavenumber
e Zonal Mean
(a) Zonal Mean
(b) Energy Spectrum
(c) Error of Energy Spectrum
Figure 4: Further investigation of long-term stability of 2 x 10 steps (10x the training length) rollout. (a): Time-averaged zonal mean of vorticity comparing ground truth (dashed black lines) with emulator runs. Only our systems with L = 2 and L = 3 accurately follow the ground truth trend, with the reference truth lies well within error bars, demonstrating that our emulator accurately captures the system's mean statistics. In contrast, other methods significantly deviate, indicating substantial drift and failure to preserve structure. (b,c): Spectrum of long-term rollout for normalized data.
Navier-Stokes. We focus on the dimensionless vorticity-ttreamfuuntion (w v/) formulation of the incompressible Navier-Stokes equations in a 2D x y domain [47, 48]:
Dw Ot
+ +W(,, = Re
V x + f + Buy:
(13)
where 2 = - v (U2, v2) is velocity, and w = V x v. N(w, +) is the Jacobian and represents non-linear advection. The flow is defined by a Reynolds number Re 10 time-constant forcing If at a given wavenumber, and Rayleigh drag = 0.1. The Coriolis parameter, 3 = 20, induces zonal jets characteristic of geophysical turbulence, mimicking the influence of Earth's rotation on atmospheric and oceanic flows [51]. The domain is doubly periodic with length L = 2n. Training
data is generated using "py2d" [52] on a 512 x 512 grid, with 18,000 w snapshots. For analysis in this paper, the data is downsampled to 256 x 256 grid, except in Section 5.3, which discusses ablation studies on different datasets and resolutions.
We adopt the UNet design following Dhariwal and Nichol [53] as our base model, following the standard setup with using a structured dataset. Specifically, our model has an encoder composed of 5 downsampling blocks that reduce the input to an 8 x x 256 feature map at the bottleneck. The decoder mirrors the encoder structure reconstruct the output. To improve the model's 's ability to capture high-frequency details, we incorporate Fourier layers [22] into each block to better handle high-frequency information (see details in Appendix A.3).
Architecture.
Given that our design must accommodate images at multiple resolutions, we align with the UNet's hierarchical information flow: lower-resolution inputs are introduced at the earlier encoder stages, while decoder blocks operating at matching resolutions are used to produce the final output. This structure ensures that the most abstract features are captured at the encoder's bottleneck, promoting a coherent and structured flow of information throughout the network.
Model configurations and baseline setup. Our experiments primarily foeus designs with two- level and three-level hierarchies, indicated by L = 2 and L = 3. For the input with 256 x 256 resolution, choose r' 8 and , 32 as default downsampling rate for latent variables (1)) (2)
242 respectively. To rigorously assess our model's performance, we also propose several baseline methods that diverge from our formulation, including different setups of input/output hierarchy configurations, as detailed in Table 1.
Hierarchical autoregressive rollout. To enable autoregressive rollout given only the initial condition u, we introduce specially-designed corner-case inputs, corresponding to --1 higher-level latent slots, during both training and inference. During the prediction, with small probability p (see details in Appendix A.3), we sample training instances which the model receives partially missing hierarchy of latent inputs and is tasked to reconstruct the missing parts. Concretely, for our I. = 2
7

---

10C
BUNEELL Can
O
(a) PCA autocorrelation along (b) PCA autocorrelation errors
(c) Averaged energy spectrum (d) Energy spectrum error over
rollout steps
along rollout steps
over wavenumber
wavenumber
Figure 5: Ablation studies on the design of the hierarchy. We evaluate our system across different hierarchy levels L and visualize the results using: (a) PCA autocorrelation and (c) energy spectrum averaging over 200 rollout steps, with corresponding deviations from the ground truth shown in (b) and (d), respectively. Our hierarchical design with L outperforms the baseline L 1, achieving better alignment with temporal structure and high-frequency components.
25
50
75
100
04
0.865e-02
(2.890c-02)
2.103c-01
4.470=-01 (1.763c-01)
7.738e-01
(2.520e-01)
1.747e-03 (3.942e006 04)
1.140e001 (6.328e-02)
3.774e-01
(1.797e-01)
7.082e-01
(3.088c401)
1.023e (3.355ee-0
7.735e-04
(1.954e-00 04)
5.651e-02 (3.199g-022 02)
1.998-001
(1 (1104e-01)
4,225;-01 01 (1.979c-01)
7.280=-01
(2.522e-01)
5.248e-04
(1.148e-04)
4.024c-02
(1.920e-02)
16066-01
(6.496e-42)
3.731e-01
(1.509c-01)
6.547e-01
(2.440c-01)
16
5.155c-04 (4.138e-04)
4.898c-02 (2.01le-02)
2.039c-01 (8.009e 02)
4.416e-01 (1.789e-01)
7.399e-01 (2.658e-01)
Table 2: Ablation study on downsampling ratio (r*) for our L = 2 model on 256 x 256 resolution data with jet. We evaluate model performance by computing roll-out mean squared error (MSE) for downsampling ratios r' = 1, 2, 4, 8, 16, and present the mean (standard deviation) over 100 trials with varied initial conditions. Notably, r¹ exhibits the minimal error, indicating an ideal balance between preserving contextual information and mitigating estimation challenges.
levels design, we present the model with [un, 0j. where 0 denotes the zero-filled tensor matching the shape of the higher hierarchy latent 2(1), and train it to predict a This ensures that the model learns to infer higher-level latent states from just the initial condition. At the evaluation time, each method is provided only with the initial condition u, In the hicrarchy L = 2, we first input [un, 0] to obtain 2 Then, using the full ppatial-temporty hierarchy states [u,, 2ntl we continue with autoregressive rollout.
5.1 Short-term Accuracy and Long-term Stability
We assess our model's performance under two scenarios: (1) Short-term predictions, up to 100 steps using mean squared error (MSH); and (2) Extremely long-term predictions, generating up to 2 x 10 steps to evaluate the alignment of predicted and true physical statistics.
Evaluation short-term estimation. Figure 3(a) presents MSE results comparing method with various baselines. All methods show similar 1-step MSE, yet the 2-Step Ahead Baseline performs poorly (see Table 6 in Appendix A.I). Our models with L 3 and L 2 consistently provide accurate estimates, surpassing all baselines despite challenges in recursive prediction error accumulation. Baselines using multiple input or output frames (e.g., 3-Step History or 2-Step Ahead) perform better for mid-term rollout, highlighting the importance of processing temporal sequences in autoregressive models. In comparison, basic hierarchical baselines (e.g., Spatial-Hierarchy) provide minimal gains, suggesting that temporal hierarchy is crucial.
or
Evaluation on long-term estimation. Over very long rollouts, mean squared error (MSE) becomes an unreliable metric for evaluating autoregressive systems. Instead, we verify if the predictions align with the true physical statisties, such as the system's 's energy given by E = (v + v.)). For statistically stationary turbulent system, like the one studied here, this energy stays around a constant mean due to the balance between energy input by the deterministic forcing and its dissipation via viscosity and the linear drag in equation 13 [54].
Stability rate. We use a simple binary energy-based approach to check if prediction adheres to correct physical patterns. By calculating the mean and standard deviation of energy from a 200K-step true trajectory, we ensure all true trajectory energy values fall within 4 standard deviations of the mean.
8

---

1425
1426
12g 225
4225
M25
140
4550
M450
M450
1goo
14.0
255 256, No Jet 256 x56xx 255, Jat
GT
Baseline: I,
Ours: L
3
GT
Baseline: I,
Ours: =
6450
M4
Mko
úoo
úooo
450
512 512x Jet
GT
Baseline:
Ours:
Figure 6: Visualization of rollout estimation across multiple datasets. We apply our approach to three different flows: (1) Re 10 256 256 resolution featuring zonal jets, (2) Re
10 256 x 256 resolution without zomal and (3) Re 10ah 512 x 512 resolution with zonal jets. Our method (L = gives more accurate predictions with lower associated residuals in all three scenarios.
A prediction is deemed stable if its energy remains within 5 standard deviations of the true mean (see Appendix A.2). Appendix B.2 includes examples with varying energy deviations, illustrating that predictions exceeding this threshold exhibit incorrect physical patterns.
Figure 3(b) presents the stability rates for all methods, based predictions over 200K steps from 100 initial conditions, evaluated every 100 emulation steps. Our approach with L = 3 consistently sustains a high stability rate during rollout, achieving 93.0% stability through the 10 steps. This greatly surpasses the closest competitor, the 2-Step History method, which only reaches 38.0%, with 2-Step Ahead and our their limitations in capturing long-term dynamics. These findings highlight the need for our methods that can handle longer temporal dependencies ensure robust long-teni performance.
following. Most other methods fail the end of rollouts, indicating
Zonal mean and energy spectrum. We extend our evaluation using (1) the zonal mean and (2) the energy spectrum (see definitions in Appendix A.2). The zonal mean reveals the flow's large-scale organization, namely, jet formation and maintenance. In geophysical turbulence, jets are crucial structures resulting from nonlinear advection, planetary rotation, and dissipation. Capturing the zonal mean is vital for verifying that the emulator preserves these structures over time. As seen in Fig. 4(a), and L 3) produce zonal mean profiles within reproducing the jets, unlike baseline (L 1) and other methods. The energy spectrum illustrates kinetic energy distribution over wavenumbers. In Fig. 4(b), L and L exhibit minimal bias, aligning closely with the reference spectrum, even at small wavenumbers. These findings confirm our model's robustness, maintaining long-term stability while accurately reflecting the system's 's statistical properties.
methods (L
bands, successfully
5.2 Ablation Studies
Ablation model hierarchy. A key aspect of our model is its ability process and estimate hicrarchical representations the data. To further examine this design, we analyze PCA autocorrela- tion to capture temporal structure and the energy spectrum evaluate frequency-domain accuracy, Computation details are provided in the Appendix A.3. The results are shown in Fig. 5 and Fig. Our model with L = preserves temporal patterns more effectively and provides more accurate estimates
components compared to = and L =
Ablation on downsampling factor. Our system eases the future prediction tasks by initially assessing a downsampled version the future frame. The downsampling ratio plays an important role in balancing estimation challenges and available information for next frame prediction. We

---

Rollout Step
25
75
100
7 (8 do loa)
5.600-04
1.8e-02 (4 (7)-(0)
8.04e- 02
-50-00)
5.12=-01 (1. (1.15e-01)
60e-00 (2.08e-01)
9.38x-01
(2.58-011)
7-
(Sobolee lood
7.01c-04 D4 (4 07)-00)
1644-01 (8 (1)
8.76--01 01 45c-013
1.40+000 (2. (2.00-01)
1666400
(2.56e-01 01)
1.90400 (2.51e001)
5.50e-04
20g-04)
6,5 6.5%
05 (0)
5.3%
02 (1.41e
1.0g-00 (6.16-002)
224--01 (1.38
5.92- (1.34g-01)
2- 3 (Soboler loes)
6.89e-04 04 (303e- 04)
1,45e
(6386- (3)
7.38
02 (3800
264c 01 (1.42-001)
074--01 (2.78e
883 (3.31--01)
Table 3: Comparison of rollout MSEs with varied training objectives. We train the models using (1) standard & + Eg loss functions; (2) Sobolev loss from et al. [25].
Method
Seconds per iter.
1-step MSE
25 25-step MSE
50-step MSE
Zenal mean ceror
Raseline:
0,2067
5.60e 04 (1.15e 04)
8.04e (3.45e 02) 3.12e-01 (1.15e-01 01)
4.20 (2.36)
Pushforward [18)
0,2834
1.08--03
(2.26e-01)
9.162-02 (4.62c-02)
3.19c-01 (1 (1334-01)
1.07 (0.44)
Ours: =
0.2204
5.50e-04 (1.20e-04)
3.37e-22(1.44e-222
1.40e-01 (6.16e-02)
0.63 (0.58)
Table 4: Comparison of training efficiency and multi-step forecasting errors. Our T. = 3 method
yields better short-term accuracy and preserves better structure coherence, compared to the T. baseline and the pushforward unrolling [18].
512 512 resolution, wjet
256 256 resolution wio
Rollout Step
Baseling: & =1 1 2.16c-01 (8. 77c-12)
6.82x-01 (2.55c001)
1.13.400 (3. (3.54c-01)
2.59c-01 (1 (1.09x00)
8.44-001 (2.85c001) 5.12c-01 (1.76e-01)
1.424400 (3 92c-01)
Ours: J. = 3
8.75c-02 (3.04c-02)
1.81e111 (9.02-012)
3896-011 (1.83e-01)
1.28--01 (5.136-112)
9, 9.98e-01 (2.55e-01)
Table 5: Mean squared error on different datasets for short-term rollout. We report mean (standard deviation) across 100 trials of short-term rollout with different initial conditions.
ablate our two-level system across various downsampling rates of r [2, 8, 32]. Results presented in Table 2 demonstrate that r¹ 8 provides the best performance for nearly all prediction steps.
A Alation on training objectives. In addition the standard /1 and 12 losses, we further provide results for experiments with the Sobolev loss, proposed in Li et al. [25]. As shown in Table 3, across different training objectives, our I. = 3 method consistently yields significant improvements compared to the L = baseline.
Comparison on training efficiency and rollout method. Unlike multi-step unrolling methods [18, 20), our approach achieves greater computational efficiency while enhancing both short-term accuracy and long-term coherence. As illustrated in Table 4, it incurs only about 10% higher cost per iteration relative to the L = baseline, yet substantially reduces errors across both short- and long- long-horizon metrics. In contrast, the pushforward method [18] mitigates distributional shift at the expense of higher computational cost and degraded short-term accuracy.
5.3 Performance on Diverse Datasets
We further evaluate the robustness of our approach across multiple datasets and spatial resolutions.
In addition the jet-forming regime with 20 that mimics the influence of Earth's 's rotation, we
consider an eddy configuration (3 = 0). As shown in Fig. 6, our approach (L = 3) consistently produces realistic and accurate flow evolutions.
contrast, the baseline model (I. = 1) shows large residuals after a short rollout. This performance gap is further quantified in Table 5, where the MSE of I. grows rapidly over time, while I, maintains low and stable errors.
Conclusion
In this work, we present a novel framework for autoregressive modeling of dynamical systems, inspired by implicit time-stepping numerical solvers. Diverging from prior approaches, our method uniquely integrates future-state insights to refine next-step predictions. To enable autoregressive operation, we introduce a temporal hierarchy via spatial compression, coupled with an encoder- decoder architecture. This design organizes multiscale encoder to distill critical features and the decoder to adaptively address prediction challenges across scales. Evaluated on turbulent flow systems, our framework achieves substantial error reduction while maintaining computational efficiency. While our current focus centers on deterministic modeling, the architecture's 's flexibility suggests promising extensions to generative paradigms (e.g., diffusion models) for probabilistic forecasting and uncertainty quantification.
tempoall information, empowering the
10

---

Broader Impacts. While better emulators for dynamical systems could be used in a wide range of applications, foresee no direct negative societal impacts.
Acknowledgements
RJ. PI., and RW gratefully acknowledge the support of AFOSR 18-10068--- the Eric and Wendy Schmidt Al in Science Fellowship program, and the Margot and Tom Pritzker Foundation. RJ, PL., RW, and MM gratefully acknowledge the support of the NSF-Simons Al-Institute for the Sky (SkAI) via grants NSF AST-2421845 and Simons Foundation MPS-AI-00010513, KJ and PH were supported by NSF RISE-2425898 and Schmidt Sciences, LLC.
References
[1] Jaideep Pathak, Shashank Subramanian, Peter Harrington, Sanjeev Raja, Ashesh Chattopadhyay. Morteza Mardani, Thorsten Kurth, David Hall, Zongyi Li, Kamyar Arizzadeneshelii et al. Fourcastnet: A global data-driven high-resolution weather model using adaptive fourier neural operators. arXiv preprint arXiv:2202.1 1214, 2022.
[2] Ruoxi liang, Peter Y. Lu, and Rebecca Willett. Embed and emulate: Contrastive representations for simulation- bbased inference. arXiv preprint arXiv:2409. 0224099.84400 /8402, 2024.
[3] Remi Lam, Alvaro Sanchez-Gonzalez, Matthew Willson, Peter Wirnsberger, Meire Fortunato, Ferran Alet, Suman Ravuri, Timo Ewalds, Zach Eaton-Rosen, Weihua Hu, et al. Learning skillful medium-range global weather forecasting. Science, 382(6677):1416-1222) 2023.
[4] Dmitrii Kochkov, Janni Yuval, Ian Langmore, Peter Norgaard, Jamic Smith, Griffin Mooers, Milan Klöwer, James Lottes, Stephan Rasp, Peter Düben, al. Neural general circulation 2024.
models for weather and climate. Nature, 632(6002)
[5] Oliver Watt-Meyer, Gideon Dresdner, Jeremy McGibbon, Spencer K Clark, Brian Henn, James Duncan, Noah D Brenowitz, Karthik Kashinath, Michael S Pritchard, Boris Bonev, et al. Ace: A fast, skillful learned global atmospheric model for climate prediction. arXiv preprint arXiv:2310.02074, 2023.
[6] Albert Musaelian, Simon Batzner, Anders Johansson, Lixin Sun, Cameron J. Owen, Mordechai Kornbluth, and Boris Kozinsky. Learning local equivariant representations for large-scale atomistic dynamics. Nature Communications, 14(1):579, Feb 2023. ISSN 2041-1723. doi:
10.1038/41467:0333363339
[7] Raidel Martin Martin-Barrios, Edisel Navas-Conyede Xuyi Zhang, Yunwei Chen, and Jorge Gulín- González. An overview about neural networks potentials in molecular dynamics simulation. International Journal of Quantum Chemistry, (2<(11):,2738, 2024.
[8] Raimondas Galvelis and Yuji Sugita. Neural network and nearest neighbor algorithms for enhancing sampling of molecular dynamics. Journal of chemical theory and computation, 13 (6):2489-250) 2500, 2017.
[9] Drew Jamieson, Yin Li, Renan Alves de Oliveira, Francisco Villaescusa-Navaroo Shirley Ho, and David N Spergel. Field-level neural network emulator for cosmological n-body simulations. The Astrophysical Journal, 952(2):145, 2023.
[10] Drew Jamieson, Yin Li, Siyu He, Francisco Villaescusa-Navaree Shirley Ho, Renan Alves de Oliveira, and David N Spergel. Simple lessons from complex learning: what neural network model learns about cosmic structure formation. PNAS nexus, 2(4):pgac250, 2023.
[11] Siyu He, Yin Li, Yu Feng, Shirley Ho, Siamak Ravanbakhsh, Wei Chen, and Barnabás Póczos. Learning to predict the cosmological structure formation. Proceedings of the National Academy of Sciences, 116(28):13825-1382 2019.
[12] Jeongjin Park, Nicole Yang, and Nisha Chandramoorthy, When are dynamical systems learned from time series data statistically accurate? arXiv preprint arXiv:2411.06311. 2024.
11

---

[13] Massimo Bonavita. On some limitations of current machine learning weather prediction models. Geophysical Research Letters, 5((1)):2023GL107777, 2024.
[14] Ashesh Chattopadhyay, Y. Qiang Sun, and Pedram Hassanzadeh. Challenges of learning multi- scale dynamics with ai weather models: Implications for stability and one solution, 2024. URL https://arxixs
orgooog//ao ht//bb//2300 07029.
[15] Ching-Yar Yao Lai, Pedram Hassanzadeh, Aditi Sheshadri, Maike Sonnewald, Raffaele Ferrari, and Venkatramani Balaji. Machine learning for climate physics and simulations. Annual Review of Condensed Matter Physics, 16, 2024.
[16] Niraj Agarwal, Dmitri Kondrashov, Peter Dueben, Eugene Ryzhov, and Pavel Berloff. A comparison of data-driven approaches build low-dimensional ocean models. Journal of Advances in Modeling Earth Systems, 13(0)02202 2021.
[17] Jayesh K. Gupta and Johannes Brandstetter. Towards multii-pptiooeeporal-ss scale generalized pde modeling, 2022. URL https://aariiii org/otaaa hts///0b/22/ 15616.
[18] Johannes Brandstetter, Daniel Worrall, and Max Welling. Message passing neural pde solvers. arXiv preprint 27Xiv:202..03376, 2022.
[19] Makoto Takamoto, Timothy Praditia, Raphael Leiteritz, Dan MacKinlay, Francesco Alesiani, Dirk Pflüger, and Mathias Niepert. Pdebench: An extensive benchmark for scientific machine learning, 2024. URL
org/ barrg/absi httpa/ab///2200////
[20] Bjoern List, Li-Wei Wei Chen, Kartik Bali, and Nils Thuerey. Diffcrentiability in unrolled training of neural physics simulators on transient dynamics. Computer Methods Applied Mechanics and Engineering. 433: 17441, 2025.
[21] Salva Rühling Cachay, Bo Zhao, Hailey Joren, and Rose Yu. Dyffusion: dynamics-informed diffusion model for spatiotemporal forecasting. Advances in neural information processing systems, 01445599-42888 2023.
[22] Zongyi Li, Nikola Kovachki, Kamyar Azizzadenesheli, Burigede Liu, Kaushik Bhattacharya, Andrew Stuart, and Anima Anandkumar. Fourier neural operator for parametric partial differen- tial equations. arXiv preprint arXiv:2010.08899, 2020.
[23] Tapas Tripura and Souvik Chakraborty. Wavelet neural operator: neural operator for parametric partial differential equations. arXiv preprint arXiv:2205.02111, 2022.
[24] Peiyan Hu, Rui Wang, Xiang Zheng, Tao Zhang, Haodong Feng, Ruiqi Feng, Long Wei, Yue Wang, Zhi-Ming Ma, and Tailin Wu. Wavelet diffusion neural operator. arXiv preprint arXiv:2412.0488 2024.
[25] Zongyi Li, Miguel iiu-Schiaffini Nikola Borislavov Kovachki, Kamyar Azizzadeneshelii Burigede Liu, Kaushik Bhattacharya, Andrew Stuart, and Anima Anandkumar. Learning chaotic dynamics in dissipative systems. In Alice H. Oh, Alekh Agarwal, Danielle Belgrave, and Kyunghyun Cho, editors, Advances in Neural Information Processing Systems, 2022. URL net//orummidd
attps://openrewe
[26] Ruoxi Jiang, Peter Y. Lu, Elena Orlova, and Rebecca Willett. Training neural operators to preserve invariant measures of chaotic attractors. Advances in Neural Information Processing Systems, 36, 2024.
[27] Yinnian He. Unconditional convergence of the euler semi-implicit scheme for the three- dimensional incompressible mhd equations. IMA Journal of Numerical Analysis, 35(2):767-801, 2015.
[28] Olaf Ronneherger, Philipp Fischer, and Thomas Brox. U-Net: Convolutional networks for biomedical image segmentation. In MICCAI, 2015.
[29] Tsung-Yi Yi Lin, Piotr Dollár, Ross Girshick, Kaiming He, Bharath Hariharan, and Serge Belongie. Feature pyramid networks for object detection. In Proceedings of the IEEE conference on computer vision and pattern recognition, pages 2117-2125, 2017.
12

---

[30] Haoqi Fan, Bo Xiong, Karttikeya Mangalam, Yanghao Li, Zhicheng Yan, Jitendra Malik, and Christoph Feichtenhofer. Multiscale vision transformers. In IEEE/CVF International Conference Computer Vision, 2021.
[31] Michael Maire and Stella X Yu. Progressive multigrid eigensolvers for multiscale spectral segmentation. In Proceedings of the IEEE International Conference on Computer Vision, pages 2184-2199 2191, 2013.
[32] Tsung-Wei Ke, Michael Maire, and Stella X Yu. Multigrid neural architectures. In Proceedings of the IEEE Conference on Computer Vision and Pattern Recognition, pages 6665-6673, 2017.
[33] Tri Huynh, Michael Maire, and Matthew R. Walter. Multigrid neural memory. In Proceedings of the 37th International Conference on Machine Learning, 2020.
[34] Arash Vahdat and Jan Kautz. NVAE: deep hierarchical variational autoencoder. Advances in neural information processing systems, 33: 13611696699979 2020.
[35] Xiao Zhang, Ruoxi Jiang, Rebecca Willett, and Michael Maire. Nested diffusion models using hierarchical latent priors. arXiv preprint arXiv:2412.05984, 2024.
[36] Ilan Price, Alvaro Sanchez-Gonzalez, Ferran Alet, Tom R Andersson, Andrew El-Kadi, Dominic Masters, Timo Ewalds, Jacklynn Stott, Shakir Mohamed, Peter Battaglia, al. Gencast: Diffusion-bart based ensemble forecasting for medium-range weather. arXiv preprint arXiv:2312.15796. 2023.
[37] Lu Lu, Pengzhan Jin, and George Em Karniadakis. Deeponet: Learning nonlinear operators for identifying differential equations based on the universal approximation theorem of operators. arXiv preprint arXiv: 1910.0319, 2019.
[38] John W Ruge and Klaus Stüben. Algebraic multigrid. In Multigrid methods, pages 73-130. SIAM, 1987.
[39] Zongyi Li, Nikola Kovachki, Kamyar Arizzadenesheli, Burigede Liu, Andrew Stuart, Kaushik Bhattacharya, and Anima Anandkumar. Multipole graph neural operator for parametric partial differential equations. Advances in Neural Information Processing Systems, 33:6755-6766, 2020.
[40] Jie Chen, Yousef Saad, and Zechen Zhang. Graph coarsening: from scientific computing to machine learning. SeMA Journal, :9 2022.
[41] Phillip Lippe, Bas Veeling, Paris Perdikaris, Richard Turner, and Johannes Brandstetter. Pde- refiner: Achieving accurate long rollouts with neural pde solvers. Advances in Neural Informa- tion Processing Systems, 36:67398-67433 2023.
[42] Qi Wang, Yuan Mi, Haoyun Wang, Yi Zhang, Ruizhi Chengze, Hongsheng Liu, Ji-Rong Wen, and Hao Sun. Multipdenet: Pde-embedded learning with multi-time-stepping for accelerated flow simulation. arXiv preprint arXiv:2501. 2025.
[43] Uni M Ascher, Steven J Ruuth, and Brian TR Wetton. Implicit-explicit methods for time- dependent partial differential equations. SIAM Journal on Numerical Analysis, 32(3):797-823, 1995.
[44] Sigal Gottlieb, Chi-Wang Shu, and Eitan Tadmor. Strong stability high-order time discretization methods. SIAM review, 43(1):89-111 12, 2001.
[45] Jason Frank, Willem Hundsdorfer, and Jan Verwer. On the stability of implicit-explicit linear multistep methods. Applied Numerical Mathematics, 25(2-3): 205, 1997.
[46] Eugene Isaacson and Herbert Bishop Keller. Analysis of numerical methods. Courier Corpora- tion, 2012.
[47] Karan Jakhar, Yifei Guan, Rambod Mojgani, Ashesh Chattopadhyay, and Pedram Hassanzadeh. Learning closed-form equations for subgrid-scale closures from high-fidelity data: Promises and challenges. Journal of Advances in Modeling Earth Systems, 16(7):20223M5003874, 2024.
13

---

[48] Yifei Guan, Ashesh Chattopadhyay, Adam Subel, and Pedram Hassanzadeh. Stable posterioni les of 2d turbulence using convolutional neural networks: Backscattering analysis and gener- alization to higher re via transfer learning. Journal of Computational Physics, 458: 1 11090, 2022.
[49] Chris Pedersen, Laure Zanna, and Joan Bruna. Thermalizer: Stable autoregressive neural emulation of spatiotemporal chaos. arXiv preprint arXiv:2503.18731 2025.
[50] Romit Maulik, Omer San, Adil Rasheed, and Prakash Vedula. Subgrid modelling for two- dimensional turbulence using neural networks. Journal of Fluid Mechanics, 122:122,,44, 2019.
[51] Geoffrey K Vallis, Atmospheric and oceanic fluid dynamics. Cambridge University Press, 2017.
[52] Karan Jakhar, Rambod Mojgani, Moein Darman, Yifei Guan, and Pedram Hassanzadeh. py2d: High-performance 2D) Navier-Stokes solver in python (version 0.1), 2024. URL https:// gi github. htmm//nnfllidd/pppp....
[53] Prafulla Dhariwal and Alexander Nichol. Diffusion models beat gans on image synthesis. In NeurlPS, 2021.
[54] Stephen B Pope. Turbulent flows. Measurement Science and Technology, 12(11):2020-202, 2001.
[55] Zongyi Li, Nikola Borislavov Kovachki, Kamyar Azizzadenesheli, Burigede liu, Kaushik Bhattacharya, Andrew Stuart, and Anima Anandkumar. Fourier neural operator for parametric partial differential equations. In International Conference on Learning Representations, 2021. URL https: stttps://pperrvieess wtwp:eewwwnnefoouuneeeeee ottp:///wo./ooooooooooooooooaoaa
14

---

NeurIPS Paper Checklist
The checklist is designed to encourage best practices for responsible machine learning research, addressing issues of reproducibility. transparency, research ethics, and societal impact. Do not remove the checklist: The papers not including the checklist will be desk rejected. The checklist should follow the references and follow the (optional) supplemental material. The checklist does NOT count towards the page limit.
Please read the checklist guidelines carefully for information on how to answer these questions. For each question in the checklist:
You should answer [Yes], [No], or [NA].
[NA] means either that the question is Not Applicable for that particular paper or the relevant information is Not Available.
Please provide short (1-2 sentence) justificattio right after your answer (even for NA).
The checklist answers
integral part of your paper submission. They visible the reviewers, area chairs, senior area chairs, and ethics reviewers. You will be asked to also include (after eventual revisions) with the final version of your paper, and its final version will be published with the paper.
The reviewers of your paper will be asked to use the checklist one of the factors in their evaluation. While "[Yes] is generally preferable to "[No]" ", is perfectly acceptable to answer "[No] provided a proper justification is given (e.g., "error bars not reported because it would be too computationally expensive" " or "we were unable to find the license for the dataset we used"). ) In general, answering "[No] " or "[NA] is not grounds for rejection. While the questions are phrased in binary way, acknowledge that the true answer often more nuanced, so please just use your best judgment and write justification elaborate. All supporting evidence can appear either in the main paper the supplemental material, provided in appendix. If you answer [Yes] to a question, in the justification please point to the section(s) where related material for the question can be found.
IMPORTANT, please:
Delete this instruction block, but keep the section heading "NeurIPS Paper Checklist",
Keep the checklist subsection headings, questions/answers and guidelines below.
Do not modify the questions and only use the provided macros for your answers.
1. Claims
Do the main claims made in the abstract and introduction accurately reflect the paper's contributions and scope?
Question:
Answer: [Yes]
Justification:
In the abstract, we described our goal building a powerful neural emulator.
Guidelines:
The answer NA means that the abstract and introduction do not include the claims made the paper.
The abstract and/or introduction should clearly state the claims made, including the contributions made the paper and important assumptions and limitations. A No or
NA answer to this question will not be perceived well by the reviewers.
The claims made should match theoretical and experimental results, and reflect how much the results can be expected to generalize to other settings.
It is fine to include aspirational goals as motivation as long as it is clear that these goals are not attained by the paper.
2. Limitations
Question: Does the paper discuss the limitations of the work performed by the authors?
Answer: [Yes]
Justification:
We discuss the limitation in the conclusion.
15

---

Guidelines:
The answer NA means that the paper has no limitation while the answer No means that the paper has limitations, but those are not discussed in the paper.
The authors are encouraged to create a separate "Limitations" section in their paper.
The paper should point out any strong assumptions and how robust the results are to violations of these assumptions (e.g., independence assumptions, noiseless settings, model well-specification, asymptotic approximations only holding locally). The authors should reflect on how these assumptions might be violated in practice and what the implications would be.
The authors should reflect on the scope of the claims made, c.g., if the approach was only tested few datasets or with few runs. In general, empirical results often depend on implicit assumptions, which should be articulated.
The authors should reflect on the factors that influence the performance of the approach. For example, facial recognition algorithm may perform poorly when image resolution is low or images are taken in low lighting. Or speech-to-text system might not be used reliably to provide closed captions for online lectures because it fails handle technical jargon.
The authors should discuss the computational efficiency the proposed algorithms and how they scale with dataset size.
If applicable, the authors should discuss possible limitations of their approach to address problems of privacy and fairness.
While the authors might fear that complete honesty about limitations might be used by reviewers as grounds for rejection, worse outcome might be that reviewers discover limitations that aren't ' acknowledged in the paper. The authors should use their best judgment and recognize that individual actions in favor of transparency play an impor- tant role in developing norms that preserve the integrity of the community. Reviewers will be specifically instructed to not penalize honesty concerning limitations.
3. Theory assumptions and proofs
Question: For each theoretical result, does the paper provide the full set of assumptions and complete (and correct) proof?
Answer: [NA]
Justification:
We don't provide theoretical results.
Guidelines:
The answer NA means that the paper does not include theoretical results.
All the theorems, formulas, and proofs in the paper should be numbered and cross- referenced.
All assumptions should be clearly stated or referenced in the statement of any theorems.
The proofs can either appear in the main paper the supplemental material, but if they appear in the supplemental material, the authors are encouraged to provide a short proof sketch provide intuition.
Inversely, any informal proof provided in the core of the paper should complemented by formal proofs provided in appendix supplemental material.
Theorems and Lemmas that the proof relies upon should properly referenced.
4. Experimental result reproducibility
Question: Does the paper fully disclose all the information needed to reproduce the main ex- perimental results of the paper to the extent that it affects the main claims and/or conclusions the paper (regardless whether the code and data are provided or not)?
Answer: [Yes]
Justification:
We describe our detailed setup in our experiment section.
Guidelines:
The answer NA means that the paper does not include experiments.
16

---

If the paper includes experiments, a No answer to this question will not be perceived well by the reviewers: Making the paper reproducible is important, regardless of whether the code and data are provided or not.
If the contribution is dataset and/or model, the authors should describe the steps taken to make their results reproducible or verifiable.
Depending on the contribution, reproducibility can be accomplished in various ways. For example, if the contribution is a novel architecture, describing the architecture fully might suffice, or if the contribution is specific model and empirical evaluation, it may be necessary to either make possible for others to replicate the model with the same dataset, or provide access to the model. In general. releasing code and data is often one good way to accomplish this, but reproducibility can also be provided via detailed instructions for how to replicate the results, access to a hosted model (e.g., in the case of a large language model), releasing of a model checkpoint, or other means that are appropriate to the research performed.
While NeurIPS does not require releasing code, the conference does require all submis- sions to provide some reasonable avenue for reproducibility, which may depend on the nature of the contribution. For example
(a) the contribution is primarily a new algorithm, the paper should make clear how reproduce that algorithm. (b) If the contribution is primarily new model architecture, the paper should describe the architecture clearly and fully.
(c) If the contribution is new model (e.g., a large language model), then there should either be way to access this model for reproducing the results or way to reproduce the model (e.g.. with an open-source dataset instructions for how to construct the dataset).
(d) We recognize that reproducibility may be tricky in some cases, in which case authors welcome describe the particular way they provide for reproducibilityy In the case of closed-source models, it may be that access to the model is limited in some way (e.g., registered users), but it should be possible for other researchers to have some path to reproducing or verifying the results.
5. Open access to data and code
Question: Does the paper provide open access to the data and code, with sufficient instruc- tions to faithfully reproduce the main experimental results, as described in supplemental material?
Answer: [No] for the code.
Justification: We have disclosed details of experimental setups for reproducing the results. Our code will be released upon acceptance.
Guidelines:
The answer NA means that paper does not include experiments requiring code.
Please see the NeurlPS code and data submission guidelines (https://nips cc/ helgiddee/ooeeuubiisseeeee
for more details.
While we encourage the release of code and data, we understand that this might not be possible, "No" is an acceptable answer. Papers cannot be rejected simply for not including code, unless this is central to the contribution (c.g., for a new open-source benchmark).
The instructions should contain the exact command and environment needed to run to reproduce the results. See the NeurlPS code and data submission guidelines (https: httpccccppplccgggiiceeee
/CodeSubmii ssionPolicy) for more details.
The authors should provide instructions on data access and preparation, including how to access the raw data, preprocessed data, intermediate data, and generated data, etc.
The authors should provide scripts to reproduce all experimental results for the new proposed method and baselines. If only subset of experiments are reproducible, they should state which ones are omitted from the script and why.
At submission time, to preserve anonymity, the authors should release anonymized versions (if applicable).
17

---

Providing as much information as possible in supplemental material (appended to the paper) recommended, but including URLs to data and code permitted.
6. Experimental setting/details
Does the paper specify all the training and test details (e.g., data splits, hyper- parameters, how they were chosen, type of optimizer, etc.) necessary to understand the results?
Question:
Answer: [Yes]
Justification:
We provide details in the experiment section.
Guidelines:
The answer NA means that the paper does not include experiments.
The experimental setting should be presented the core the paper to level of detail that is necessary to appreciate the results and make sense of them.
The full details can be provided either with the code, in appendix, supplemental material.
7. Experiment statistical significance
Question: Does the paper report error bars suitably and correctly defined or other appropriate information about the statistical significance of the experiments?
Answer: [Yes]
Justification:
We report std in Figure 3, 4 Table 2, 3.
Guidelines:
The answer NA means that the paper does not include experiments.
The authors should answer "Yes" if the results accompanied by error bars, confi- dence intervals, or statistical significance tests, at least for the experiments that support the main claims of the paper.
The factors of variability that the error bars are capturing should be clearly stated (for example, train/test split, initialization, random drawing of some parameter, or overall run with given experimental conditions).
The method for calculating the error bars should be explained (closed form formula, call to a library function, bootstrap, etc.)
The assumptions made should be given (e.g.. Normally distributed errors).
It should be clear whether the error bar is the standard deviation or the standard error of the mean.
OK to report 1-sigma error bars, but one should state The authors should preferably report a 2-sigma error bar than state that they have 96% CI, if the hypothesis of Normality of errors is not verified.
For asymmetric distributions, the authors should be careful not to show in tables or figures symmetric error bars that would yield results that out of range (e.g. negative error rates).
If error bars reported in tables or plots, The authors should explain in the text how they were calculated and reference the corresponding figures or tables the text.
8. Experiments compute resources
For each experiment, does the paper provide sufficient information on the com- puter resources (type of compute workers, memory, time of execution) needed to reproduce the experiments?
Question:
Answer: [Yes]
Justification:
We describe those details in the appendix.
Guidelines:
The answer NA means that the paper does include experiments.
The paper should indicate the type of compute workers CPU or GPU, internal cluster, or cloud provider, including relevant memory and storage.
18

---

The paper should provide the amount of compute required for each of the individual experimental runs as well estimate the total compute.
The paper should disclose whether the full research project required more compute than the experiments reported in the paper (e.g., preliminary or failed experiments that didn't make into the paper).
9. Code of ethics
Question: Does the research conducted in the paper conform, in every respect, with the NeurIPS Code of Ethics https: httpp:///neurippe wwwpsurpp...pppupieeeee hccpppuiic/thhiccsueeeelleleelll
Answer: [Yes]
Justification: None
Guidelines:
The answer NA means that the authors have not reviewed the NeurIPS Code of Ethics.
If the authors answer No, they should explain the special circumstances that require a deviation from the Code of Ethics.
The authors should make sure to preserve anonymity (e.g.. if there special consid- eration due to laws or regulations in their jurisdiction).
10. Broader impacts
Question: Does the paper discuss both potential positive societal impacts and negative societal impacts of the work performed?
Answer: [Yes]
Justification: Yes, provide the Broader impacts in the paper.
Guidelines:
The answer NA means that there no societal impact of the work performed.
If the authors answer NA or No, they should explain why their work has no societal impact why the paper does not address societal impact.
Examples of negative societal impacts include potential malicious unintended uses (e.g., disinformation, generating fake profiles, surveillance), fairness considerations (e.g., deployment of technologies that could make decisions that unfairly impact specific groups), privacy considerations, and security considerations.
The conference expects that many papers will be foundational research and not tied to particular applications, let alone deployments. However, if there is a direct path to any negative applications, the authors should point it out. For example, it is legitimate to point out that an improvement in the quality of generative models could be used to generate deepfakes for disinformation. On the other hand, it is not needed to point that a generic algorithm for optimizing neural networks could enable people to train models that generate Deepfakes faster.
The authors should consider possible harms that could arise when the technology is being used as intended and functioning correctly, harms that could arise when the technology being used as intended but gives incorrect results, and harms following from (intentional or unintentional) misuse of the technology.
If there are negative societal impacts, the authors could also discuss possible mitigation strategies (e.g., gated release models, providing defenses in addition attacks, mechanisms for monitoring misuse, mechanisms to monitor how system learns from feedback over time, improving the efficiency and accessibility ML).
11. Safeguards
Question: Does the paper describe safeguards that have been put in place for responsible release of data or models that have a high risk for misuse (e.g., pretrained language models, image generators, scraped datasets)?
Answer: [NA]
Justification: Our work doesn't introduce any security concerns.
Guidelines:
The answer NA means that the paper poses such risks.
19

---

Released models that have high risk for misuse or dual use should be released with necessary safeguards to allow for controlled use of the model. for example by requiring that users adhere to usage guidelines or restrictions to access the model or implementing safety filters.
Datasets that have been scraped from the Internet could pose safety risks. The authors should describe how they avoided releasing unsafe images.
We recognize that providing effective safeguards challenging, and many papers do not require this, but we encourage authors to take this into account and make best faith effort.
12. Licenses for existing assets
Are the creators original owners of assets (e.g., code, data, models), used in the paper, properly credited and are the license and terms of use explicitly mentioned and properly respected?
Question:
Answer: [NA]
Justification: We don't use existing assets.
Guidelines:
The answer NA means that the paper does not use existing assets.
The authors should cite the original paper that produced the code package dataset.
The authors should state which version of the asset used and, if possible, include URL.
The name of the license (e.g., CC-BY 4.0) should be included for each asset.
For scraped data from a particular source (e.g., website), the copyright and terms of service of that source should be provided.
If assets are released, the license, copyright information, and terms of use in the package should be provided. For popular datasets, paperswi thcode. com datasets has curated licenses for some datasets. Their licensing guide can help determine the license of dataset.
For existing datasets that are re- re-packaged, both the original license and the license of the derived asset (if it has changed) should be provided.
If this information is not available online, the authors are encouraged to reach out to the asset's 's creators.
13. New assets
Question: Are new assets introduced in the paper well documented and is the documentation provided alongside the assets?
Answer: [NA]
Justification: paper does not release new assets
Guidelines:
The answer NA means that the paper does not release new assets.
Researchers should communicate the details of the dataset/code/model as part of their submissions via structured templates. This includes details about training, license, limitations, etc.
The paper should discuss whether and how consent obtained from people whose asset is used.
At submission time, remember to anonymize your assets (if applicable). You can either create an anonymized URL or include an anonymized zip file.
14. Crowdsourcing and research with human subjects
Question: For crowdsourcing experiments and research with human subjects, does the paper include the full text of instructions given to participants and screenshots, if applicable, as well details about compensation (if any)?
Answer: [NA]
Justification: the paper does not involve crowdsourcing nor research with human subjects.
20

---

Guidelines:
The answer NA means that the paper does not involve crowdsourcing nor research with human subjects.
Including this information the supplemental material is fine, but if the main contribu- tion of the paper involves human subjects, then as much detail as possible should be included the main paper.
According to the NeurIPS Code of Ethics, workers involved in data collection, curation, other labor should be paid at least the minimum wage in the country of the data collector.
15. Institutional review board (IRB) approvals or equivalent for research with human subjects
Question: Does the paper describe potential risks incurred by study participants, whether such risks were disclosed the subjects, and whether Institutional Review Board (IRB) approvals (or an equivalent approval/review based on the requirements of your country or institution) were obtained?
Answer: [NA]
Justification: the paper does not involve crowdsourcing nor research with human subjects.
Guidelines:
The answer NA means that the paper does not involve crowdsourcing nor research with human subjects.
Depending on the country in which research is conducted, IRB approval (or equivalent)
may be required for any human subjects research. If you obtained IRB approval, you should clearly state this in the paper.
We recognize that the procedures for this may vary significantly between institutions and locations, and expect authors to adhere to the NeurIPS Code of Ethics and the guidelines for their institution.
For initial submissions, do not include any information that would break anonymity (if applicable), such as the institution conducting the review.
16. Declaration of LLM usage
Question: Does the paper describe the usage of LLMs if it is an important, original, or non-standard component of the core methods this research? Note that if the L.M used only for writing, editing, or formatting purposes and does not impact the core methodology. scientific rigorousness, or originality the research, declaration is not required.
Answer: [NA]
Justification: the core method development in this research does not involve LLMs any important, original, or non-standard components.
Guidelines:
The answer NA
that the core method development in this research does involve LLMs as any important, original, or non-standard components.
Please refer to our LLM policy (https / //tp/r/psn...//eee for what should or should be described.
222L/LLLM
21

---

A Experimental details.
A.1 Further experiments.
Rollout steps
75
100
Baseline: : L
5.60e-04 (1.15e-04)
8.04e-02 (3.45e-02)
3.12e-01 (1.15e001)
6.19c-01 (2.08e-01)
9 (2.88e-01)
FNO [22) (4x Params.)
1.58-033
(2.82x004)
5.22e-02
(3.122002)
1.Hc-01 (9.78c-02)
4.22o-01
(1.95e011)
7.12=-01 (2.24c001
Spatial Hicrorchical
5.15c-04 (1.14-099
6.59c-02
(3,59-002)
2.69c-01
(1.31c-01)
5.77c-01 (2.39c-01)
9 (3 (3.08<001)
History Hierarchy
5.11e-04 (1,14e-04)
6.41e-02
(4.03e-02)
2.36e-01 (1.23e-01)
4.99c-01 (2.29;-01)
7.79c-01 (2.70e011)
2-step Abead
1.06e-03 (2.18e-04)
4.23e-02 (1.87e-02)
1.79e-01 (7.70e402)
4.11e.01 (1.69e-01)
7.17e-01 (2140e-01)
2-step History [19]
4.84e-04 (1)09e-04)
5.21e-02 (2,76e 02)
2.02e-01
(8.99e402)
4.42+-01
(1.86e01)
7.67e-01 (2.71e01)
3-step History [19]
4.53e-04 (9.65-085)
4.28c-02 (2.10e.02)
1.81e-01 (8.71e-02)
3.96e-01
(1.57c401)
6.81e-01 (2.39e-01)
Ours: F. 2
5.25e-04 (.255-0441.11e-004 (1))
4.02e-02 (1.92.-02)
1.61e-01 (6.50c-02)
3.73e-01 (1.51e-01)
6.55e-01 (2.44c001)
Ours: L 3
5.50c-04 (1,20c-04)
3.37e-02 (1.41e-02)
1,40-011 (6.16e-02)
3.24e-01 (1.18c001)
5,92e-01 (1.84c-01)
Table 6: Comparisons to other methods. We compare to other methods using mean squared error (MSE) for autoregressive roll-out across different lengths. We report the results using average and standard deviation, in parentheses, and demonstrate that ours (L = 3) outperforms others in all steps above -tsepp MSE. We also show that building models with both spatial and temporal hierarchies enhances the stability of the estimation.
Rollout steps
75
1000
Baseline: L=]
8.02e-07 (4.89e-07)
9.25e-06 (4.33e-06)
1.49c-05
(5.56e-06)
2.16c-05 (7.40x-06)
2.93c-05 (1.05e-45)
PNO (22) (4x Params.)
4.90e-07
(1.91e-77)
4.93e-06 (1.6le-06)
9.68e-06
(2.59e-06)
0.52e-05 (3.71e-06)
2.14e-05 (4.97e-06)
Sputial Hierarchy
6.28e-07 (3.96e-07)
1.03-05 (4.83 06)
1.87e055
(7.54e-06)
2.65e-05 (9.94g.06)
3.48e-05 05 (1.14e 05)
History-Hien Hierarchy
8.52e-07 (6.56ie(7)
8.95--06 (5. 06)
1.33e-05
(6.66c-06)
1.84e-05 (8.33e-06)
2.42e-05 50e-000
2-step Abead
0.89e-07 (4.10e-07)
5.88c-06
(3.29c-06)
9.94e-06
(4.42=006)
I 1.49e-05 (5.71e-06)
2.10e-05 (7.75e 06)
2-step History [19]
9.83e-07 (6, (6.22--07)
11c-06 (4.29
1.13c-05
(4.90c-06)
I 1.88--05 (5.87c-06)
2.18c-05 (? (6066-46)
3-stop History [19]
5.65c-07 (4.40c-07)
7.05c-06 (3.93e-06)
1.09c-05
(5.08c-06)
1.55c-05 (6.30c-06)
2.12c-05 (8.01e-06)
Ours: L=2
7.86e-07 (4.75c-07)
9.42e-06 (4. (4.91c-06)
1.24c-05 (5,59c-06)
1.62c-05 (6. (4-344))
2,10c-05 (6.24c-06)
Ours: L=3
6.37e-07 (3.80-00)
6.28c-06 (2.90s-06)
9.29e-06 (3.72e-06)
1.28e-05 (4.52e-46)
1.76e-05 (5.44c-06)
Table 7: Energy spectrum error comparison. Our L 3 model achieves significant reduction in
energy spectrum error during short-term rollouts (up to 200 steps), outperforming all comparison methods- including baseline with 4x more parameters for all prediction horizons beyond single- step forecasting.
Experimenting with FNO. All experiments thus far have utilized a UNet architecture with Fourier layers, as described in Section 5. To provide additional comparison, we run = with Fourier Neural Operator (FNO) [55]. Key adaptations include: (1) Using the e + fz loss (consistent with our other experiments in Table 6) instead of the default Soboley-norm objective [25], as the latter caused rapid prediction divergence (even within 200 steps): (2) For 256 x 256 resolution data, choose the Fourier mode number through grid search over {64, 96, 128 )2 (3) We use 3 FNO layers, which results in model with 130M parameters 4x larger than L = 3 configuration.
Model Efficiency. Our model utilizes the UNet's hierarchical framework to efficiently process multi-scale data. As shown in Table 8, compared to the baseline L = 1. ours L = 3 adds only few convolutional heads for handling and outputting latent variables z, resulting minimal parameter and inference time overhead.
Param Counts (M) Forward Time (seconds) 34.88
Baseline: L=
0.030
Ours: T. 3
36.67
0.031
Table 8: Model Running Time and Parameters comparison. On input with 256 x 256 resolution, our design (L- = 3) leverages the inherent hierarchical representation of UNet to process hierarchical latent variables, leading to minimal computational overhead than the baseline (I 1).
A.2 Evaluation metrics. Mean squared error (MSE). We use
22

---

MSE( Uh, ú,)) [u,. ü, 1
as our primary metric quantify the prediction error for short-term predictions.
Long-term stability. To evaluate the long-term stability of predictions, we leverage the system's s conserved energy, defined E { vy). As discussed in Section 1, calculate the standard deviation of the ground truth energy values. A trajectory is deemed stable if its energy remains within standard deviations this reference. This threshold informed by the observation that all training data lie within 4 standard deviations, while even the true dynamics extrapolated over ten times the training horizon (simulating extreme long-term behavior) do not exceed 5 standard deviations, as shown in Fig. 7. Appendix B.2 provides detailed visual comparison when the predicted dynamics exceed +5 standard deviation range.
1e3
Sampled Reference Truth
Energy Magitude
-5 stud
0
25000 50000
75000 100000 125000 150000 175000 200000 Rollout time steps
Figure 7: Energy of true dynamics along extended timescale. We compute the energy evolution of the true dynamics over 10-times the training dataset length (2 x 10 time steps). We show that energy along the sampled true dynamics remains within the +5 standard deviation of its mean.
Fourier spectrum error. We follow the implementation in py2d compute the energy spectrum' For short-term evaluation, we compute the mean absolute error of the energy spectrum the spatial FFT F(un), averaged over N timesteps:
N
F unl -(ūn]ll
N
10a), where a ground truth trajectory unavailable for every new initial condition, and the system loses memory of initial conditions showing invariant spectrum due to its ergodic properties, we instead compare against a fixed reference spectrum [[[[u]] == N E Flun], computed from an extremely long true trajectory. The error is then:
For long-term evaluation (e.g.,
Σ Fluml [F][u]]]]
N
==1
Zonal mean. The zonal mean of vorticity is computed by averaging vorticity perpendicular to jets (i.e., along the x-direction), shown as:
Un
NN
Un E Rdd denotes the vector of the vorticity æ-axis. N denotes the number discretization points.
23

---

A.3 Experimental setup.
Model Architecture We adopt the UNet architecture following the design in [53], with several customizations. For the experiments on 256 x 256 resolution data, our encoder consists of 5 groups with latent channel sizes of [16, 32, 64, 128, 128, 128). The spatial resolution is halved after each encoder group, resulting in 8 x resolution at the bottleneck for 256 x 256 input images. Each encoder group contains two convolutional blocks, which is made up of two convolutional layers with group normalization and residual connections. The decoder mirrors the encoder and incorporates skip connections from corresponding encoder layers. To improve the model's 's ability to capture long-range spatial dependencies and uulti-frrqqency signals, we add Fourier layer to each convolutional block. To balance between accuracy and computational efficiency, we limit the number of frequency modes to 96 for 256 x 256 data and apply Fourier convolutions in a depth- wise manner, i.e., without inter-channel communication.
For 512 x 512 resolution inputs, we use a similar architecture but add an extra encoder group, resulting in latent channel sizes of [16, 32, 64, 128, 128, 128). In the = setup, we set r² and , = 32, the same rates used for our 256 x 256 resolution experiments. Using these rates, we inject latent codes z(") and 22
(2) into encoder groups with feature resolutions of 64 x 64 and 32 x 32.
Data preprocessing. To stabilize training, we normalize the data using a constant scaling factor such that the resulting values have an approximate standard deviation of one. Specifically, for the dataset with Reynolds number 10000, we apply a dividing factor of 10, while for the dataset with Reynolds number 5000, we use dividing factor of 6.
Corner cases. For models that take multiple temporal frames as input, we simulate the initial rollout setting during training by randomly zeroing out early history frames with probability of 15%, approximating the absence of pre-initial frames.
Upsampling projections. Prior to injection, we upsample the latent variables and apply convolutional layers match the corresponding encoder channels. Decoding is performed at the same spatial resolutions as the injected latent codes. The same architectural setup as L = 3 is used to construct baseline models for Spatial Hierarchy and History Hierarchy. In the = case, we use r' 32 and inject x1( into the encoder group at the 64 64 resolution level.
Training details. We process the raw data generated from "py2d" [52] solver. Our emulator predicts vorticity at 0.05 time intervals, representing 500x coarser temporal resolution than the numerical solver's 's 0.0001 timestep. We normalize the simulation data to approximate standard normal distribution. Since the data is naturally zero-centered, we simply scale it by constant factors (10 for jet-containing flows and 6 for jet-free flows in Section 5.3) to achieve unit standard deviation.
For all experiments, we use the AdamW optimizer with learning rate at 10-4. We conduct experiments on NVIDIA A100, H100, and L40S GPUs.
B
Additional visualizations.
B.1
Further visualizations.
We provide additional visualization of short-term rollout in Fig. 8, Fig. 9, Fig. 10, and Fig. 11. The comparison across various datasets and methods show consistent improvements of our methods (L and L 3).
B.2 Stability.
We present four trials of long-term rollouts across all methods, each with distinct initial conditions. The energy evolution over 2 x 10 timesteps is shown in Figures 12 and 15, while the corresponding dynamic visualizations appear in Figures 13 through 17.
Our analysis reveals that deviations beyond the +5 standard deviations range of the ground truth energy distribution consistently correlate with unstable, exploding dynamics or overly smoothed, averaged patterns. While baseline methods exhibit persistent failures once instability occurs, our I, = 2 model demonstrates unique recovery pattern, where the dynamics (and associated energy values) can return to physically plausible states after temporary deviations. This robustness stems
24

---

GT 1125
L
L ==3 3
GT U50
L
1
== 3
425
14g6
1425
1600
4.50
Figure 8: We apply our approach flow dataset of Re 10 256 256 resolution without zonal jets. Our method (L = 3) gives more accurate predictions with lower associated residuals.
25

---

25-Step Prediction Difference
50-Step
75-Step
100-Step Prediction Difference
Prediction Difference
Prediction Difference
Ground Truth
==1
Spatial Hier.
History Hierr
2 2-Step Aleaad
2-Step History
3etep History
Figure Visualization of prediction and residual ground truth across methods and prediction steps. Our methods L 2, L = 3 consistently outperform all compared methods.
from the guidance provided by the top-level compressed variables, which help correct errors in fine-grained details.
For stability calculations, we conservatively classify a trajectory as unstable if it ever exceeds the predefined +5 standard deviations' bounds, regardless of subsequent recovery. This ensures fair comparison across methods.
26

---

25-Step Prediction Difference
50-Step
75-Step
100-Step Prediction Difference
Prediction Difference
Prediction Difference
Ground Truth
==1
Sputial Hier
History Hier.
2-Step Alead
2-Step History
3-sep History
Figure 10: Visualization of prediction and residual to ground truth across methods and prediction steps. Our methods L = 2, L = 3 consistently outperform all compared methods.
27

---

25-Step
50-Step
75-Step
100-Step Prediction Difference
Prediction Difference
Prediction Difference
Prediction Difference
Ground Truth
==1
Spatial Hier
History Hier.
2-Step Alead
2-Step History
3-sepp History
Figure 11: Visualization of prediction and residual to ground truth across methods and prediction steps. Our methods L = 2, L = 3 consistently outperform all compared methods.
28

---

Trial
Trial 2
1e2
Sampled
le2
Sampled Predictions
Eneogy Maginude
Enengy Magitudk
1=]
25000
50000 75000 100000 125000 150000 175000 200000
25000
50000
75000 100000 125000 159000 175000 200000
Sampled Predictions
1n3
Sampled Preaictions
Sputial Riee
Envegy Magitaal
Energy Magitadle
25000
50000
75000
100000 325006
55000
375000 200000
25000
S0000
75000
100000
125000
15000 375000 200000
1e3
Sampled Precictions
1e3
Sampled Precictions
History Hilrr.
Enwwgy Magiude
Enowyy Magitade
25000 50000 75000
100000
325000
150000 175000 200000
25000
50000
75000
100000
125000
150000
175000
200000
Sampled Predictions
Sampled Predictions
History Hier.
Energy Magibusl
Erergy Magibusla
175000
Sampled Presictions
Sampled Presictions
2-step History
Enang Magitate
Enangy Magitutle
50000
50000
1o2
Sampled Predictions
3elep History
Energy Huggh flgguude
Enwrgy Ma Rugiuude
25000
50000
75000
100000
325000
150000
375000 200000
25000
50000
75000
100000
125000
150000
175000
200000
1n3
Sampled
1+5
Sampled
Ours L=22
Erwroy Magiuude
Erwrgp Maginude
50000 25000
75000 100000
125000 15000D 375000 200080
50000 25000
75000 100000
125000 15000D 375000 200000
A03
Sampled Predictions
103
Sampled Predictions
Ours: L-3
Spargy Magitule
Snargy Maginude
25000 50000 75000 100000 125000 150000 175000 200000 Railout time steps i
25000 50000 75000
100000 125000 150000 175000 200000 Railout time steps
Figure 12: Energy V.S. rollout time steps the predicted dynamics. The dashed horizontal line shows the maximum and minimum values computed using standard deviations from the mean of true dynamics. Our L method able to maintain energy along the long-term predictions.
29

---

Step:
25k
50k
75k
100k
125k
150k
175k
Sputial Hier.
History Hier.
2-step Ahend
2-step Histury
dan Hastryy A
OunsL=22
Ours:La3
Figure 13: Long-Term Rollout Visualization. Long-tern term dynamics corresponding to Trial 1 in Fig. 12. The rightmost column shows the frame with the maximum deviation-exeee exceeding +5 standard deviations from the reference truth statistics, with a blank block indicating that all states along the sampled trajectory remain within this range. When energy predictions fall outside the reference distribution, non-physical dynamics emerge texploding/averaging artifacts; blank blocks indicate predictions within the +5 std. range). Our L 3 model demonstrates superior stability across all comparisons, while exhibits deviations correlated with energy prediction errors.
30

---

Step: 1
25k
50k
75k
100k
125k
150k
175k
Most Devinted
Sputial Hier.
History Hiee..
2-tep Altead
2-sepp History
3-step History
Ours:L=2
Ours:e3
Figure 14: Visualization long-term rollout. We visualize the long-tern dynamics corresponds to trial 2 in Fig. 12. The rightmost column shows the frame with the maximum deviation- exceeding +5 standard deviations from the reference truth statisties, with blank block indicating that all states along the sampled trajectory remain within this range. The results show that when the energy is deviated from the reference truth distribution, the dynamies exhibit exploding or averaging non- physical dynamics. Among all comparison methods, our L = 3 gives the most stable predictions.
31

---

Trial
Trial
162
Sampled
le2
Sampled
Enwigy Maginude
Enengy Magitudk
1=1
25000
50000 75000 100000 125000 150000 175000 200000
25000
50000
75000 100000 125000 159000 175000 200000
1-3
Sampled Predictions
1+3
Sampled Preaictions
Sputial Hier.
Envegy Magitatee
Enaegy Magitadle
25000
50000 75000 160000 125006 150000 375300 200000
25000
S0000
75000
100000 125000
15000 375000 200000
1n3
Sampled Precictions
ln1
Sampled Predictions
History Hilrr.
Enwnyy Magiude
Enowyy Magitad
25000 50000 75000 100000 325000 150000 175000 200000
25000
50000
75000
100000
125000
150000
175000
200000
Sampled Predictions
Sampled Predictions
History Hier.
Energy Magibusl
Energy Magibusl
Sampled Presictions
Sampled Precictions
2-step History
Energy Magitate
Snargy Magitutle
50000
50000
102
Sampled Predictions
3elep History
Enwagy Hu gluude
Enwrgy Ma Ruugtude
25000
50000
75000
100000
125000
150000
175000 200000
25000
50000
75000
100000
125000
150000
175000
200000
1+3
Sampled
1+3
Sampled
Ours L=22
Erwroy Magiuude
Erwrgy Maginude
25000
50000
75000 100000 125000 15000D 375000 203080
50000 25000
75000 100000
125000
15000D 375000 200000
A03
Sampled Predictions
103
Sampled Predictions
Ours: L-3
Spargy Magitule
Snargy Maginude
25000 50000 75000 100000 125000 150000 175000 200000 Railout time steps i
25000 50000
75000
100000 125000 150000 175000 200000 Railout time steps
Figure 15: Energy V.S. rollout time steps the predicted dynamics. The dashed horizontal line shows the maximum and minimum values computed using standard deviations from the mean of true dynamics. Our L method able to maintain energy along the long-term predictions.
32

---

Step: 1
25k
50k
75k
100k
125k
150k
175k
Most Devinted
Sputial Hier.
History Hierr.
2-tep Altead
2-sepp History
3-step History
Ours:L=2
Ours:=33
Figure 16: Visualization long-term rollout. We visualize the long-term term-temm dynamics corresponds to trial 3 in Fig. 15. The rightmost column shows the frame with the maximum cettaatinn-xxeeee +5 standard deviations from the reference truth statisties, with blank block indicating that all states along the sampled trajectory remain within this range. The results show that when the energy is deviated from the reference truth distribution, the dynamies exhibit exploding or averaging non- physical dynamics. Among all comparison methods, our L = 3 gives the most stable predictions.
33

---

Step: 1
25k
50k
75k
100k
125k
150k
175k
Most Devinted
Sputial Hierr
History Hier.
2-tep Altead
2-sepp History
3-step History
Our Ours:Laa =2
Ours:e3
Figure 17: Visualization long-term rollout. We visualize the long-temn dynamics corresponds to trial 4 in Fig. 15. The rightmost column shows the frame with the maximum cettaatinn-xxeeeee +5 standard deviations from the reference truth statisties, with blank block indicating that all states along the sampled trajectory remain within this range. The results show that when the energy is deviated from the reference truth distribution, the dynamies exhibit exploding or averaging non- physical dynamics. Among all comparison methods, our L = 3 gives the most stable predictions.
34


Reviews:
[{"review_id": "9QDKR9mloP", "reviewer": "Reviewer_CHmV", "summary": "Authors propose a novel autoregressive multiscale architecture for stable and accurate integration of partial differential equations. Consider a physical field with several resolutions $\\phi_h(t), \\phi_{2h}(t), \\phi_{3h}(t)$ where subscript is a grid spacing. The approach of the authors is, roughly, to predict $\\phi_h(t+\\Delta t), \\phi_{2h}(t + 2\\Delta t), \\phi_{3h}(t + 3\\Delta t)$ from $\\phi_h(t), \\phi_{2h}(t + \\Delta t), \\phi_{3h}(t + 2\\Delta t)$. This is done by incorporating fields with low resolution into the hidden layers of UNet architecture.\n\nNumerical experiments indicate that implicit neural emulators proposed by authors lead to improved stability, improved accuracy for longer rollouts and better statistical measures (energy spectrum, zonal mean).", "questions": "**Main questions:**\n1. Suppose we have the scheme with three levels $u\\_{n}, z^{(1)}\\_{n+1}, z^{(2)}\\_{n+2}, z^{(3)}\\_{n+3}$. To start iterations we need $u\\_0$ (known) and $z^{(1)}\\_{1}, z^{(2)\\}_{2}, z^{(3)}\\_{3}$ (unknown). How coarser state from future time steps are supplied on the first step?\n\n2. A flow of information looks a little peculiar. Why is it important/beneficial to use coarser information from the future to predict the state of the system? For example, to predict $u\\_{n+1}$ authors suggest to use not only $z^{(1)}\\_{n+1}$ which seems natural but also $z^{(2)}\\_{n+2}$, $z^{(3)}\\_{n+3}$, etc, which should not be directly related to $u\\_{n+1}$. Can the authors please comment on that choice?\n\n3. Related to the previous question. Coarser states are predicted but sometimes not used directly for prediction of corresponding fine state. For example, on step $n-1$, network predicts $u_{n}$ (and also $z^{(1)}\\_{n+1}, z^{(2)}\\_{n+2}, z^{(3)}\\_{n+3}$) using $u\\_{n-1}, z^{(1)}\\_{n}, z^{(2)}\\_{n+1}, z^{(3)}\\_{n+2}$ but on previous two time steps $z^{(2)}\\_n$ and $z^{(3)}\\_n$ were predicted. Presumably, these state are directly related to $u_n$ but not used in its prediction. Did the authors try to evaluate the architecture that uses $u\\_{n}, z^{(1)}\\_{n+1}, z^{(2)}\\_{n+1}, z^{(3)}\\_{n+1}$ to predict $u\\_{n+1}$ along with coarse latents for future time?\n\n4. Authors show results for $L=1$, $L=2$, $L=3$. Did authors try $L=4$, $L=5$ and so forth? What do authors expect if one keeps increasing the number of coarse levels?\n\n5. Once tacit assumption is that it is easier to predict coarse state for longer lead times. Can the authors support this assumption by theoretical or empirical evidence? If one runs $L=1$ on a coarser version of data is it expected that stability of rollout improves?\n\n6. (Figure 4) Why is the error of the energy spectrum for the $L=3$ case larger than for $L=2$? Simultaneously $L=3$ shows better stability than $L=2$. Similarly, for zonal mean the results seem to indicate that $L=2$ surpasses $L=3$. These results in combination make it slightly unclear whether the increase in the number of levels is beneficial or not. Can the authors please comment on that?\n\n**Minor issues**\n1. (line 1, abstract) \"Neural PDE solvers offer a powerful tool for modeling complex dynamical systems, ...\" I am not sure if modeling is an accurate word in this context. It seems to me solvers are used to approximate a solution when the model is already fixed.\n\n2. (lines 5-8, abstract) \"Drawing inspiration from the stability properties of numerical implicit time-stepping methods, our approach leverages predictions several steps ahead in time at increasing compression rates for next-timestep refinements.\" Sounds ambiguous. Please consider adding \"we developed an approach that\" after the first comma or replacing \"Drawing inspiration from\" with \"Inspired by\".\n\n3. (line 12, abstract) \"..., significantly outperforming autoregressive baselines, ...\" The method by the authors is autoregressive too, so I would suggest adding the clarification, e.g., \"standard autoregressive baselines\" or \"non-hierarchical autoregressive baselines\".\n\n4. (line 102, Section 3) \"Consider a nonlinear dynamical system ...\" The proposed method uses standard geometrical coarsening so it is suitable for PDEs on spatial domains. All examples the authors provide are for Navier-Stokes equation (equation (13)). I suggest adding to equation (1) some partial derivatives with respect to spatial coordinates. Otherwise it is not clear what is the fundamental difference between (1) and (2) and why one needs a partial derivative in (1). Besides that, line 104 contains an explanation of $\\boldsymbol{u}_{t}$ that was never used in (1) and only later appears in (2).\n\n5. (line 212, Section 5) \"... (see details in Appendix A.2)\" Should be Appendix A.3.\n\n6. (line 255, Section 5.1) It is not explained what the zonal mean is. Please, consider providing reference to the appendix where zonal mean is defined.", "limitations": "yes", "rating": 6, "confidence": 4, "soundness": null, "presentation": null, "contribution": null, "strengths": "", "weaknesses": "", "quality": 4, "clarity": 4, "significance": 4, "originality": 4, "strengths_and_weaknesses": "**Strengths**\nAuthors suggest an interesting original research idea and thoroughly confirm it on a set of reasonable experiments. Overall, the principle idea is clearly explained and central claims are well-supported.\n\n**Weaknesses**\nMostly minor typos, unclear parts and several undiscussed issues.\n\nSee the section below for details.", "was_revised": true, "final_justification": "My view does not change after rebuttal, so I maintain my original score.\n\nI believe that the authors propose highly original idea that can be useful for many domains. The proposed architecture modifications are supported by thoughtful experimental study with many ablations and comparison with related techniques.", "rebuttal": "We greatly appreciate your feedback and your insightful comments on our work! We provide our response as follows: \n\n### Q1 \n> Suppose we have the scheme with three levels $z_{n+1}^{(1)}, z_{n+2}^{(2)}, z_{n+3}^{(3)}$. To start iterations we need $u_0$ (known) and $z_{1}^{(1)}, z_{2}^{(2)}, z_{3}^{(3)}$ (unknown). How coarser state from future time steps are supplied on the first step?\n\nDetails on handling the initial cases were briefly provided in Appendix A.3 (see the \"Corner Cases\" subsection).\nSpecifically, during training, with a small probability, we allow the model to receive partially empty latent variables as input and predict the missing ones in the hierarchy.\nFor example, in the $L=2$ setting, with a 15% probability, we ask the model to take $[u_n, \\mathbf{0}]$ as input and output $z_{n+1}^{(1)}$. During evaluation, we only need a single extra forward run to produce the full latent state $[u_n, \\hat{z}_{n+1}^{(1)}]$ to continue the autoregressive rollout.\nTherefore, all methods are evaluated using only a single-frame initial condition for fair comparison. As an example, to generate a sequence trajectory of length $T=1000$, our $L=3$ model performs 1002 forward passes, where the extra $L - 2 = 3$ passes handle these corner cases.\nWe will update this discussion in the revised version!\n\n### Q2\n> A flow of information looks a little peculiar. Why is it important/beneficial to use coarser information from the future to predict the state of the system? For example, to predict $u_{n+1}$ authors suggest to use not only \n$z_{n+1}^{(1)}$ which seems natural but also $z_{n+2}^{(2)}, z_{n+3}^{(3)}$, etc, which should not be directly related to $u_{n+1}$. Can the authors please comment on that choice?\n\nOur key innovation lies in leveraging future information, for example, $z_{n+2}^{(2)}, z_{n+3}^{(3)}$ in $L=3$ method. \nThis design can be interpreted in two complementary ways: \n(1) As discussed in Section 4.2, it aligns with the philosophy of multi-step implicit methods. Though not in its canonical form for steps  > 2. As an instance, a rearranged two-step Adams–Moulton $u_{n+1} \\approx  u_{n+2} - \\Delta t \\left(\\frac{5}{12} f(u_{n+2}) + \\frac{8}{12} f(u_{n+1})-\\frac{1}{12}f(u_{n})\\right)$ illustrates how future terms contribute to current prediction. And the iterative refinement of $\\hat{u}_{n+1}$ through the chain of $\\hat{z}_{n+1}^{(1)}$ and  $\\hat{z}_{n+1}^{(2)}$ resembles the root-finding algorithms of the implicit solvers.\n(2) From a temporal modeling standpoint, the use of long-range guidance from $z_{n+2}^{(2)}, z_{n+3}^{(3)}$ helps stabilize $\\hat{u}_{n+1}$ and improves temporal consistency in long-term forecasts. And with the abstraction strategy of downsampling, our method balances between the model's stability and its capacity to look multi-step ahead.\n\n### Q3\n\n> Related to the previous question. Coarser states are predicted but sometimes not used directly for prediction of corresponding fine state. For example, on step $n-1$, network predicts $u_n$ (and also $z_{n+1}^{(1)}, z_{n+2}^{(2)}, z_{n+3}^{(3)}$) using $u_{n-1}, z_n^{(1)}, z_{n+1}^{(2)}, z_{n+2}^{(3)}$ but on previous two time steps $z_n^{(2)}$ and $z_n^{(3)}$ were predicted. Presumably, these state are directly related to $u_n$ but not used in its prediction. Did the authors try to evaluate the architecture that uses $u_n, z_{n+1}^{(1)}, z_{n+1}^{(2)}, z_{n+1}^{(3)}$ to predict $u_{n+1}$ along with coarse latents for future time?\n\nThank you for suggesting this new baseline experiment. Compared to our design, where latent variables $z_{n+1}^{(1)}, z_{n+2}^{(2)}, z_{n+3}^{(3)}$ are distributed across future time steps, the configuration using $z_{n+1}^{(1)}, z_{n+1}^{(2)}, z_{n+1}^{(3)}$ concentrates on all latent variables at the same time step to predict target $u_{n+1}$. \n\nThis configuration would likely make the latent states more directly correlated with $u_{n+1}$, potentially simplifying learning. However, we expect that such a design would converge toward would converge into the configuration resembling spatial hierarchy design evaluated in Table 1, with input on $u_n, z_{n+1}^{(1)}, z_{n+1}^{(2)}, z_{n+1}^{(3)}$, where our empirical results find that the hierarchical variables of the same time frame only has minimal improvements compared to the $L=1$ baseline.\n\n### Q4\n>Authors show results for $L=1, L=2, L=3$. Did authors try $L=4,L=5$ and so forth? What do authors expect if one keeps increasing the number of coarse levels?\n\nThank you for the suggestion! We conducted additional experiments with \n$𝐿=4$ on the 512×512 resolution turbulent flow dataset. Specifically, we added an extra coarse level $z^{(4)}$ corresponding to an 8×8 resolution, on top of our $L=3$ model. We followed the same training and evaluation protocols as used for $L=3$. The results are summarized below:\n\n| Model | 1 Step | 10th Step | 25th Step | 50th Step | \n|-------|---------|-----------|-----------|-----------|\n| L=3   | 1.41e-03 (3.07e-04) | 1.56e-02 (9.92e-03) | 5.75e-02 (3.04e-02) | **1.81e-01** (9.02e-02) | \n| L=4   | **1.36e-03** (2.47e-04) | **1.40e-02** (4.62e-03) | **5.72e-02** (2.10e-02) | 2.20e-01 (7.93e-02) |\n\nTable 1: Comparison between $L=3$ and $L=4$ for $512\\times512$ data.\n\nWe find that the $L=4$ model performs slightly better on short-term predictions (≤ 25 steps) but achieves similar results on longer rollouts, when compared to $L=3$. We suspect that the 8×8 resolution may contain too limited information for a $512 \\times 512$ input for further significant improvements.\nTo conclude, we think that improving the long-term performance of deeper hierarchies may require revisiting the resolution choices of intermediate levels $z^{(1)}, z^{(2)}, z^{(3)}$ to ensure more balanced information flow across scales. And a more complex dataset might lead to more choices and a higher hierarchy.\n\n### Q5\n>Once tacit assumption is that it is easier to predict coarse state for longer lead times. Can the authors support this assumption by theoretical or empirical evidence? If one runs $L=1$ on a coarser version of data is it expected that stability of rollout improves?\n\nWe thank you for suggesting this ablation experiment! Yes! We agree that the ease of predicting coarse states, especially for $L$-step ahead, could partly contribute to our long-term stability.\nOne related example was the observation we made in Appendix B.2 regarding a unique recovery pattern seen in our $L=2$ model.\nWe will consider adding this ablation study in the final version of the submission.\n\n ### Q6\n> (Figure 4) Why is the error of the energy spectrum for the $L=3$ case larger than for $L=2$? Simultaneously, $L=3$ shows better stability than $L=2$. Similarly, for the zonal mean, the results seem to indicate that $L=2$ surpasses $L=3$. These results, in combination, make it slightly unclear whether the increase in the number of levels is beneficial or not. Can the authors please comment on that?\n\nFigure 4 primarily benchmarks long-term rollout stability. Since the MSE metric is less informative in this setting, we instead assess the model's ability to preserve physical properties. Given the extended rollout horizon, the variance is significant (see the shaded regions in Figure 4(a)). Both $L=2$ and $L=3$ models show good alignment with the ground truth and perform substantially better than the baseline methods. In this figure, the comparison between $L=2$ and $L=3$ is less critical than their overall improvement over the baselines. In contrast, Figure 3(b) offers a more stringent criterion for long-term evaluation, under which the performance gap between $L=2$ and $L=3$ becomes more pronounced."}, {"review_id": "4Edf3fGJ4D", "reviewer": "Reviewer_PLCM", "summary": "This paper introduces a neural PDE solver, with a claim of the implcit manner, leading to a method for a sliding window prediction. The results presented show the effectiveness of the proposed method. The hierarchical design is not news but play a great role in the method. The analysis of the method with the energy spectrum is insightful and interesting (at least for me).", "questions": "1. I wonder where the GT is from. Maybe the analytic solution should be listed for the canonical flow to provide a good view for the readers.\n2. Some repeated experiments for the results to provide the bar of the curves.\n3. how the latents $z$ should be, I think some more theoretical insights for more complex standard PDEs should be provided.\n4. how the disantce metirc $d$ infleuence the problem, what if it is not the l1/l2 norm", "limitations": "N/A", "rating": 4, "confidence": 3, "soundness": null, "presentation": null, "contribution": null, "strengths": "", "weaknesses": "", "quality": 2, "clarity": 2, "significance": 2, "originality": 2, "strengths_and_weaknesses": "Pros:\n1. the analysis for the NS results are great, maybe it is more great if it will be visualized in a sphere.\n2. the method shows some insight in the implcit manner in the Neural PDEs.\n\nCons:\n1. the latent state for $z$ can not be theoretically guaranteed. Therefore for some PDEs, there might be not latent state existing to work as the explicit prediction manner.\n2. I am not sure for the training time. I guess from the loss function design it should take long.\n3. the resolution for the results are not large to fully show the effectiveness of the proposed method.", "was_revised": true, "final_justification": "Most of my problems about the empirical setting are resolved, including the metric, the cases for the GT reconstruction and training time. But I still concern abut the method from the theoretical prespective. Therefore I just raise my score to 4.", "rebuttal": "We appreciate your feedback!\n\n### W1\n\n> the latent state for $z$ can not be theoretically guaranteed. Therefore for some PDEs, there might be not latent state existing to work as the explicit prediction manner.\n\nThe notion of the latent space $z$ in our framework is intentionally designed to be flexible, as long as it provides a compressed representation of the original data. While more sophisticated latent spaces—such as semantic-aware representations for video data—can be beneficial for specific domains, we adopt a deliberately simple design where $z$ corresponds to spatially downsampled inputs. In the context of PDEs, this form of representation is quite general and has been widely used, for instance, in multigrid methods. We demonstrate that even with this generic choice of latent space, our model achieves significant improvements over strong baselines.\n\n\n### W2\n\n>I am not sure for the training time. I guess from the loss function design it should take long.\n\nAs shown in Table 6 (Appendix A.1), the configuration with $L=3$ exhibits nearly identical runtime and parameter count compared to the baseline with $L=1$. This efficiency arises from our model’s ability to handle multi-resolution data effectively. Specifically, the multi-level architecture reuses most of the U-Net components from the single-level baseline, and we exploit the inherent multi-resolution structure of the U-Net to process inputs at their corresponding feature resolutions. Please refer to Section 5 for additional implementation details. \n\nBenefited from our design, the model only requires a single forward pass per training iteration. With H100 GPU and training mode with batch size 32, single iteration for $L=1$ takes 0.2067 (0.0017) second and $L=3$ takes 0.2204 (0.0021) second. This demonstrates that our method introduces minimal overhead while enabling more expressive multi-level representations.\n\n| Model                        |Seconds per training step| 1 Step                | 10th Step            | 25th Step            | 50th Step            | 75th Step            | 100th Step           |\n|------------------------------|------------------------------|-----------------------|----------------------|----------------------|----------------------|----------------------|----------------------|\n| L=1                 | **0.2067** (0.0017) | 9.17e-02 (7.96e-03)  | 3.87e-01 (6.17e-02)  | 8.21e-01 (1.16e-01)  | 1.17e+00 (8.22e-02)  | 1.29e+00 (6.15e-02)  | 1.34e+00 (5.62e-02)  |\n| L=2  (ours)                        |0.2134 (0.0017) | **5.25e-04** (1.15e-04)  | 8.11e-03 (4.15e-03)  | 4.02e-02 (1.92e-02)  | 1.61e-01 (6.50e-02)  | 3.73e-01 (1.51e-01)  | 6.55e-01 (2.44e-01)  |\n| L=3  (ours)                        |0.2204 (0.0021)| 5.50e-04 (1.20e-04)  | 6.59e-03 (2.50e-03)  | **3.37e-02** (1.41e-02)  | **1.40e-01** (6.16e-02)  | **3.24e-01** (1.18e-01)  | **5.92e-01** (1.84e-01)  |\n \nTable 1: Comparison of training time and short-term accuracy.\n\n### W3\n> the resolution for the results are not large to fully show the effectiveness of the proposed method.\n\nOur experiments are conducted on a highly chaotic turbulent flow dataset with a Reynolds number of $10^4$, where we evaluate long-term stability over an extended horizon of more than $2 \\times 10^5$ rollout steps. Setting up experiments at a resolution of $256 \\times 256$ is already extremely challenging due to the complexity of the dynamics. To further demonstrate the robustness and effectiveness of our approach, we also conduct experiments on higher-resolution data at $512 \\times 512$, showing that our proposed method scales well and maintains strong performance in this more demanding setting.\n\n### Q1\n> I wonder where the GT is from. Maybe the analytic solution should be listed for the canonical flow to provide a good view for the readers.\n\nWe use an off-the-shelf numerical solver to generate the ground truth data; please refer to Section 5 for implementation details. Our experiments are conducted on highly challenging cases, many of which, such as those governed by the Navier–Stokes equations, do not have known closed-form solutions. \n\n### Q2\n> Some repeated experiments for the results to provide the bar of the curves.\n\n\nWe report standard deviations for most of our results to capture the variability across predicted trajectories—for example, in Figure 3(a), Figure 4(a, c), and Table 2 in the appendix.  We acknowledge the importance of statistical robustness and will include results from additional repeated runs in the final version to further strengthen our findings.\n\n### Q3\n> how the latents $z$ should be, I think some more theoretical insights for more complex standard PDEs should be provided.\n\nAn effective latent space $z$ should provide a compressed yet informative representation of the data, making it easier to estimate than directly predicting the original data, while still retaining sufficient information for accurate prediction. This principle is closely related to several theoretical frameworks, such as the information bottleneck and predictive coding. We will update the final version of the paper to include a discussion of these connections.\n\n### Q4\n> how the disantce metirc $d$ infleuence the problem, what if it is not the l1/l2 norm\n\nIn addition to the standard $l1$ and $l2$ losses, we here provide results for experiments with the Sobolev loss, proposed in [r1].\nWe show across different training objectives, our $L=3$ method consistently yields significant improvements compared to the $L=1$ baseline.\n\n| Model                        | 1 Step                | 10th Step            | 25th Step            | 50th Step            | 75th Step            | 100th Step           |\n|------------------------------|-----------------------|----------------------|----------------------|----------------------|----------------------|----------------------|\n| L=1 ($l1+l2$ loss)                         | 5.60e-04 (1.15e-04)  | 1.29e-02 (4.87e-03)  | 8.04e-02 (3.45e-02)  | 3.12e-01 (1.15e-01)  | 6.19e-01 (2.08e-01)  | 9.38e-01 (2.88e-01)  |\n| L=1 (Sobolev loss)           | 7.01e-04 (4.07e-04)  | 1.64e-01 (8.54e-02)  | 8.76e-01 (2.45e-01)  | 1.40e+00 (2.30e-01)  | 1.64e+00 (2.56e-01)  | 1.69e+00 (2.51e-01)  |\n| L=3  ($l1+l2$ loss)                       | **5.50e-04 (1.20e-04)** | **6.59e-03 (2.50e-03)** |**3.37e-02 (1.41e-02)**  | **1.40e-01 (6.16e-02)**  | **3.24e-01 (1.18e-01)**  | **5.92e-01 (1.84e-01)**  |\n| L=3 (Sobolev loss)           | 6.49e-04 (3.03e-04) | 1.45e-02 (6.35e-03) | 7.28e-02 (3.80e-02) | 2.64e-01 (1.42e-01) | 5.74e-01 (2.78e-01) | 8.83e-01 (3.21e-01)|\n\nTable 2: Comparison of short-term accuracy when trained with different metrics\n\n| Model                        | 1k steps               | 10k steps           | 100k steps            |\n|------------------------------|-----------------------|----------------------|----------------------| \n| L=1 (Sobolev loss)           | 3% | 0% | 0% |\n| L=3 (Sobolev loss)           | **100%** | **99%** | **97%** |\n\nTable 3: Comparison of long-term stability when trained using the Sobolev norm\n\n[r1] Learning Dissipative Dynamics in Chaotic Systems, Zongyi Li, Miguel Liu-Schiaffini, Nikola Kovachki, Burigede Liu, Kamyar Azizzadenesheli, Kaushik Bhattacharya, Andrew Stuart, Anima Anandkumar"}, {"review_id": "qu2veYTphw", "reviewer": "Reviewer_HWeY", "summary": "This paper introduces a Multiscale Implicit Neural Emulator to address the long-term instability and error accumulation problems common in neural PDE solvers. Inspired by the stability of implicit numerical time-stepping methods, the proposed approach leverages a hierarchy of coarse-grained future state predictions to guide next-step forecasting. By conditioning on compressed representations of multiple future time steps, the model enforces long-range temporal coherence and significantly improves stability over extended rollouts. Finally authors demonstrated performance of the proposed method on 2D turbulent fluid simulations, achieving better performance than autoregressive baselines.", "questions": "Have you evaluated the performance of your method on other types of PDEs?", "limitations": "Yes.", "rating": 3, "confidence": 3, "soundness": null, "presentation": null, "contribution": null, "strengths": "", "weaknesses": "", "quality": 3, "clarity": 3, "significance": 2, "originality": 3, "strengths_and_weaknesses": "Strengths:\n\n・In my understanding, implicit methods have been a cornerstone of success in numerical analysis. Incorporating such techniques into neural PDE solvers could open new directions for advancement in this field.\n\n・The overall presentation is clear and well-written.\n\nWeaknesses:\n\n・In the experiments, the method is only evaluated on a single PDE — the Navier-Stokes equations. Investigating its scalability and performance on other types of PDEs would strengthen the validation of the approach.\n\n・In numerical analysis, the superiority of implicit methods over explicit ones is supported by various theoretical insights (e.g., order of accuracy, stability criteria). It is not discussed whether such theoretical understanding can be carried over or adapted to the neural setting. Addressing this question would significantly strengthen the technical foundation of the proposed method.", "was_revised": true, "final_justification": "While the result is good, the rebuttal does not solve my concerns enough, especially, concern about technical soundness. I have fix my score 3.", "rebuttal": "We appreciate your feedback!\n\n### W1 and Q1\n> In the experiments, the method is only evaluated on a single PDE — the Navier-Stokes equations. Investigating its scalability and performance on other types of PDEs would strengthen the validation of the approach.\n\nWe would like to emphasize that emulating the Navier-Stokes equations under high Reynolds numbers (e.g., $10^4$), with high-resolution data ($256 \\times 256$ and $512 \\times 512$), and over very long rollout horizons (up to $2 \\times 10^5$ steps) presents an extremely challenging task. Achieving substantial improvements in such a demanding setting is notably difficult and highlights the strength of our approach.\n\nFor reference, we have also conducted experiments on a relatively simpler 1D Kuramoto–Sivashinsky PDE. Results are shown below:\n| L   | 1-step             | 5-step             | 10-step            | 15-step            | 20-step            | 25-step            | 50-step            | \n|-----|--------------------|--------------------|--------------------|--------------------|--------------------|--------------------|--------------------|\n| 1   | 5.43e-06 (1.67e-06) | 2.11e-04 (6.25e-05) | 4.16e-03 (4.30e-03) | 3.88e-02 (4.31e-02) | 1.44e-01 (1.02e-01) | **4.07e-01** (1.76e-01) | **2.03e+00** (1.68e-01) | \n| 3   | **5.21e-06** (1.51e-06) | **1.31e-04** (6.06e-05) | **3.78e-03** (5.09e-03) | **3.28e-02** (4.04e-02) | **1.40e-01** (9.70e-02) | **4.03e-01** (1.88e-01) | **2.00e+00** (1.88e-01) | \n\nTable: Comparison of mean squared error (MSE) on Kuramoto-Sivashinsky\n\n### W2\n\n> In numerical analysis, the superiority of implicit methods over explicit ones is supported by various theoretical insights (e.g., order of accuracy, stability criteria). It is not discussed whether such theoretical understanding can be carried over or adapted to the neural setting. Addressing this question would significantly strengthen the technical foundation of the proposed method.\n\nOur design is inspired by implicit methods and shares a conceptual similarity: to produce the final estimate $\\hat{u}_n$, the model performs multiple stages of estimation and refinement through the latent variables $\\hat{z}^{(L-1)}_n$, ..., $\\hat{z}^{(1)}_n$. This process loosely approximates the behavior of an implicit numerical solver that iteratively refines solutions via root-finding steps. However, since these operations involve transitions across multiple latent spaces, whose structure may not admit a simple closed-form, in contrast to classical implicit solvers, it becomes challenging to derive a rigorous theoretical analysis. We hope that the strong empirical performance demonstrated by our method will motivate future work toward a deeper theoretical understanding of this framework."}, {"review_id": "71koEdUaIo", "reviewer": "Reviewer_vohk", "summary": "The authors present a new method for predicting downsampled future states along with the full-resolution, next state. Additionally, the architecture takes in the full-resolution current state, and a downsampled prediction of a future state. The method performs well on the given benchmarks in short and long term rollouts.", "questions": "- It would be interesting to see how the improvement of the implicit scheme varies with the dynamics of the PDE system. For simple systems, explicit schemes are often good enough so would the performance gain be less?\n\n- In Figure 6, the GT for the jet system (top row) for time 25 and 50 look nearly identical (u_25, u_50). Is this the desired behavior?", "limitations": "The authors don’t seem to present any limitations of their work, despite claiming to do so in the conclusion. First, the method seems to be constrained to only Unets, since it doesn’t seem clear how to inject hierarchical information into other architectures. Since the method is constrained to Unets, it doesn’t seem applicable to mesh-based problems. In addition, mesh-based problems don’t admit a consistent downsampling operator for evaluating hierarchical labels. Also, the method would rely on a numerical solver to “kick-start” the surrogate by solving L steps to get the initial, downsampled future states. These aren’t problems with the work, but should be mentioned for clarity.", "rating": 3, "confidence": 3, "soundness": null, "presentation": null, "contribution": null, "strengths": "", "weaknesses": "", "quality": 2, "clarity": 3, "significance": 2, "originality": 2, "strengths_and_weaknesses": "Strengths\n\n- Predicting multiple future states and taking in representations of a future state is interesting.\n\n- The baselines are done well to test different parts of the proposed method.\n\nWeaknesses\n\n- Perhaps there should be a baseline:  $ u_{n+1}',  u_{n+2}' = f_\\theta(u_n, u_{n+1}) $  , to see the effect of the downsampler on the implicit framework. It seems likely that downsampling may be necessary to prevent the network from learning the identity map between $u_{n+1}'$ and $u_{n+1}$.\n\n- There are also a number of stability-enhancing modifications that the authors cite: pushforward/unrolled training, spectral refinement/PDE-refiner. Unrolled training can improve Unet stability by quite a lot [1], and it would be worthwhile to see how implicit methods compare to spectral refinement or unrolled training.\n\n1. Bjoern List, Li-Wei Chen, Kartik Bali, Nils Thuerey, Differentiability in unrolled training of neural physics simulators on transient dynamics, https://www.sciencedirect.com/science/article/pii/S0045782524006960", "was_revised": true, "final_justification": "- Results are good, but limited to a single PDE and architecture\n- Extensions to other architectures are not obvious, and the improvements in other PDEs is not clear\n- The work has an interesting perspective and with more verification or improvements, it could be a promising method.", "rebuttal": "We appreciate your feedback!\n\n### W1\n> Perhaps there should be a baseline: $u_{n+1}', u_{n+2}' = f_{\\theta}(u_{n}, u_{n+1})$ , to see the effect of the downsampler on the implicit framework. It seems likely that downsampling may be necessary to prevent the network from learning the identity map between $u_{n+1}'$ and $u_{n+1}$\n\nExactly! Without downsampling or any form of information reduction, the model could exploit shortcuts by directly copying high-resolution latent variables. In the proposed setup, where $u_{n+1}', u_{n+2}' = f_{\\theta}(u_n, u_{n+1})$, the model could trivially map $u_{n+1}$ to $u_{n+1}'$ via identity, effectively reducing the function to $u_{n+2}' = f_{\\theta}(u_n, u_{n+1})$. This is equivalent to our 2-Step History baseline reported in Table 1.\nTo verify the empirical effects, we provide results for models without downsampling. We observe that both (1) $L=2$ (No Downsampling) and (2) $L=3$ (No Downsampling) consistently underperform compared to our $L=2$ and $L=3$ models.\n\n| Model                        | 1 Step                | 10th Step            | 25th Step            | 50th Step            | 75th Step            | 100th Step           |\n|------------------------------|-----------------------|----------------------|----------------------|----------------------|----------------------|----------------------|\n| L=1                          | 5.60e-04 (1.15e-04)  | 1.29e-02 (4.87e-03)  | 8.04e-02 (3.45e-02)  | 3.12e-01 (1.15e-01)  | 6.19e-01 (2.08e-01)  | 9.38e-01 (2.88e-01)  |\n| L=2 (No downsample)          | 9.90e-04 (2.14e-04)  | 1.29e-02 (7.02e-03)  | 5.87e-02 (2.89e-02)  | 2.10e-01 (8.89e-02)  | 4.47e-01 (1.76e-01)  | 7.74e-01 (2.52e-01)  |\n| L=3 (No downsample)          | 2.51e-03 (4.41e-04)  | 1.56e-02 (5.11e-03)  | 7.58e-02 (3.15e-02)  | 2.98e-01 (1.31e-01)  | 6.05e-01 (2.30e-01)  | 9.20e-01 (2.94e-01)  |\n| L=2 (ours)                   | **5.25e-04** (1.15e-04)  | 8.11e-03 (4.15e-03)  | 4.02e-02 (1.92e-02)  | 1.61e-01 (6.50e-02)  | 3.73e-01 (1.51e-01)  | 6.55e-01 (2.44e-01)  |\n| L=3 (ours)                   | 5.50e-04 (1.20e-04)  | **6.59e-03** (2.50e-03)  | **3.37e-02** (1.41e-02)  | **1.40e-01** (6.16e-02)  | **3.24e-01** (1.18e-01)  | **5.92e-01** (1.84e-01)  |\n\nTable: Comparison of short-term accuracy on downsampling.\n\n### W2\n\n> There are also a number of stability-enhancing modifications that the authors cite: pushforward/unrolled training, spectral refinement/PDE-refiner. Unrolled training can improve Unet stability by quite a lot [1], and it would be worthwhile to see how implicit methods compare to spectral refinement or unrolled training.\n\nThank you for your suggestions. We have run experiments with pushforward training.\nWe note that our initial attempt with pushforward training, following the codebase of [r2], did not converge during training—especially in chaotic regimes—as also observed in [r3,r4]. To address this, we introduced a warm-up phase of 50 epochs without pushforward, which significantly improved convergence and overall performance. Consistent with the discussion in [r3], we observe that while pushforward training improves long-term stability, it comes at the cost of short-term accuracy.\n\n| Model                        | 1 Step                | 10th Step            | 25th Step            | 50th Step            | 75th Step            | 100th Step           |\n|------------------------------|-----------------------|----------------------|----------------------|----------------------|----------------------|----------------------|\n| L=2    (ours)        | **5.25e-04** (1.15e-04)  | 8.11e-03 (4.15e-03)  | 4.02e-02 (1.92e-02)  | 1.61e-01 (6.50e-02)  | 3.73e-01 (1.51e-01)  | 6.55e-01 (2.44e-01)  |\n| L=3     (ours)        | 5.50e-04 (1.20e-04)  | **6.59e-03** (2.50e-03)  | **3.37e-02** (1.41e-02)  | **1.40e-01** (6.16e-02)  | **3.24e-01** (1.18e-01)  | **5.92e-01** (1.84e-01)  |\n| Pushforward    | 1.08e-03 (2.26e-04)  | 1.70e-02 (8.63e-03)  | 9.16e-02 (4.62e-02)  | 3.19e-01 (1.34e-01)  | 6.37e-01 (2.71e-01)  | 9.80e-01 (3.36e-01)  |\n\n| Model       | 1k steps               | 10k steps           | 100k steps            |\n|------------------------------|-----------------------|----------------------|----------------------| \n| Pushforward   | 98% | 96% | 88% | 73% |\n| L=3           | **100%** | **100%** | **100%** | **93%** |\n\nTable: Comparison of long-term stability rate with pushforward \n\nRegarding the unrolled training [r1], it introduces higher memory consumption and longer per-epoch training time due to additional forward passes and the requirement of differentiability across the rollouts (e.g., 2-step unrolling requires approximately 2 times larger computational time than our $L=3$ model). We will report these results in the next discussion phase and include the discussions in our refined version.\n\nAt last, we want to clarify that both unrolled training [r1] and pushforward training [r2] address the challenge of error accumulation in autoregressive rollouts, while our hierarchical implicit method focuses on exploiting the spatial-temporal structure and iteratively refining latent states using future frames' information.\nThese are orthogonal strategies, and combining them (e.g., unrolled training with our method) remains a promising direction, albeit with increased computational cost.\n\n### Q1\n\n>It would be interesting to see how the improvement of the implicit scheme varies with the dynamics of the PDE system. For simple systems, explicit schemes are often good enough, so would the performance gain be less?\n\nThank you for the suggestions. \nOur goal is to develop a robust framework for long-term simulation of complex and chaotic systems.. For reference, we have also conducted experiments on a relatively simpler 1D Kuramoto–Sivashinsky PDE.\n\n| L   | 1-step             | 5-step             | 10-step            | 15-step            | 20-step            | 25-step            | 50-step            | \n|-----|--------------------|--------------------|--------------------|--------------------|--------------------|--------------------|--------------------|\n| 1   | 5.43e-06 (1.67e-06) | 2.11e-04 (6.25e-05) | 4.16e-03 (4.30e-03) | 3.88e-02 (4.31e-02) | 1.44e-01 (1.02e-01) | **4.07e-01** (1.76e-01) | **2.03e+00** (1.68e-01) | \n| 3   | **5.21e-06** (1.51e-06) | **1.31e-04** (6.06e-05) | **3.78e-03** (5.09e-03) | **3.28e-02** (4.04e-02) | **1.40e-01** (9.70e-02) | **4.03e-01** (1.88e-01) | **2.00e+00** (1.88e-01) | \n\nTable: Comparison of mean squared error (MSE) on Kuramoto-Sivashinsky\n\n### Q2\n> In Figure 6, the GT for the jet system (top row) for time 25 and 50 look nearly identical (u_25, u_50). \n\nThank you for pointing this out! We accidentally uploaded the same figure for $u_{50}$. We will fix this in the final version. We will correct this in the final version. For additional visualizations, please refer to Appendix B, Figures 8–17, which cover various setups beyond Figure 6, including different resolutions, Reynolds numbers, short-term accuracy, and long-term stability.\n\n### Limitations\n> The authors don’t seem to present any limitations of their work. First, the method seems to be constrained to only Unets, since it doesn’t seem clear how to inject hierarchical information into other architectures... These aren’t problems with the work, but should be mentioned for clarity.\n\nWe would first like to clarify a key misunderstanding: our method does not require a numerical solver at inference time to \"kick-start\" the surrogate. Instead, the model is trained to estimate the full latent hierarchy from a single initial frame. Additional implementation details are provided in Appendix A.3.\n\nSpecifically, during the training, with a small probability, we let the model receive partially empty latent variables as input and predict the missing ones in the hierarchy latents.\nFor example in $L=2$, with a probability of 15%, we ask the model to take $[u_n, \\mathbf{0}]$ as input and output $z_{n+1}^{(1)}$. So that in the evaluation, with one extra forward run, we will have the full latent states $[u_n, \\hat{z}_{n+1}^{(1)}]$ to continue the autoregressive rollout.\nAll methods are evaluated given only one time frame initial condition. E.g., for generating an $T=1000$ sequence trajectory, our $L=3$ model will forward $1002$ times, where the extra $L-2=3$ forward passes are used to handle these corner cases.\nWe will update it in the refined version.\n\nRegarding the architectural scope, we do not consider the reliance on U-Net a limitation. U-Net is a highly flexible and widely adopted architecture across various domains and data modalities, and its hierarchical structure aligns naturally with our framework. That being said, our method is not inherently restricted to U-Nets and can potentially be extended to other architectures with hierarchical capabilities.\n\nFor mesh-based applications, while it is true that designing hierarchical representations requires additional considerations (e.g., defining consistent downsampling operators), prior work has proposed effective strategies such as vertex clustering, remeshing, and spectral downsampling. While adapting our framework to these settings would require additional architectural development, we view this as a promising direction with many exciting developments [r5] rather than a fundamental limitation.\n\n\n[r1] Differentiability in unrolled training of neural physics simulators on transient dynamics, Bjoern List, Li-Wei Chen, Kartik Bali, Nils Thuerey\n\n[r2] Message passing neural PDE solvers, Johannes Brandstetter, Daniel Worrall, and Max Welling\n\n[r3] Pdebench: An extensive benchmark for scientific machine learning, Makoto Takamoto, Timothy Praditia, Raphael Leiteritz, Dan MacKinlay, Francesco Alesiani, Dirk Pflüger, and Mathias Niepert.\n\n[r4] A dynamics-informed diffusion model for spatiotemporal forecasting, Salva Rühling Cachay, Bo Zhao, Hailey Joren, and Rose Yu\n\n[r5] Hao et al., High-Fidelity, Artist-Like 3D Mesh Generation at Scale."}]

Author rebuttal:

Meta-review: The paper introduces a contribution to the field of neural PDE solvers by addressing two challenges: long-term stability and error accumulation.  The proposed hierarchical implicit neural emulator is original and effective.
 
While some reviewers state that theoretical aspects could have been explored, they also found that the current empirical evidence is robust and convincing. Rebuttals have been key in the discussion and we are happy to propose to accept the paper.
Decision: Accept (poster)