Neural Emulator Superiority: When Machine Learning for PDEs Surpasses its Training Data
Authors: Felix Koehler
Nils Thuerey
Venue: NeurIPS 2025 poster ()
Keywords: PDEs
Neural Operators
Neural Emulators
Physics
Surrogate
Optimization
Primary area: machine_learning_for_sciences
TL;DR: We show that neural operators/neural emulators/neural simulators can outperform the numerical method that produced their training data.

Abstract: Neural operators or emulators for PDEs trained on data from numerical solvers are conventionally assumed to be limited by their training data's fidelity. We challenge this assumption by identifying "emulator superiority," where neural networks trained purely on low-fidelity solver data can achieve higher accuracy than those solvers when evaluated against a higher-fidelity reference. Our theoretical analysis reveals how the interplay between emulator inductive biases, training objectives, and numerical error characteristics enables superior performance during multi-step rollouts. We empirically validate this finding across different PDEs using standard neural architectures, demonstrating that emulators can implicitly learn dynamics that are more regularized or exhibit more favorable error accumulation properties than their training data, potentially surpassing training data limitations and mitigating numerical artifacts. This work prompts a re-evaluation of emulator benchmarking, suggesting neural emulators might achieve greater physical fidelity than their training source within specific operational regimes.

Project Page: https://tum-pbs.github.io/emulator-superiority

Neural Emulator Superiority: When Machine Learning for PDEs Surpasses its Training Data
Felix Koehler & Nils Thuerey
Technical University of Munich Munich Center for Machine Learning (MCML) {f wookkoeee thuerey/Stum.de de
Abstract
Neural operators or emulators for PDEs trained on data from numerical solvers are conventionally assumed be limited by their training data's fidelity. We challenge this assumption by identifying "emulator superiority." " where neural networks trained purely low-fidelity solver data achieve higher accuracy than those solvers when evaluated against a higher-fidelity reference. Our theoretical analysis reveals how the interplay between emulator inductive biases, training objectives, and numerical error characteristics enables superior performance during multi- step rollouts. We empirically validate this finding across different PDES using standard neural architectures, demonstrating that emulators can implicitly learn dynamies that are more regularized or exhibit more favorable error accumulation properties than their training data, potentially surpassing training data limitations and mitigating numerical artifacts. This work prompts a re-evaluation of emulator benchmarking, suggesting neural emulators might achieve greater physical fidelity than their training source within specific operational regimes.
Project Page: tum-pbs.github. wttthu..oo/muulttrr-suppp
Introduction
The numerical simulation of physical systems governed by partial differential equations (PDEs) involves an inherent tradeoff between computational cost and accuracy. Traditional numerical methods such as finite difference, finite element, and spectral approaches have well-understood error characteristics stemming from modeling, discretization, lincarization, and iterative solver truncation. These errors can be systematically reduced by refining the physical model, increasing resolution. employing higher-order schemes, or tightening solver tolerances but always at the cost of increased computational burden.
Machine learning approaches to PDE solving have emerged as promising alternatives, with neural net- works trained emulate the behavior of numerical solvers while offering computational advantages (Li et al., 2021a: la; Kochkov et al., 2024). The conventional wisdom suggests that neural emulators inherit the limitations of their training data: neural network trained data from low-fidelity solver should, at best, reproduce the accuracy of that solver.
In this paper, we demonstrate that this conventional wisdom does not always hold. We identify and analyze a phenomenon we call "emulator superiority." in which neural networks trained on data from low-fidelity numerical solvers can outperform those same solvers when tested against high-fidelity references. For this counterintuitive result, we provide both theoretical proof for linear PDEs (advection, diffusion, Poisson) and empirical evidence across linear advection and the nonlinear Burgers equation using diverse neural network architectures (convolutional networks, dilated ResNets, FNOs, UNets, Transformers). Crucially, this observation is different from previous work that learned improved coarse solvers from high-fidelity data (Bar-Sinai et al., 2019), mixed low-fidelity data with
39th Conference on Neural Information Processing Systems (NcurIPS 2025).

---

PDE 0ea = 7:0ga
ResNet, UNeel, FNO, Transformer, etc. (wit inductive bisc)
u
u iz
(if..
Emulator fo (leameds reguarized, more
&M P
coarse @ Fast, but Inascurate
Train Trajectories (contain structured numerical errors)
Superiority
physical model)
MiSoPP) <1
One-Step Supervised Traiing
MIP. P PIR.
M(P,P)
fine P Ground slow
Figure I: A PDE can be solved using different numerical algorithms. A coarse solver P produces data trajectories with inherent numerical errors (a). Neural emulators have inductive biases that enable them to learn regularized, more physical model fo (b). When we measure the accuracy of both emulator fo and coarse simulator P against higher fidelity ground truth P (c), the networks can exceed the training reference, giving rise to "emulator superiority" (d).
high-fidelity data (Lu et al., 2022) enhanced coarse data using physics-informed residuum soft penalizers (Li et al., 2021b). Instead, we purely train the neural emulator on low-fidelity data and still observe this superiority.
We do not suggest that emulator superiority is a universal guarantee; its occurrence and intensity are dependent on the specific PDE, the characteristics of the low-fidelity solver's errors, the chosen emulator architecture, the training objective, the evaluation metric, and the time horizon. However, demonstrating that such superiority exists, even within specific regimes, fundamentally challenges the notion that neural emulators are irrevocably bound by the accuracy limitations of their training data.
This finding suggests that neural emulators might implicitly learn regularized or corrected dynamics, potentially exceeding the fidelity of their low-fidelity training data. Furthermore, prompts re- evaluation of benchmarking practices, as emulators could potentially mitigate numerical artifacts even high-fideitty reference solvers. Since imperfections, e.g., due to insufficient resolution. numerical diffusion, or non-converged iterative schemes are often practically unavoidable or specifically intro- duced to reduce compute cost (Turek, 1999; Müller al., 2007), imperfections permeate into data trajectories used to train machine learning models. If the same data source is used for benchmarking, this creates the paradox that a better emulator might be unfairly penalized.
Our contributions are:
We formalize the differences between training objective and inference task of autoregressive neural emulators to identify the concept of superiority.
We provide theoretical proof for the existence of state-space superiority for the three major linear protypical PDEs using Fourier analysis. For this, we identify both forward and backward superiority based on how superiority occurs spectrally.
We confirm state-space and autoregressive superiority in emulator learning for the advection and Burgers equation using a wide range of nonlinear neural architectures.
Classical and Neural Solutions to PDEs
Numerical methods approximate solutions to continuous PDEs by transforming them into discrete problems solvable by computers. general PDE describing the evolution of field u(t, x) over time t and space æ E S C [R can be written as
L(u),
where L is spatial differential operator, which may be nonlinear and potentially inhomogencous.

---

2.1 Numerical PDE Solvers and Error Sources
Classical numerical solvers are based on frst-principled symbolic manipulations of the continuous description. Discretization in space (using methods like finite differences, finite elements, or spectral methods) and time (e.g,, Euler steps, Runge-Kure PatgeeKuttaaao methods) leads numerical schemes. Many such schemes can be expressed as discrete-time, first-order Markovian processes:
[++1|
(n(
h
where u. is the discrete state vector at time step it] on a grid characterized by h, and PA, is the discrete evolution operator encapsulating the numerical method and respecting boundary conditions.
The accuracy of Pa is limited by several error sources: (1) spatial discretization errors, (2) time discretization errors, (3) linearization errors when approximating nonlinear terms, and (4) iterative solver errors from terminating iterative methods prematurely. These errors accumulate interact over time, potentially leading to instability (often characterized by exponential error growth).
Different numerical schemes exhibit varying trade-off's regarding stability, accuracy, and compu- tational cost. For example, implicit methods often allow larger time steps (better stability) at the cost of increased numerical diffusion (a form of spatial error), while high-order methods reduce discretization errors but can be computationally expensive or are prone to oscillations Veque, 2007). Spectral methods restricted to simpler geometries (Boyd, 2001). The choice of domain, boundary conditions, and specific physical parameters further influences solver performance and error behavior.
offer high spatial accuracy for smooth solutions but are typically
2.2 Machine Learning Approaches for PDE Solving
The potential for computational speedups and tackling unsolved problems has motivated the use of neural networks (NNs) for helping solve PDEs. We are interested in approaches that learn operators that map function spaces, analogously to the numerical evolution operator Ps. in equation However, we relax the strict requirement of discretization invariance common in the neural operator literature (Kovachki et al., 2023; Gao et al., 2025) and instead focus on networks trained and applied at a fixed discretization h. We refer such a neural surrogate fe s P/, neural emulator (Koehler et al., 2024). Specifically, we consider autoregressive neural emulators that operate similarly to equation 2, predicting the next state based only on the current state: u, [++1] fo(u))
[t], While one can motivate non-Markovian approaches for states with truncated spectra (Ruiz et al., 2025), the Markovian assumption allows for a direct comparison with traditional numerical time steppers.
Supervised Training for Autoregressive Emulators common paradigm for training an autore- gressive emulator fo involves supervised learning on data generated by chosen numerical solver Pi. First, trajectories are generated by rolling out the solver from various initial conditions u, Ini
[0]
Tp. (n))
[0]
[0]
[S]
[0] Pr(
P. (M,
11,
Pairs of consecutive states {ug] un emulator Se is then trained to minimize one-step prediction error objective:
[ are sampled from these trajectories to form a dataset. The
L(O)= E (ulll Bool",.tt" " (1))), n/!0 ~I,, (s( (fo u "]
(3)
where is metric, e.g., a mean squared error (MSE). Variations exist, such as unrolled training, which minimizes errors over multiple steps (Brandstetter al., 2022; List et al., 2025), or incorporating physical constraints (Li et al., 2021b). However, the fundamental principle remains learning to mimic the behavior of Pr.
Generalization Challenges. Autoregressive emulators are evaluated by recursively rolling out the model for multiple time steps and comparing a reference trajectory,
M = E M.Li [ (f6(4))) ) (u.
(4)
Here, Z denotes applying the emulator times autoregressively. The test setup (denoted by tildes) might differ from the training setup in several ways, probing different aspects of generalization.
3

---

Temporal/Autoregressie generalization involves evaluating long-term accuracy and stability by testing for > i.c., for more autoregressive applications than during training. Exponential error growth el usually signals instability. State-space generalization examines performance on different initial conditions In A In or over extended time horizons beyond those used for training data generation, where the physics might enter a regime that is not well-represented by the training data. Metric generalization refers to evaluating the emulator using different metric Ç than the one used for training Ç. This allows for probing specific aspects of performance, such as accuracy in higher frequencies/modes Solver generalization evaluates against a different numerical simulator Pk. * Ph, where both solvers are consistent with the same underlying PDE equation 1 but employ different numerical implementations, varying in aspects such as scheme order or solver tolerance.
Emulators can be tested for a single type of generalization or a combination of many. This paper focuses on combining solver generalization with temporal or state-space generalization. We particularly interested in the scenario where the emulator fo is trained on data from a low-fidelity solver Ph. but evaluated against a high-fidelity reference solver Pn- This sets the stage for investigating whether Sa can outperform its own training data source Pn. when compared the more accurate reference Pr in different regimes of physics for multiple autoregressive evaluation steps.
2.3 Comparing Emulators and Emulator Superiority
Given two neural emulators f and in order to determine which one is better we could evaluate their rollout error e! and |00, respectively. Then, if g( = =eli /ell) < 1 we know that f is better than J at rollout step [t] (under given Č. In and PP.). Now, assume f was not another baseline neural emulator, but the training simulator PA used to inform fa. Then, we can define the superiority ratio as
[3 (f6(( 10) Pk 10
E 1 IIn
<0
(5)
10 wtn 5 ((% (u,
(0) PA (u (0)) )
Intuitively, one might expect el6 2 1; the emulator, trained to mimic Ph. should not more accurate than Pn, itself when compared better reference Pr. However, this intuition overlooks crucial differences between the training process and the evaluation scenario. The emulator fo optimizes the specific loss function of equation (often a one-step MSE) over a training distribution Ia. while superiority is measured via multi-step rollouts using potentially different initial conditions In. Furthermore, the function class of the emulator (e.g.. specific neural network architecture) imposes its own inductive biases. These differences create an opportunity for the emulator to find a parameterization 0 that, while faithfully reproducing Pa, in terms of the training objective, exhibits more favorable error accumulation properties during autoregressive rollouts or might generalize better to different regimes, leading to glt) < 1. We term this phenomenon emulator superiority.
This superiority is unlikely to hold universally (4[4, VI,, VỘ), but its existence, even within specific regimes (c.g., for some time horizon, certain initial conditions or eгror metrics), challenges the notion that emulators are strictly bound by their training data's fidelity. Based on whether superiority is due to better state-space generalization or autoregressive generalization, we define state-space superiority or autoregressive superiority, respectively. Note that while state-space superiority is already observable in the first step, i.e., e(L < autoregressive superiority requires at least two steps, i.., e(] and g(?) 1.
3 A Case Study for Superiority under Linear Conditions
In the following subsections, we demonstrate that state-space superiority can be achieved closed- form for three prototypical linear PDEs which are the basis for modelling many advanced phenomena, e.g., in fluid mechanics (Ferziger et al., 2020). Crucially, these derivations show how a simple linear emulator ansatz, by virtue of its own distinct error characteristics (stemming from its simple functional form, iee., inductive bias) relative to higher and lower modes in Fourier space, can outperform the more complex (but numerically imperfect) low- fwdfiidlity solver it was trained on; already after the first step, i.c., [L4) 1.

---

(a) Magnitude Error
(b) Phase Error
Amplifiying
Leading
Pelatiu freer
3.0
Dampening
Lagging
Scheme
Explict
Impict
-0.4
0.005
0.01
0.05 0.1
0.5
0.005
0.01
0.05 0.1
0.5
Relative Mode m/)
Relative Node - m/
Figure 2: Both explicit é, and implicit in schemes introduce an error in the modes' magnitude either in dampening (reducing) unstably amplifying (increasing) it. The former is visible as numerical diffusion in state space (a). This error more pronounced for higher modes and larger |r1]- For low Ø < 0.25, the explicit scheme has less numerical diffusion. For higher é, the implicit scheme is better. The explicit scheme incurs a stability restriction In < to be stable for all modes. Similarly, the explicit scheme is better in terms of phase enror within stability limits (b).
3.1 Advection Equation
The advection equation Ōzu + c0,u describes a simple linear and hyperbolic transport process, e.g., the movement of concentration of a species uff, via the homogeneous velocity c. For simplicity, let us consider the one-dimensional spatial domain f (0, 1) with periodic boundary conditions. We discretize the domain into N intervals of equal length Ar 1/N and denote by un E R
the vector nodal values of u at the left ends of each interval.
Explicit Finite Difference Schemes in State Space The class of finite difference schemes cquidistant grids can be represented by cross-correlations' with a kernel kn, i.c., they advance the state in time via i [++1] k, u where * denotes the cors-correlation with periodic boundary conditions (i.e., circular "same" padding). For the advection equation, the simplest possible scheme is a first-order upwind method, which is given by the kernel
K- 10 1-71 71
n1
71 0,
k, = [-
1 1+ 71
0]
71
< 0,
(6)
where we used 71 = - cN E R to absorb all parameters into one variable. Based on the sign of n, the scheme is either forward backward differencing. The absolute value of this variable, i.e. h1l. is also referred to as the Courant-Friedrichs--lewy (CFL. number (Le 2007).
Explicit Schemes via Fourier Multipliers Working under periodic boundary conditions has the merit that we can equally represent our state space a vector of complex-valuee Fourier coefficients ú, = (u,)) ) E N/2+1 where we used the real-valued Fourier transform Fn. In this case, the time-
[++1] Fi. (kn) ) üü where @ denotes the element-wise product. Without loss of generality, we restrict ourselves to 711 < 0. Then, the Fourier coefficients of the explicit first-order upwind method are given by Fz (k,,))))
stepping can equivalently be done in Fourier space via !
ê,. 1 + m "e ow where e o denotes the element-wise exponential applied to the vector scaled roots of unity
W, = - [0. 1, 1,,,, N/2 E CN2++1, which are based on the imaginary unit ;2 -1.
Implicit and Analytical Schemes in Fourier Space the Fourier domain, we can easily define two more numerical schemes: An implicit one that would require the solution to a linear system of equations in state space, and an analytical scheme. The implicit scheme is given by the Fourier multiplier i = 1 2 (I - 71 ow)) where Z denotes the element- division. The analytical scheme given by the Fourier multiplier å = ws More details on their derivation are presented in the appendix section A.1.I.
w..
Numerical Errors of the Schemes Since all numerical schemes, the explicit ê₁, the implicit ins and the analytical on diagonalize in Fourier space, we can analyze their mode-wise errors individually. The analytical time stepper serves as the reference since it does not introduce any
In deep learning context, this is often loosely referred convolution (LeCun al., 1989).
5

---

(a) Magnitude Su Superiority
(b) Phase Superiority
2.0
CFUys
3,0
-0.3
Superirity
-0.03
1.0
Relative
Made
y- MMV
0.0
0.005
0.01
0.05
0.1
0.5
0.005
0.01
0.05 0.2
0.5
Relative Mode d-- miN
Rolative Mode: mmn)
Figure 3: The superiority ratio for simple two-parameter ansatz when trained on data generated by an implicit first-order upwind fitted at relative mode 1b = M/N E [0, 1/2] and tested across relative mode Φ = m/N E [0, 1/2]. M is the one mode at which we fit the emulator, m is the one mode we test it at, and N is the spatial resolution. A value below 1 indicates that the learned method is superior to its training simulator. Superiority occurs in both magnitude and phase error "forwardly" in that we are superior for region with & 2.
numerical errors. To express the errors more conveniently, let us introduce the relative mode = m/ E [0, 1/2] Then, we can analyze the magnitude error ê - - and 7ūx 10 and the phase error (larg(
arg((x)) ))/ (rg ) and (larg(*)]] -argg ())) [a/rrg((x)) These can be expressed in closed-form and are functions of Φ and m We summarize an evaluation at representative values of & and 71 in figure 2. The key insight is that these errors are not random perturbations but
rather systematic deviations that follow well-defined patterns in Fourier space.
Training a Linear Emulator Let us return to the setting of training an autoregressive emulator fe for the mapping u [++1] fo (u). [t) Many architectures from mmage-to-image tasks in computer vision (Gupta and Brandstetter, 2023) or the more targeted field of neural operators (Kovachki et al., 2023) natural fit for this task. For simplicity, first consider simple linear cross- correlation/convllution. which equivalently diagonalizes in Fourier space
[++1] y
00
Oo
0]
[ * u,
[e+1] (6u + w
[ ū
If the train metric function is the mean squared error C(a, b) equation 3 diagonalizes. Let us further consider the scenario in which the initial condition distribution only produces states with solely mode M, i.e., un = csin..21) M - d) with some meaningful distributions over and In this case, the ansatz only informed at mode M, i.e. we get
-1l bll the training objective of
arg min (ΩM - fM) > 0o + 0 N - fM
with any of the aforementioned references fn. (ên, in, J. For a proof, see section 1.4.1. This complex-valued equation for the two real parameters 0, 0j R that can be solved in closed-form.
Superiority of a Trained Linear Emulator Let us use the abbreviation &i M/N [o, 1/2] for the relative mode at which we informed the ansatz. Then ar.. shall denote the linear ansatz with the parameters found at training mode & using any reference A, E {er, C,, ā, Again using the diagonalization properties and further assuming that the test distribution I, also only produces states with one mode m present (potentially differrnt from the training mode A2), we can express the superiority in magnitude simply gld = 1 1 1)
As such, the superiority ratio is a function of relative train mode 4b, relative test mode Φ, system properties and the choice of training reference ńy (that here also serves as the haseline) and testing reference Pn. Using the implicit scheme for training and the analytic scheme for testing, we find the one-step superiority {!l) depicted in figure 3. When training on an implicit first-order upwind method, can become better than this method for a different part of the Fourier space because we have an ansatz with a favorable generalization property.
Our theoretical analysis reveals several crucial conditions for achieving superiority. First, we must train a low-capacity linear ansatz on high-capacity (non-analytical) reference method. Second,

---

Diffusion
2.0
1.5
ls 10.00.00
q=1
Vo 0.3
4.0
1.0
Relarive Rishgg
Midi p
0.5
0.01
0.1
0.0
0.02
9.55
0.5
0.095
Figure 4: Superiority of an emulator fitted on relative mode v = M/N E [0, 1/2] for the diffusion (across different difficulties/ineensttess 72) and Poísson equation (across different numbers of iterations q of the Richardson iteration). Similar to the advection equation (figure 3), diffusion displays superiority forwardly, i.c., ó > v whereas Poisson has backward superiority, i.c., o < vb.
superiority emerges when testing at higher modes than those used during training, demonstrating a "forward" generalization property. Third, the difficulty parameter Y1 must be sufficiently large in absolute value, where numerical errors are most pronounced. discussion to how the superiority develops over multiple test rollout steps (see A. and we discuss that the results equivalently hold for physics-informed training (see section A.4.2).
the appendix, we extend the
3.2 Diffusion Equation
Similar reasoning can be applied to the linear diffusion equation, 0.u = vV again under periodic boundaries on the unit interval where a Fourier analysis diagonalizes standard finite difference schemes. Similar to the advection equation, we learn explicit functional form based data from an implicit scheme. In figure 4a we present the results of the closed-fome superiority expression. is again function of the fitting mode d M/N 1/2]. the testing mode m/N [0, 1/2]. and a new combined parameter "/2 = 2vAN > 0.
The diffusion equation demonstrates similar superiority pattern to that of advection. Specifically, superiority emerges in regions where Φ > v. indicating that emulators generalize effectively to higher modes than those used during training. The magnitude of this superiority effect increases with higher values of the difficulty parameter 7/2. Fundamentally, this reinforces our core finding: when emulator fo with favorable inductive biases (in this case, structure based simple cross- comelltionncoomollution mimicking a forward-in in ti-t central-in-spatt scheme) is trained on data from numerically imperfect solver Pn, (here, backward-in-tiee representations that outperform its training data when evaluated against a high-fidelity analytical solution Pa. i.e., is state-spac space superior, here already after the first time step.
central-ing central-in-space ipace scheme), it can develop
3.3 Poisson Equation
The superiority principle also be observed when solving stationary elliptic problems like the Poisson equation, -V = f. again under periodic boundary conditions. For this, we consider the low-fidelity scheme an unconverged iterative solver. For simplicity, we here use Richardson iteration. Such solvers act as smoothers: they efficiently damp high-frequency errors (modes 6 m/ near 1/2) but converge slowly for low-frequency errors (ó near 0). Let the resulting Fourier multiplier after q iterative steps (mapping source fm to approximate/rrunateed solution úm) be (q).
If we train an ansatz with a second-order accurate discretization of the second derivative in Fourier space, but with a free parameter 0, this trained ansatz can perform superior to the low-fidelity solver when evaluated against a fully analytic solution to the Poisson equation. This superiority can again be derived in closed-form and is also function of the fitting mode v = M/N E [0. 1/2] and the testing mode o m/N N E [0, 1/2). Here, it also depends on the number of iterations of the solver. We display configurations in figure 4b. Crucially, this setting we see backward superiority: the emulator outperforms the iterative baseline largely for á < v/. This contrasts with the forward superiority patterns seen earlier and lets us conclude that a state-spuce superiority is possible the direction where the reference produces more errors and the ansatz provides favorable inductive bias.
7

---

GondNet
Bilated Rosblet
FNG
Varilla Transformer
Autoregressive Gen erneliiization
Au Autoregressive Generalization Explicit Reference (6.. 0.9)
State Space Autoregr. Ceneralization Implicit Reference 3.0)
Implicit Reference 3,0)
2.0
(b)
Superont
1.5
30
Time
Time
Figure 5: Almost all considered emulator architectures achieve autoregressive superiority when the IC distribution does not change from training to test (a) & (b). Hence, superiority can also be achieved against the better explicit scheme within its stability limit (b). State-space superiority (under varied IC distribution) is only possible with the favorable local inductive bias of a ConvNet (c).
4 Extension to Nonlinear Settings
While our analysis in section 3 rigorously established the principle of state-space superiority for linear emulators under specific conditions, practical scientific machine learning often involves complex, nonlinear neural network emulators. We now shift our focus to these settings, investigating whether the core nasig, that an emulator's inductive biases and training process leads to dynamics that are implicitly more regularized or exhibit more stable error ieraagationextetln to nonlinear architectures and complex PDE systems. This transition also allows to explore the emergence of autoregressive superiority, phenomenon less readily observed with the linear ansatz of before.
4.1 Nonlinear Emulators for the Linear Advection Equation
Building our linear analysis of state-space superiority, we now investigate if nonlinear emula- tors Sprcispccifi ConvNets (Fukushima, 1980; LeCun et al., 1989), Dilated ResNets (He al., 2016; Yu and Koltun, 2016; Stachenfeld al. 2021), Fourier Neural Operators (FNOs) (Li et al., 2021a), and Transformers (Vaswani et al., 2017)-can also achieve superiority for the linear advection equation, potentially through different mechanisms, with key findings illustrated in figure 5.
Our experiments focus on two primary settings designed to probe different aspects of superiority. First, to assess autoregressive superiority, we train and test emulators using identical initial condition (IC) distributions. This is examined for emulators trained on data from an implicit scheme (6, at n -300, figure 5a) and an explicit first-order upwind scheme operating within its stability limit (7 -0.9, figure 5b). This setup tests if emulators can achieve more stable or accurate multi-step rollouts than their training solver. We observe clear autoregressive superiority. Nearly all tested architectures achieve superiority ratio g(() < 1 after some initial time steps. As expected, there cannot be superiority on the first time step (G(1) 2 1) because the emulators are trained to mimic the one-step behavior of the training solver these specific ICs. Remarkably, this superiority manifests even when training on data from the relatively accurate explicit scheme (figure 5b) for which we cannot prove superiority in the linear case (see A.1.5). This indicates that the nonlinear emulators learn dynamics that accumulate errors more favorably during rollouts compared to the solver they were trained on, phenomenon not readily achievable with simple linear emulators.
Second, to evaluate state-space superiority, we train emulators on data from the implicit scheme (6 Y1 3.00, figure 5c) and then test their generalization to ICs containing broader range of modes, thereby probing generalization to unseen physical states. We find that only the local feed-forward ConvNet able generalize the same extent the linear case. In contrast, architectures designed with global receptive fields, such as FNOs and Transformers, show limited ability to generalize in this partieular task of extrapolating from single-mode training multi-mode testing.
4.2 Nonlinear Emulators for the Burgers equation
We now extend our investigation to a scenario where both the PDE and the reference solvers are nonlinear, using the Burgers equation as a test case. The data generated by an implicit first-order upwind scheme. The key difference between our fow-ficclity training solver and the high-fidelity
8

---

Low-Fidelity (Picard-Truncater)
Neural Emulator (UNet)
High-Fidelity PPicar--Coveergel
UNet
1.00
(a)
(b)
(c)
0.75
1.2
(d)
Space
0.75
1.1
0.50
1.0
0.25
0.9
-0.75
0.00
0.8
Time Step
Time Step
Figure 6: When an autoregressive neural emulator is trained on the first two snapshots of a simulator that truncates the nonlinear solve of an implicit Burgers integrator (a), it learns a better representation of the shock propagation phase than the simulator (b), thereby more closely matching simulator which is nonlinearly converged (c). As such, becomes superior over its training simulator for certain number of time steps (d).
reference lies in the treatment of the nonlinearity: the low-fidelity solver employs Picard iteration that is truncated (i.e., not fully converged), while the high-fidelity reference uses a converged nonlinear solve (the schemes derived in appendix section B.2.1).
Not fully resolving the nonlinearity results in an incorrectly propagated shock (i.c., a wrong shock angle) depicted in figure 6a and figure 6c. The Burgers equation, courtesy of being nonlinear PDE, has the property that it develops a rich spectrum over the trajectory. The richer the spectrum, i.c., once are in the shock formation and shock propagation phase, the more errors nonlinearly truncated scheme makes (see figure 13 for more details). Hence, for this experiment, we train a vanilla UNet (Ronneberger et al., 2015) only the initial "smooth" frames from this low-fidelity solver. Despite training only on these smooth initial states, the UNet, when rolled out, learns to propagate shocks more accurately than its training data, developing shock angle that closely matches the high-fidelity reference (figure 6b,d) and thus demonstrating superiority. This superiority is hypothesized to stem from the UNet implicitly correcting the truncated nonlinear solve and its convolutional architecture generalizing effectively in Fourier space: training on smooth, low-frequency data enables it to better handle the emergent high-frequency shocks, mirroring mechanisms from section 3.
5 Limitations & Outlook
The phenomenon of emulator superiority, while intriguing, is primarily observed within specific operating regimes rather than a universal guarantee. This characteristic, however, aligns with the needs of many scientific and engineering applications, such as Model Predictive Control (MPC) (Rawlings al., 2017), where systems often operate within well-defined boundaries, over receding horizons, near particular trajectories. In such scenarios, an emulator's accuracy-cost proposition, even if confined to a specific regime, offers significant practical value. Future work should focus on precisely identifying these operational boundaries for superiority and understanding their dependence on the interplay between the PDE, solver characteristics, and emulator architecture.
Furthermore, the intensity of the superiority effect appears linked to the low-fidelity training solver being sufficiently inaccurate relative to the high-fidelity reference; highly precise training solver naturally reduces the margin for emulator improvement. While our examples often use solvers with clear error differentials to illustrate the principle, real-world scenarios frequently involve suboptimal or computationally constrained solvers (Turek, 1999; Müller et al., 2007). Thus, an important research avenue is to investigate the sensitivity of superiority to varying degrees of training data fidelity and to explore whether emulators can still offer advantages when trained moderately accurate data.
There is interesting relationship of emulator superiority and adversarial examples in image classification (Goodfellow et al., 2015). One could use equation and an optimization over it to find inputs u, (or corresponding distributions Σ,), metrics Ç. time horizones [t] that give rise to superiority. We believe this could be a valuable tool further understand this intriguing phenomenon.
[0]
Beyond these points, broader agenda includes deeper investigation into how different training strategies (e.g., unrolled training, physics-informed losses) and diverse architectural choices influ- ence superiority. Extending this analysis to more complex, multi multi-physics problems will also be

---

erucial. Ultimatlly, understanding and strategically hamesssing emulator superiority could lead to the development of more efficient and physically realistic emulation tools.
6
Related Work
The potential for computational speedups and tackling unsolved problems has motivated the use of neural networks (NNs) for helping solve PDEs. Physics-Informed Neural Networks (PINNs) approximate the continuous solution function x(t,x) directly (Lagaris et al., 1998; Raissi al., 2019). Machine learning models can also be integrated into traditional numerical methods to replace heuristic components, such turbulence models (Duraisamy al., 2019) or flux limiters (Discacciati et al., 2020). On a different abstraction level, NNs can be used to correct coarse solver outputs (Um al., 2020: Kochkov al.. 2021). A machine learning model can also be used to fully replace the numerical solver. for time-dependen problems typically in an autoregressive manner (Brandstetter al., 2022: Lippe al., 2023). The Neural Operator subfield poses stricter constraints on the architecture (Kovachki et al., 2023). A recent trend is to design emulators capable of solving multiple PDEs zero-shot or as foundation models (McCabe al., 2024; Herde et al., 2024; Holzschuh et al., 2025).
Oftentimes, the neural architectures contain inductive biases showing similarities to numerical solvers (Alt et al., 2023; McCabe et al., 2023), albeit them not being strictly bound to their numerical constraints (Koehler et al., 2024). This indicates why can use them compute-efficiet surrogates. While McGreivy and Hakim (2024) emphasize that some speedups might not hold to the strictest mathematical scrutiny, from a more practical perspective, there have been extensive positive results for neural approaches, c.g., in weather prediction (Kochkov et al., 2024).
Training these networks is oftentimes based on synthetic data produced by existing numerical schemes. Mohan et al. (2024) find that training on numerical simulation data containing structural errors can lead neural models learn these artifacts, potentially limiting their generalization despite appearing to learn. In the domain of particle physics, generative models have demonstrated "amplification, where high-fidelity distributions are recovered from limited or noisy training data (Butter al., 2021). To the best of our knowledge, our work is the first to identify how aspects of the emulator training pipeline can lead networks outperforming data with significant (deterministic) structural errors.
7 Conclusion & Impact
We introduced "emulator superiority." " where neural emulators trained on low-fidelity solver data outperform those solvers against high-fidelity reference in specific regimes. Our analysis reveals this stems from the interplay of emulator inductive biases, training objectives, solver error characteristics, and multi-step rollout dynamics, allowing emulators to implicitly correct structured errors. This work invites a reconsideration of the relationship between neural emulators and their training data, suggesting that the emulators' architectural inductive biases and learning dynamics can, under specific conditions, enable them transcend the accuracy limitations of their training simulators.
Our findings critically highlight the importance of high-quality testing references in emulator eval- uation, as flawed references with numerical artifacts obscure true performance, potentially penalizing emulators that achieve greater physical fidelity. This paradox challenges conventional benchmarking, urging careful interpretation of metrics considering the reference's 's own limitations. Consequently, developing robust evaluation frameworks that account for superiority effects is essen- tial for accurately assessing neural emulators in scientific computing. Embracing these insights can therefore foster the development of more powerful and reliable emulation tools, opening new avenues for scientific discovery and engineering innovation.
Acknowledgements
The authors thank Bjoern List, Patrick Schnell, Rene Winchenbach, Michael McCabe, and our anonymous reviewers for fruitful discussions. Felix Koehler acknowledges funding from the Munich Center for Machine Learning (MCML).
10

---

References
Tobias Alt, Karl Schrader, Matthias Augustin, Pascal Peter, and Joachim Weickert. Connections between numerical algorithms for pdes and neural networks. Journal of Mathematical Imaging and Vision, 65(1):185-200) 2023.
Yohai Bar-Sinai, Stephan Hoyer, Jason Hickey, and Michael P Brenner. Learning data-driven discretizations for partial differential equations. Proceedings of the National Academy of Sciences, 116631):15344-15549 2019.
John P Boyd. Chebyshev and Fourier spectral methods. Courier Corporation, 2001.
James Bradbury. Roy Frostig. Peter Hawkins, Matthew James Johnson, Chris Leary, Dougal Maclaurin, George Necula, Adam Paszke, Jake VanderPlas, Skye Wanderman- and Qiao Zhang. JAX: composable transformations of Python+NumPy programs, 2018. URL
http://gil htt://githun github. com/ homm/gooolejaaa
Johannes Brandstetter, Daniel E. Worrall, and Max Welling. Message passing neural PDE solvers. In The Tenth International Conference Learning Representations, ICLR 2022, Virtual Event, April 25-29, 2022. OpenReview.net, 2022. URL VS1x3HPYKSU.
net/ wwweeeeooormttttee id=
Anja Butter, Sascha Diefenbacher, Gregor Kasieczka, Benjamin Nachman, and Tilman Plehn. GAN- plifying event samples. SciPost Phys., 10(6))139, 2021, doi: 10.21468/SciPosstte
10.6.139.
Niccolo Discacciati, Jan S Hesthaven, and Deep Ray. Controlling oscillations in high-order discontinu- ous galerkin schemes using artificial viscosity tuned by neural networks. Journal of Computational Physics, 409: 109304, 2020.
Karthik Duraisamy, Gianluca laccarino, and Heng Xiao. Turbulence modeling in the age of data. Annual Review of Fluid Mechanics, 51:357-377, 2019.
Joel H Ferziger, Milovan Perić, and Robert Street. Computational methods for fluid dynamics, volume 4. Springer, 2020.
Kunihiko Fukushima. Neocogmitron: A self-organizing neural network model for mechanism of pattern recognition unaffected by shift in position. Biological cybernetics, 36(4):193-2)) 1980.
Wenhan Gao, Ruichen Xu, Yuefan Deng, and Yi Liu. Discretization-nnaariinn
the discretization mismatch errors in neural operators. In The Thirteenth International Conference on Learning Representations, 2025.
lan J. Goodfellow, Jonathon Shlens, and Christian Szegedy. Explaining and harnessing adversarial examples. In Yoshua Bengio and Yann LeCun, editors, 3rd International Conference on Learning Representations, ICLR 2015, San Diego, CA, USA, May 7-9, 2015, Conference Track Proceedings, 2015. URL http:://ax:iii org htrg//bb//ti 1112/1112..00 6572.
Jayesh K Gupta and Johannes Brandstetter. Towards multi-spatiotemporl -ccal scale generalized PDE modeling. Transactions on Machine Learning Research. 2023. ISSN 2835-8856. URL https: //openreview het//forum::dddPSDDDrr
Kaiming Hc. Xiangyu Zhang, Shaoqing Ren, and Jian Sun. Deep residual learning for image recognition. In 2016 IEEE Conference Computer Vision and Pattern Recognition, CVPR 2016, Las Vegas, NV, USA, June 27-30, 2016, pages 770-778. IEEE Computer Society, 2016. doi: 10.1 10.1100//VP..2006.990 URL https: //woi./////oo org/10 10. 1119/CCVR..0 www.C.2........... 2016.90.
Maximilian Herde, Bogdan Raonic, Tobias Rohner, Roger Käppeli, Roberto Molinaro, Emmanuel de Bézenac, and Siddhartha Mishra. Poseidon: Efficient foundation models for pdes. Advances in Neural Information Processing Systems, 37:72525-72624, 2024.
Benjamin Holzschuh, Qiang Liu, Georg Kohl, and Nils Thuerey. PDE-transformer: Efficient and versatile transformers for physics simulations. In Aarti Singh. Maryam Fazel, Daniel Hsu, Si- mon Lacoste-Julien, Felix Berkenkamp. Tegan Maharaj, Kiri Wagstaff, and Jerry Zhu, editors, Proceedings of the 42nd International Conference on Machine Learning, volume 267 of Pro- ceedings of Machine Learning Research, pages 23562-23602. PMLR, 13-19 Jul 2025. URL. https: /ttp://eerccee...e///e mlr, www.ppressvvvee
htmll.

---

Diederik P. Kingma and Jimmy Ba. Adam: A method for stochastic optimization. In Yoshua Bengio and Yann LeCun, editors, 3rd International Conference on Learning Representations, ICLR 2015, San Diego, CA, USA, May 7-9, 2015, Conference Track Proceedings, 2015. URL http:///axivv org/abs///.. //////////// ftbs/1412.690 6980.
Dmitrii Kochkov, Jamie A Smith, Ayya Alieva, Qing Wang, Michael P Brenner, and Stephan Hoyer. Machine learning-acces accelerated computational fluid dynamies. Proceedings of the National Academy of Sciences, 118(21), 2021.
Dmitrii Kochkov, Janni Yuval, Ian Langmore, Peter Norgaard, Jamie Smith, Griffin Mooers, Milan Klöwer, James Lottes, Stephan Rasp, Peter Düben, et al. Neural general circulation models for weather and climate. Nature, 632(8027):1000-1000 2024.
Felix Koehler, Simon Niedermayr, Rüdiger Westermann, and Nils Thuerey. Apebench: A benchmark for autoregressive neural emulators of pdes. In Amir Globersons, Lester Mackey, Danielle Belgrave, Angela Fan, UIrich Paquet, Jakub M. Tomczak, and Cheng Zhang, editors, Advances in Neural Information Processing Systems 38: Annual Conference on Neural Information Processing Systems 2024, NeurIPS 2024, Vancouver, BC, Canada, December 10 15, 2024, 2024. URL http:// nips. cw/paper_fiiiiecccceeeeeeeeeeepe wwl.fffleeapppee htppp/pppp///////// 2atr///222//aa/////// html.
hww.oodbbbbeeee--baaaaaatcc..a.e Abstract- Daattt--taatttttte and _aaddd_Beemaaaeee.e Track.
Nikola B. Kovachki, Zongyi Li, Burigede Liu, Kamyar Azizzadeneshelii Kaushik Bhattacharya, Andrew M. Stuart, and Anima Anandkumar. Neural operator: Learning maps between function spaces with applications to pdes. J. Mach. Learn. Res., 24:89:1-89:97, 2023. URL http:
html.
Isaac E. Lagaris, Aristidis Likas, and Dimitrios Fotiadis. Artificial neural networks for solving ordinary and partial differential equations. IEEE Trans. Neural Nenworks, 965)9877-1000, 1998. doi: 10.1109/72.712177. URL https: //te::////o... httt::://doi.org/111/
Yann LeCun, Bernhard Boser, John S Denker, Donnie Henderson, Richard E Howard, Wayne Hubbard, and Lawrence D Jackel. Backpropagation applied to handwritten zip code recognition. Neural computation, 1(4):541-551) 1989.
Randall Le Finite difference methods for orlinary and partial differential equations: steady- state and time-dependent problems. SIAM, 2007.
Zongyi Li, Nikola Borislavov Kovachki, Kamyar Azizzadenesheli, Burigede Liu, Kaushik Bhat- tacharya, Andrew M. Stuart, and Anima Anandkumar. Fourier neural operator for parametric partial differential equations. 2021a. URL https: /tttp::///penrrvv// mntt/foruurrd--ccPPXytttmmmm
Zongyi Li, Hongkai Zheng, Nikola Kovachki, David Jin, Haoxuan Chen, Burigede Liu, Kamyar Azizzadenesheli, and Anima Anandkumar. innsore-informed neural operator for learning partial differential equations. arXiv preprint arXiv:2111.03799 2021b.
Phillip Lippe. Bas Veeling, Paris Perdikaris, Richard E. Turner, and Johannes Brandstetter. Pde-refiner: Achieving accurate long rollouts with neural PDE solvers. 2023. URL http://papers... www..iccppppp..eeeeeeeee caaperffile//pppeeeee Cawr:ao----sttttt-----oaoaa
.ttml.
Bjoern List, Li-Wei Chen, Kartik Bali, and Nils Thuerey. Differentiability in unrolled training of neural physics simulators on transient dynamics. Computer Methods in Applied Mechanics and Engineering, 433:1 17441, 2025.
Lu Lu, Raphaël Pestourie, Steven Johnson, and Giuseppe Romano. Multifidelity deep neural operators for efficient learning of partial differential equations with application to fast inverse design of nanoscale heat transport. Physical Review Research, 4(2):023210, 2022.
Michael McCabe, Peter Harrington, Shashank Subramanian, and Jed Brown. Towards stability of autoregressive neural operators. Transactions on Machine Learning Research, 2023. ISSN 2835-8856. URL https: http://prrreeeeeeeeeeeeoorrrerea
12

---

Michael McCabe, Bruno Régaldo-Sain Blancard, Liam Holden Parker, Ruben Ohana, Miles D. Cranmer, Alberto Bietti, Michael Eickenberg, Siavash Golkar, Géraud Krawezik, François Lanusse, Mariel Pettee, Tiberiu Tesileanu, Kyunghyun Cho, and Shirley Ho. Multiple physics pretraining for spatiotemporal surrogate models. Amir Globersons, Lester Mackey, Danielle Belgrave, Angela Fan, Ulrich Paquet, Jakub M. Tomczak, and Cheng Zhang, ed- itors, Advances in Neural Information Processing Systems 38: Annual Conference on Neu- ral Information Processing Systems 2024, NeurIPS 2024, Vancouver, BC, Canada, Decem- ber 10 - 15, 2024, 2024, URL /ttt::://pppprse nips. cc/ cc/paper_aape httpppaapppleee/a/llaaaaaaaaeeee d7 2tb.aa5aade22771114000000eeeeee0 wtt.aaeessssccoo--aaaatt------- Censsracc-Conferenee www.fffllllpppeee hallee/paperr2024//aaat html..
Nick McGreivy and Ammar Hakim. Weak baselines and reporting biases lead to overoptimism in machine learning for fluid-related partial differential equations. Nat. Mac. Intell, 6(10): 1256-1269, 2024. doi: 01108//44255-0024-0887777 URL https://dol. http:///...rrg//.... ww..ooo./..... 1038/ s42256-024
Aaron Meurer, Christopher P. Smith, Mateusz Paprocki, Ondfej Čertík, Sergey B. Kirpichev, Matthew Rocklin, Amit Kumar, Sergiu Ivanov, Jason K. Moore, Sartaj Singh, Thilina Rathnayake, Sean Vig, Brian E. Granger, Richard P. Muller, Francesco Bonazzi, Harsh Gupta, Shivam Vats, Fredrik Johansson, Fabian Pedregosa, Matthew J. Curry, Andy R. Terrel, Štěpán Roučka, Ashutosh Saboo, Isuru Fernando, Sumith Kulal, Robert Cimrman, and Anthony Scopatz. Sympy: symbolic computing in python. PeerJ Computer Science, 3:3103, January 2017. ISSN 2376-5992. doi: 10.7717/peerj-cs.100 103. URL https://dol org/ 10. 7717/71117/// c 1033
Arvind Mohan, Ashesh Chattopadhyay, and Jonah Miller. What you see is not what you get: Neural partial differential equations and the illusion of learning. arXiv preprint arXiv:2411.1510/ 2024.
Matthias Müller, Bruno Heidelberger, Marcus Hennix. and John Ratcliff. Position based dynamics. Journal of Visual Communication and Image Representation, 18(2):109-118 2007.
Maziar Raissi, Paris Perdikaris, and George E. Karniadakis. Physics-informed informed neural networks: A deep learning framework for solving forward and inverse problems involving nonlinear partial differential equations. J. Comput. Phys., 378:686-707, 2019. doi: 10.1016/..CP.201... 10.004.
URL httph:tt://// http::://ddoo.ooo org//..0....1.... 1016/j j- jep 2018 10.045.
James Blake Rawlings, David Q Mayne, and Moritz Diehl. Model predictive control: theory, computation, and design, volume 2. Nob Hill Publishing Madison, WT, 2017.
Olaf Ronneberger, Philipp Fischer, and Thomas Brox. U-net: Convolutional networks for biomedical image segmentation. In International Conference on Medical image computing and computer- assisted intervention, pages 234-241. Springer, 2015.
Ricardo Buitrago Ruiz, Tanya Marwah, Albert Gu, and Andrej Risteski. On the benefits of memory for modeling dependent pdes. In The Thirteenth International Conference Learning Representations, ICLR 2025, Singapore, April 24-28, 2025. OpenReview.net, 2025. URL https:
//openreview.net// http:iiiwwwnnneftuuuteeea 1d==okga5k333..
Kimberly Stachenfeld, Drummond Fielding, Dmitrii Kochkov, Miles D. Cranmer, Tobias Pfaff, Jonathan Godwin, Can Cui, Shirley Ho, Peter W. Battaglia, and Alvaro Sanchez- Learned coarse models for efficient turbulence simulation. CORR, abs/2112. 2021. URL. https: / /ttp///axvvv///////////e/////avt httgg/ab//21111 15275.
Stefan Turek. Efficient Solvers for Incompressible Flow Problems: An Algorithmic and Computational Approach, volume 6. Springer Science & Business Media, 1999.
Kiwon Um, Robert Brand, Yun (Raymond) Fei, Philipp Holl. and Nils Thuerey. Solver-i--thellop:: Learning from differentiable physics to interact with iterative pde-solvers.
2020. URL. https: htttps://proceedingg.. wwwpeeedddinnnn...nssss.. wttpppppscccppppee/200//// Abst html.
Ashish Vaswani, Noam Shazeer, Niki Parmar, Jakob Uszkoreit, Llion Jones, Aidan N Gomez, Eukasz Kaiser, and Illia Polosukhin. Attention all you need. Advances in neural information processing systems, 30, 2017.
13

---

Fisher Yu and Vladlen Koltun. sclti context aggregation by dilated convolutions. In Yoshua Bengio and Yann LeCun, editors, 4th International Conference on Learning Representations, ICLR 2016, San Juan, Puerio Rico, May 2-4, 2016, Conference Track Proceedings, 2016. URL http:// /ttxp://axxiv.../
waxxi..rg/aaiii. org.abs/ htt/a//ss/1....11 07122.
14

---

Checklist
1. Claims
Question:
Do the main claims made in the abstract and introduction accurately reflect the
paper's contributions and scope?
Answer: [Yes]
Justification: We provide the definition of emulator superiority in section 2.3 and present theoretical proof in section 3 which extend experimental evidence in section 4.
Guidelines:
The answer NA means that the abstract and introduction do not include the claims made in 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 is clear that these goals attained by the paper.
2. Limitations
Question:
Does the paper discuss the limitations the work performed by the authors?
Answer: [Yes]
Justification: We discuss the limitations of our work in section 5 and also highlight in the introduction that emulator superiority not universal guarantee but limited to specific operating regimes.
Guidelines:
The answer NA means that the paper has 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 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 the scope of the claims made, e.g., if the approach was only tested on a few datasets or with few runs. In general, empirical results often depend implicit assumptions, which should be articulated.
The authors should reflect the factors that influence the performance of the approach. For example, facial recognition algorithm may perform poorly when image resolution 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 to 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 used by reviewers as grounds for rejection, worse outcome might be that reviewers discover limitations that aren't 't acknowledged in the paper. The authors should use their best judgment and recognize that individual actions in favor of transparency play 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 prool's
Question: For each theoretical result, does the paper provide the full of assumptions and complete (and correct) proof?
15

---

Answer: [Yes]
Justification: Our theoretical results are described on a higher level in the main part of section and we provide more details the derivations section A.
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 or the supplemental material, but if they appear in the supplemental material, the authors are encouraged to provide a short proof sketch to provide intuition.
Inversely, any informal proof provided in the core of the paper should be complemented by formal proofs provided in appendix or supplemental material.
Theorems and Lemmas that the proof relies upon should be 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 affects the main claims and/or conclusions of the paper (regardless whether the code and data are provided or not)?
Answer: [Yes]
Our theoretical results self-contained with the extended description in sec- tion A. Moreover, we provide Jupyter notebooks in the supplemental material using SymPy (Meurer et al., 2017). For the experiments section 4, we provide all hyperparameters and further settings in section B. All code used to generate the results is provided in the supplemental material. will be made available upon acceptance.
Justification:
Guidelines:
The answer NA means that the paper does not include experiments. If the paper includes experiments, 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.
the contribution is dataset and/or model, the authors should describe the steps taken to make their results reproducible or verifiable.
Depending 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 be necessary to either make it possible for others to replicate the model with the same dataset, provide access 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 replicate the results, access 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.
specific model and empirical evaluation, may
While NeurlPS does not require releasing code, the conference does require submis- sions to provide some reasonable avenue for reproducibility, which may depend on the nature of the contribution. For example
(a) If the contribution is primarily a new algorithm, the paper should make it clear how reproduce that algorithm.
(b) If the contribution is primarily
model architecture, the paper should describe the architecture clearly and fully.
(c) the contribution is a new model (c.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 or instructions for how to construct the dataset).
(d) We recognize that reproducibility may be tricky in some cases, in which case authors are welcome describe the partieular way they provide for reproducibilitty In the case of closed-source models, it may be that access to the model is limited in
16

---

some way (e.g., to registered users), but it should be possible for other researchers have some path 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 faithfully reproduce the main experimental results, as described in supplemental material?
Answer: [Yes]
Justification: All code used to generate the results, including Jupyter notebooks with SymPy manipulations for theoretical results, is provided in the supplemental material. We will release this material 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/
caassinnooiccyy for more details.
While we encourage the release of code and data, we understand that this might not be possible, so "No" is an acceptable answer. Papers cannot be rejected simply for not including code, unless this is central to the contribution (e.g., for new open-source benchmark).
The instructions should contain the exact command and environment needed to run to reproduce the results. See the NeurIPS code and data submission guidelines (https: /Cooeeueoeeeeeeeeeee
/ /nipp. wwtppppllcggccieeeoocoooooco Stte:////gguieeeeeeuuuuuiiiiii sesiissssiooiooiieeeeeeoo onPolicy) for more details.
The authors should provide instructions on data access and preparation, including how 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).
Providing much information possible in supplemental material (appended to the paper) recommended, but including URLs to data and code permitted.
6. Experimental setting/details
Question: 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?
Answer: [Yes]
Justification:
All hyperparameters and further settings are provided in section B.
Guidelines:
The answer NA means that the paper does include experiments.
The experimental setting should be presented in 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, as supplemental material.
7. Experiment statistical significance
Does the paper report error bars suitably and correctly defined or other appropriate information about the statistical significance of the experiments?
Question:
Answer: [Yes]
Justification: For the theoretical results of section 3, no errorbars are needed because the arising convex optimization problem can be solved in closed-form and we test on a single mode. For the neural network experiments in section 4, provide the error bars for all runs in terms of the median out of at least 10 different re-initialization seeds together with the 50% percentage interval. More details on this can be found in section B.
Guidelines:
17

---

The answer NA means that the paper does not include experiments.
The authors should answer "Yes" if the results are accompanied by enror 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 capturing should be clearly stated (for example, train/test split, initialization, random drawing of some parameter, or overall
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 the standard error of the mean.
is OK to report I-sigma error bars, but one should state it. 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 that would yield results that are out range (e.g. negative error rates).
If error bars reported in tables 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 the com- puter resources (type of compute workers, memory, time of execution) needed to reproduce the experiments?
Question:
Answer: [Yes]
Justification: Our experiments are rather small in size, but we still provide the information the computer resources used in section B.
Guidelines:
The answer NA means that the paper does not include experiments.
The paper should indicate the type of compute workers CPU or GPU, internal cluster, or cloud provider, including relevant memory and storage.
The paper should provide the amount compute required for each of the individual experimental runs 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 it 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: / httpppp.cccppbuuee httprppps.//pppppi//hccccccoaaa wwwhhhhhhiiieeeeeeeee.
Answer: [Yes]
Justification: We have reviewed the NeurlPS Code of Ethics and confirm that our work conforms
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 deviation from the Code of Ethics.
The authors should make sure to preserve anonymity (e.g., there a special consid- eration due to laws 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: [NA]
18

---

Justification: Our work is rather foundational and does not have direct societal impacts. We work on the theoretical understanding of emulators and their use in scientific computing.
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 (c.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 is being used as intended but gives incorrect results, and harms following from (intentional unintentional) misuse of the technology.
If there are negative societal impacts, the authors could also discuss possible mitigation strategies (c.g., gated release of models, providing defenses in addition to attacks, mechanisms for monitoring misuse, mechanisms to monitor how system learns from feedback over time, improving the efficiency and accessibility of ML).
11. Safeguards
Does the paper describe safeguards that have been put in place for responsible release of data models that have high risk for misuse (e.g., pretrained language models, image generators, or scraped datasets)?
Question:
Answer: [NA]
Justification: Our work does not use data or models that have a high risk for misuse.
Guidelines:
The answer NA means that the paper poses no such risks.
Released models that have high risk for misuse or dual-use usal-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 or 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: [Yes]
Justification: When using others work (PDE data generation and neural architecture) we cite the original papers.
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.
19

---

The authors should state which version of the asset is 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 re-packaged, both the original license and the license of the derived asset (if it has changed) should provided.
If this information is not available online, the authors are encouraged to reach out to the asset'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:
There are no new assets introduced in the paper.
Guidelines:
The answer NA means that the paper does not release new assets.
Researchers should communicate the details the dataset/code/mooel as part of their submissions via structured templates. This includes details about training, license, limitations, etc.
The paper should discuss whether and how consent was obtained from people whose asset used.
At submission time, remember to anonymize your assets (if applicable). You can either create 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 as details about compensation (if any)?
Answer: [NA]
Justification: There are crowdsourcing experiments research with human subjects in the paper.
Guidelines:
The answer NA means that the paper does not involve crowdsoureing nor research with human subjects.
Including this information in 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 in the main paper.
According the NeurIPS Code of Ethics, workers involved in data collection, curation, or other labor should be paid at least the minimum wage in the country of the data collector.
15. Institutional review board (IRB) approvals equivalent for research with human subjects
Question: Does the paper describe potential risks incurred by study participants, whether such risks were disclosed to 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: There are no crowdsourcing experiments or research with human subjects in the paper.
20

---

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 important, original, or non-standard component of the core methods this research? Note that the LLM used only for writing, editing, or formatting purposes and does not impact the core methodology. scientific rigorousness, or originality of the research, declaration is not required.
Answer: [NA]
Justification: There no LLMs used in the paper.
Guidelines:
The answer NA means that the core method development in this research does not involve LLMs as any important, original, non-standard components.
Please refer to our L.LM policy http:://eearrreeec////ofeeeeeeee /ct////ofeerreccc//2/////LL for what should or should not be described.
21

---

A More Details on the theoretical considerations
The derivations presented herein are supported by symbolic computations performed with SymPy (Meurer et al., 2017). To facilitate reproducibility, Jupyter notebooks are available in our supple- mental material: github,
wtttttumm////mmaooorrroor superiorty.... This section aims to provide a more comprehensive exposition of the mathematical foundations underlying the numerical schemes and the concept of emulator superiority discussed in the main text.
A.1 Advection Equation
We begin by considering the one-dimensional linear advection equation on a periodic domain of unit length (L 1):
8,u + cd2u = 0, wd0 t) - (1,,
where E R. is the constant advection speed. To simplify the analysis by reducing the number of free parameters, we adopt the following reparametrized form
71 O2u, NA
0,u
(10)
where
(11)
a combined variable to absorb experimental parameters. Its absolute, i.e., 171/, also referred to as the Courant-Friedrich--Leew (CFL) number. The number of discretization points N and the time step size At will cancel out in the following analysis. As discussed in the main text section 3.1, we choose an equidistant discretization of the spatial domain S2 into N intervals of equal length Az 1/N. Then, we consider the left endpoint of each interval a nodal degree freedom. Collectively, we can write the vector of grid coordinates as X [0. Ar. 2Az,.... (N - 1)Ar)" E RN, The left boundary point is part of our discretization, but the right boundary point not. Let u, E RR- be the vector of nodal values of at the coordinates X;. An upper index in square brackets is used to denote time index, i.e., u.
[0] would be the initial state, // is the state time step later, etc.
A.1.1 Extended Derivations for the Advection Equation Schemes
The simplest possible discretization in space is based on a first-order approximation to the first spatial derivative. To ensure stability of the scheme, have choose between forward and backward differencing based on the direction of the advection speed c, i.e., based on the sign of 71- The stable configuration would be forward difference if negative and backward difference if c is positive. Hence, since "/, reverses the sign, as defined in equation 11, we have to use forward difference if 7/ > 0 and backward difference if 71 0. As such, we get the following two possible approximations to the advection equation
d
71
Au (ui+1 wz)
if 71 > if i < 0
dt " NAt
(12)
Ar (u; u-11)
Here we used u, to index the nodal degrees of freedom at spatial position Ię = iAz Since Ar L/N and L 1, that 1/Ar and 1/N cancel, leaving
d di u4
71 [u411 U4 if m u 4:-1
(13)
Al
if m <
Further, we can Fuler approximation discretize in time. If we evaluate the right-hand side the previous point in time, we obtain an explicit first-order upwind method.
[++]]
[4]
71
u (
if >
(14)
At
if
<
where we can easily re-arrange for the next state in time and see that At cancels giving
u [t]
[++]]
u u
if >
+
m
1
if 21 <
(15)
22

---

We see that the discrete solution is solely parameterized by combined variable M. Without loss of generality, let us assume that ?1 O. Then, we can write the update on the entire state vector as
[++1]
71 0
71
8
+
0
0
0
U
u
U2 [+1
71
+
71
(16)
:
0
71
1+ 7
0
[++1]
[e]
N-11
0
0
71
1+
71
21
=A
=u
Clearly, the matrix A is a circulant matrix. Hence, we could equivalently express the update rule via a cross-correlation on periodic boundaries as
u |++ I-n
1 71
*oo
u
(17)
where * denotes the periodic cross-correlation. It is known that circulant matrices diagonalize using the discrete Fourier transform (DFT) matrix
1
1
w
w
w
F
I
w
w
2
(18)
N-1
(NN-1)
w (N-1)((-1))
where w = e -12=NN is the N-th primitive root of unity. For our explicit upwind scheme, this gives us the Fourier multiplier
ē = 11+1 - m own
(19)
where e ow denotes the element-win wise exponential applied to the vector of scaled roots of unity
Wa 1 1. N/2 11" E CN/2+1 These allow to advance the state in Fourier space via
[++1
é Púúl
(20)
where O denotes the element- wise product.
Implicit First-Order Upwind Scheme If instead of evaluating the right-hand side at the previous point in time, we evaluate it at the next point in time, obtain implicit scheme. Assuming again T1 < 0, in state space it would read
I 1 0)
[++1]
!+]
(21)
Hence, in order to advance in time in state space, we would have to solve linear system. Gladly, the diagonalization property of the DFT matrix allows us solve the system by mode-wise division in Fourier space. This gives the Fourier multiplier
àn Ø (1 - n + www)
(22)
where 2 denotes the element-wise division.
Analytical Scheme in Fourier Space The explicit and implicit first-order upwind schemes incur an error due to the approximation of space and time. Both can be avoided by using exact solution in Fourier space which is possible for all linear PDEs under periodic boundaries with bandlimited initial conditions. First, use the exact first derivative in Fourier space given
2T
Fr(0.u)) =
0
0
1
... N/2 + 1 (um)
(23)
L
where dj E N/2+1 are the mode-wise complex-valuee Fourier multipliers with which the state in Fourier space ún is element-wiss multiplied find the Fourier representation of the first derivative. Hence, Fourier space, the semi-discrete advection equation reads
d
71 d C ú
(24)
dt
N
23

---

Working in Fourier space lets us diagonalize the system. This makes the ODE much easier to solve. We can use the matrix exponential, which is simple to calculate when matrices are diagonal. This gives us an analytical time stepper in Fourier space.
[++1]
ods NAC
h
(25)
6
where e" denotes the blement-wise exponential. Let us further reformulate the Fourier multiplier á as
d expo(d), (dr,
71
At)
(26)
NAt
2T
exp=
0 1
2
N/2 + m)
(27)
N
(28)
where we used L = 1.
A.1.2 Error of the Schemes
The diagonalization in Fourier space allows express the mode-wist error done by each scheme in terms of the error of the Fourier multiplier. For this, we use the following representations of the three schemes (again assuming 71 0) for a relative mode Φ = m/N N [0, 1/2]
è.p = 1 + m
716" 12
(29)
1
06 ās = E i'n 2nd
(30)
+ 7-e-112200 12m6
(31)
Since all multipliers are complex-valuede we can compare the errors against the analytical scheme both in terms of magnitude and phase. The magnitude of the analytical reference is always lass - 1 which is to be expected since the pure hyperbolic transport should not alter the magnitude of the solution, i.e., there should be artificial diffusion.
Then, we get first the magnitude error of the explicit scheme
|2 |åsl
as cos ( atan2 (0, 4 .7 sin 2 (N Φ) + 4.71-1 sin? (t o) +))
(32)
(4 7p sin?
(π Φ) + 4 "1 sin? (T Φ) + 1| 1,
and the magnitude error of the implicit scheme
[2 àø
+
(33)
O
4 7/ sin? (IT Φ) - m sin? (I )) +
Moreover, have the phase error of the explicit scheme
-
arg (n 71 e 2 + 1)
(34)
(arg(ā:))
arg (62- +0)
and the phase error of the implicit scheme
2ir (
|arg( -arg(a)) |arg(ás)
arg
1.
(35)
arg (c2
Those expressions have been plotted in figure 2 in the main text for representative choices of Φ and
Ti.
24

---

A.1.3 Closed-Form Emulator trained on the Implicit Scheme
As discussed in the main text, we will inform an emulator fo at one relative mode & M/N. The functional form of the ansatz will be similar to the explicit scheme and given by
[++1] 1] 2/
[0]
80
Koo
Wi [d]
(36)
in the state space, i.e., has the same upwinding bias the explicit scheme and hence should only be used for the assumed "/1 < O. We find its Fourier multiplier as
q = 0o + 0 ow..
(37)
When fitting against the implicit scheme at relative mode iyp we find the following solutions for the complex-valued equation
8u
(38)
4- sin? (át v) 4 71 sin? (TT +) +1'
-4 71 sin? (T v)) 71 +1
01 4. 71 sin (π 2) - 4 m sin" ( 2b) +1
(39)
This can be plugged back into the ansatz to get the emulator
h hig.o (-71 + (-4 7 sin ( w) + 71 + 1) ( 2-n
(40)
7.7 sin" (π v) - 4 71 sin (i v) +
A.1.4 Closed-Form Superiority when using the Implicit Scheme as Baseline
We can use the closed-form expression for the found emulator heg to formulate the closed-form superiority. For this, choose the analytical scheme ås as the reference and use the training simulator (the implicit scheme, (x) as the baseline.
This gives in terms of phase
arg(âin ) arg(ds)) (
arg((e) (á.)
arg ((-71 + (-4 71 sin" (I (+) + + 1) elim)))) .e -22m)))) arg (e2imm)))
arg - 21000
arg (e2
(41)
This equation given here to show that there are indeed closed-form expressions. Unfortu- nately, most of them did not amend themselves easy simplification which is why we present further lengthy expressions in the supplemental material in the form of a HTML table: https://tum- pbs. http////uaooouupproooiy/yoolllie
A.1.5
More Combinations of References
The main text only considered the case of the following combination of solvers:
Training on the implicit scheme
Baseline the implicit scheme
Testing the analytical scheme
While testing on the analytical scheme is reasonable since it allows to exactly assess the full error conducted, changing up the training and baseline schemes allows study additional effects.
When we fit on the explicit scheme, we get the following solutions for the complex-valued equation
0n = -
(42)
0 = + 1.
(43)
25

---

Baseline Scheme Implicit Train Scheme Implicit
Baseline Scheme Implicit
Train Scheme Analytic
1.75
1.50
1.25
1.00
0.75
Gamma
0.25
0.00
00.3
Explicit
Analytic
0.03
2.00
1.75
0.003
0.03
1.50
0.3
1.25
1.00
0.75
0.50
0.25
0.00
0.005
0.01
0.05
0.1
0.5
0.005 0.01
0.05
0.1
0.5
Figure 7: All possible combinations of training references and baseline references for the advection superiority when measured in terms of the magnitude error against the analytical scheme. The first column represents training on the implicit scheme, the second column displays training on the analytical scheme. The first row displays using the implicit scheme the baseline reference, and the second row uses the explicit scheme as the baseline. The superiority figure in the top left identical figure 3 presented in the main text in which favorable inductive bias leads to becoming better than the reference. In the bottom left we see that we use the same training reference but baseline against the stronger explicit scheme, almost superiority is achievable under the stability limit, except for 6 2 0.25 in which the implicit scheme was also better than the explicit scheme as shown in figure However, surprisingly,
explicit ansatz becomes superior beyond the stability limit hn| for regions of á around the trained tb. The second column highlights that superiority is even stronger if use the analytical scheme as the training reference.
As expected, since the functional form of the ansatz exactly the same as the training scheme, we fully "re-learn" the training scheme. As such, there can also be no superiority in terms of magnitude or phase, hence we have {(t) = 1. As such, we also will not show any figures for this case.
This leaves with the of fitting against the analytical scheme which gives the following solutions for the complex-valued equation
0, sin (2 π 71 1)) sin (2 1/)
(44)
0 sin (2 1/ (71 + 1)) sin (2 17)
(45)
This yields the fitted emulator
(e ( sin (2 (2. -b (71 + 1)) sin (2 n 7 +)) 2irz
has d
(46)
sin (2 v/)
Using either the implicit or the analytical scheme as the training reference, either the explicit or the implicit scheme as the baseline, and the analytical scheme as the testing reference, we discuss four potential combinations in figure 7.
26

---

Superiority of Implicit- Implicit-Impiicit Analysis over number of test rollout steps for & 0.01 and #= 0.003
Analysis Magnitude -mplici--mppiii
Analysis Phase-mpliictimpllic
2.0
Gamma
-3.0
-0.03
Psi
0.003
10" 10
10' 10" 10a 10"
101
T 10"
10'
10
Figure 8: The superiority achieved in the theoretical advection example shows persistence across multiple test rollout steps. For regimes within the stability limit m 1, superiority converges to after many time steps T > 1, demonstrating sustained performance parity with the baseline. Superiority terms of phase or angle appears persist until numerical issues evaluating the arctan arise. For magnitude- baseed superiority (i.e., g(dl < 1), the advantage is preserved for significant number of practical time steps before eventually equalizing with the baseline.
A.1.6 Influence of the Number of Test Rollout Steps
While the main text focused single-step evaluation scenario, here extend our analysis to multi-step test rollouts. For this, we consider this form of the superiority
10
[+]]
(47)
again with the implicit scheme is the training reference, the implicit scheme is as the baseline, and the analytical scheme áø as the testing reference. We present the results figure 8. While the results show that the superiority in terms of magnitude is lost for large t, it is maintained for a significant number steps. For faster dynamics, i.e., higher [7]., the superiority lost faster.
A.2 Diffusion Equation
We now turn to the one-dimensional linear diffusion equation, also known as the heat equation, given by
duu = vŪ u.
(48)
where 0 is the diffusion coefficient. As with the advection equation, consider unit interval domain (L = 1) with periodic boundary conditions. To streamline the analysis, we introduce the
reparametrized form of
7/2
0,u
Onzu,
(49)
2N
where the dimensionless parameter ?t2 defined as
F22 Nu At.
(50)
In the derivations to follow, the number spatial points N and the the time step At will cancel.
A.2.1 Deriving the Schemes for the Diffusion Equation
The second spatial derivative uzzu is commonly approximated using second-order central finite difference scheme: 0zz (14+1 2u, + u,-1)/Art Ar... Substituting this into the reparametrized equation equation 49, we obtain the semi-discrete form:
Ui+1 2u; + Ui-1
(51)
dt
AN
A
Given that the spatial discretization uses intervals of length Ar = 1/N (since L 1), we have Az? 1/N2. This 2/N2 term in the denominator cancels with the 2 term in the definition of 72 within the fraction, simplifying the equation to:
d
72 (ui+1 2u, + 44-1).
u;
(52)
2At
27

---

Applying forward Euler method for the time derivative, (u!+1) //A, and evaluating the right-hand side at time t, yields the explicit Forward Time Centered Space (FTCS) scheme:
[++1]
u [6]
At 24t (w 2nl u/
72
(53)
The At terms cancel, resulting in the fully discrete update rule:
[++]]
/ + 22 (w - 2u + w!_1)
(54)
This can be written in matrix form for the entire state vector un
[++
72
12/2
12 /
U [c]
[+1]
7/22
72
0
01 [1]
W1
(55)
[++1]
0
12
2
72
72/2
N--1
72/2
72/2
1
72
=AD
Note that the diagonal elements are 1 y2 because the term + (-2u!)) contributes Y2 [t] The
matrix A circulant. Its operation is equivalent convolution with kernel [/2/2, 1 72- 7/ /
(when centered at 1 -22). In Fourier space, the update rule becomes unf [+++]] ê, O ú where the Fourier multiplier é for mode (corresponding to relative mode m/N E [0. 1/2)) is:
és + 72 222 -2 2 t
e 12x6)
(56)
2
1 + 72
c 2)
(57)
1 + (cos (2 1).
(58)
The FTCS scheme for diffusion is stable if &s < 1, which implies & 1 for the worst-case mode (cos(2:)) (2 11).
Implicit Scheme for Diffusion Using a backward Euler method for time discretization (evaluating the spatial derivative at time t + 1) yields the backward-in-time central-in-space (BTCS) scheme:
[e+1] u 2
72 [t+1l
(u,41
2u,
[++1] + u [++]]
(59)
Rearranging gives:
[++1]
72 [++1] 21 (((+)]] u!!')
(u
(60)
The operator acting on up in Fourier space has the multiplier 1 72 (co0((2TT) 1). Thus, the Fourier multiplier for the implicit BTCS scheme is:
[[+1]
is 1 72(c00(2 1)*
(61)
This scheme is unconditionally stable for 7/2 > 0.
Analytical Scheme for Diffusion In Fourier space, the diffusion equation transforms to
2Tm
2
4m ) -
üm.
(62)
L
For L = 1, this simplifies to
d
dt üm uu(Tmm) Pim.
(63)
The solution after one time step At becomes
um
[++1] x|ا! expppppppp (2Tmm) 2tt).
(64)
28

---

Diffusion
Poisson
0.10
0.1
0.05
Unden
0.2
Unist-Smisthing
%--03
q-10
0.00
0.2
10 -00
4-000
-0.09
Difuuin
-0.2
Smoottiotiing
Schame
Explct
FD2
-2.4
rrpick
Ecratle
-0.10
Retatuve Madee
Reletiun Mode
Figure 9: Error analysis for diffusion (a) and Poisson (b) schemes. Similar to the advection analysis in figure the explicit scheme has smaller error when within its stability limit T/2 (a). For the Poisson equation, iterative solutions converge against the curve of the direct inversion with second-order finite difference scheme when the number of iterations g increases (b).
Using the substitution 7/2 22N226, have DÂt -22/(2N2). The exponent can then be rewritten as
- -2(2m)) =
2N2 (2xm²)²
(65)
7/2
4T 2
(66)
2N
2n m2
--12 N2
(67)
So, the analytical Fourier multiplier is
O = exp (-12 N2
2 2m²
exp
(-2
12 (2TT))))
(68)
A.2.2 Errors of the Diffusion Schemes
Using the relative mode ø = m/N, N, the multipliers are
Es = 1 + 72(cos(2TT)) - 1),
(69)
1
is (2(cos((2T)) - 1)
(70)
(-
72 (2Tp)²)
âs exp
(71)
2
The relative magnitude errors are
|å
= (72 (cos
()) + 1) 22 -
(72)
(+
and
1 |xpl
-2-72- sin? (TT Φ)) + 22.m--72.00
(73)
2 2.72 sin (IT Φ)
Since these schemes for diffusion using central differences do not introduce phase errors (their Fourier multipliers real), only magnitude errors are analyzed. These magnitude errors are depicted in figure 9a for various & and 7/..
A.2.3 Closed-form Superiority Expressions for Diffusion
Similar to the advection case, emulator (ansatz) can be trained data from low- cow-fidelity diffusion solver (e.g., the explicit FTCS scheme ég. or the implicit BTCS scheme ig) and then evaluated against the analytical solution dz- Here, we will choose the ansatz
+ 0 (20s((TT)) 1),
(74)
29

---

Baseline Scheme Implicit Train Scheme Implicit
Baseline Scheme Implicit
Train Scheme Analytic
1.75
1.50
1.25
1.00
0.75
Gamma
0.25
0.03
0.00
0.3
Explicit
Explicit
Analytic
3.0
2.00
Psi
1.75
0.003
0.03
1.50
0.3
1.25
1.00
0.75
0.50
0.25
0.00
0.005
0.01
0.05 0.1
0.5
0.005
0.01
0.05
0.1
0.5
Figure 10: Superiority analysis for the diffusion equation across different combinations of training references (columns) and baseline references (rows). The analytical scheme serves as the testing reference for all evaluations, revealing regimes where trained emulators outperform their baseline solvers.
i.e., it has functional form similar to the FTCS scheme but has the 7/2 value as free parameter O. When training on the implicit scheme at relative mode 1o M/N [0, 1/2), we find
7/2
(75)
2 72 sin (π +1
In other words, in contrast to the FTCS scheme, our fitted ansatz finds a "corrected" ?2- It converges 72 if close O. We can plug the found back into the ansatz. Since the ansatz, due to its inductive biases, better approximates ås than the training solver does (especially in regions where the training solver has large errors), superiority can be observed. Detailed expressions for de (the fitted ansatz) and the resulting superiority ratios are provided in the supplemental HTML Various combinations of training and baseline schemes, all tested against the analytical solution, are presented in figure 10.
material: https:/tum--ps. babss.iituu..ioemullaaa
A.3 Poisson Equation
The Poisson equation, - f. serves prototypical example of elliptic partial differential equation (PDE). It is found at the core of many steady-state physical descriptions, such as modeling electrical potentials structural deformations. Furthermore, it frequently appears when integrating systems with constraints, for instance, in the pressure correction step (projection method) for solving the incompressible Navier-Stokes equations.
A.3.1 Deriving the Schemes for the Poisson Equation
consider a unit interval domain (L = 1) with periodic boundary conditions.
Direct Finite Difference (FD) Scheme The second-order spatial derivative can be approximated using a second-order central finite difference scheme. For the Poisson equation, this yields:
U
2u; + Wi-1 fur
(76)
Ar
where Ar 1/N is the spatial step size. This equation represents linear system Apu, f.. which can be solved directly.
30

---

In Fourier space, the operator corresponding to the central finite difference approximation of -02x has a multiplier for mode m (with relative mode & m/ N E [0, 1/2)) given by
-(2 cos (2TT)
2 1 cos()TT))
kFD
(77)
Ar
Ar?
Thus, the equation in Fourier space is k f... Solving for im gives the Fourier multiplier for the direct FD solution:
Ar²
B
(78)
KFD.
2(1 cos(2))))
The solution for each mode is then üm = fm. This method, while using finite differences, provides direct (non-iterative) solution in Fourier space.
Iterative Scheme Based on FTCS (Richardson Iteration) Alternatively, the Poisson equation can be solved iteratively by considering it as the steady-state solution of a pseudo-time-dess time-dependent diffusion (heat) equation: 0,u ux + f- As pseudo-time т + oo, we expect 0,u - 0, which recovers the solution -r f- Applying the Forward Time Central Space (FTCS) scheme to this pseudo-time problem leads to the Richardson iteration. In Fourier space, for mode m (corresponding to relative mode Φ m/M N E 1/2),, the update is
(2(cos(2mg) 1)
(
a + 8,
(79)
A
where 6, is the pseudo-time step and q the iteration count. For optimal convergence (related the FTCS stability for the heat equation), we choose the maximum stable pseudo-time step, 6, = Az² Substituting this into equation 79 yields:
Ar
((cos(2TT) -
i[g+1] m
d.9 + m
ii + fm
(80)
Ar?
cos(2 Ar
fm.
(81)
2
This is linear affine recurrence relation. With an initial guess 600 the solution after q iterations is:
[0] - (cos (2 8 Ar²
i [4] (cos(2
cos (2 2 fm.
(82)
A common initialization is am O. Under this assumption, [o] im where the iteration- dependent Fourier multiplier for this iterative solver
Ar 1 (cos(2
16.4
(83)
2 1
cos(2TT))
In the limit the number of iterations q approaches infinity (with the important caveat that 2n must not equal 0 any integer multiple of π, which is equivalent requiring that the relative mode ø must not equal 0 or any integer multiple of 1/2), we observe that the iterative solver's Fourier multiplier ba.q converges asymptotically to the direct finite difference solution multiplier åg as defined in equation equation 78. This convergence behavior demonstrates that the iterative scheme ultimately produces a solution that matches the solution of the discretized Poisson equation obtained through direct finite difference methods.
Analytical "Pourier-Spectral Scheme While the direct FD solution 3g is exact for the discretized system equation 76, it still contains spatial discretization error compared to the true solution of the continuous PDE - = f. A fully analytical solution method, often termed a Fourier- spectral method, utilizes the exact Fourier representation of the continuous differential operator. The continuous second derivative 0.r acting on a function u(x) transforms to multiplication by (- in Fourier space, where kim 2mm. L is the dimensional wavenumber for mode m. For / = km = 2mm 2T d. Thus, for the continuous Poisson equation -dzx u = f, its Fourier transform is - k-hh )m = k22 ûm = fm. The Fourier multiplier for the "true" " analytical inverse Laplacian operator (-0xx) (using relative mode @ m/N) therefore
1
(Ar))
(2m²)² (for Φ + 0).
(84)
(2TNの)²
mm
31

---

For m 0 (the mean or DC component), fo must be zero for periodic solution to exist (compatibility condition). and úg is typically set to zero or determined by other means (c.g.. ensure zero mean solution). This multiplier d is free from spatial discretization error, assuming is sufficiently smooth and periodic such that its Fourier series converges appropriately.
A.3.2 Errors of the Poisson Schemes
We will use the following versions of the Fourier multipliers
ás
(85)
4 π 0
8
(86)
2 2 cos (2 ②)
0.5 0.5 cos (2 ())
b0.g 1 cos (2- Φ))
(87)
Note that we excluded the factor of Ax? in the definition of all schemes since it just acts as global scaling factor and does not affect the relative errors (and the superiority expressions).
For evaluating the different approximate solution methods, we use the fully analytical (Fourier- spectral) solution åø (from equation 84, with Φ = m//N as the ground truth reference. Let ég denote the direct FD solution multiplier Bg (from equation 78), and Egeg denote the iterative solver multiplier (from equation 83).
The relative magnitude error of the direct FD scheme is:
2 |
π²² x
(88)
sin (π +))
And the relative magnitude error of the iterative scheme (after iterations) is:
[ a
|åol
cos (2 Φ) x² x².
2.0
1.0
(89)
cos (2 π 6)
These errors are plotted in figure 9b. One might also be interested in the error of the iterative scheme relative to the direct FD scheme it is attempting to converge to. This "iterative truncation error" is:
10.9
Ea E
cos (2 π Φ) -
2.0 sin" (a d)
(90)
cos (2 ②)
A.3.3 Closed-form Superiority Expressions for Poisson
As an ansatz, we choose the following functional form:
i 2 cos(2
(91)
where 0 a free parameter. Its functional form is similar to the finite difference direct solver (FD) solution, but has the 0 parameter free parameter. When this emulator is trained on data from an unconverged iterative solver bw.q for a fixed & and training mode y and evaluated against the fully analytical (Fourier-spectral) solution (ā²), find the parameter value
1.0 cos (6.28318530711959 v)).
(92)
Superiority can arise if the emulator better approximates áø. The baseline for comparison is the training data source, i..q. The superiority ratio is then:
|
EPoiss
(93)
2
02
(6.283115337 (7959-9)-110 sin (T-+6)
2.0 . 0²
cney (2 (2)))())
1.0
cos (2 Φ)
32

---

Trired Rerative Andyteieeatte
Socdie MuyyicHb
Sochne
2.0
tantoms
1.22
1.00
0.3
0.25
0.00
0.00
2.01
0.05
0.1
05
0.022 0.21
0.05 0.2
0.01 0.03 0.1
PN
Figure 1I: Analysis of emulator superiority for the Poisson equation when the emulator is trained on data from the iterative (Richardson) scheme (using iterations, denoted pg.q).. The panels explore superiority under different evaluation configurations: The baseline solver (denominator of the superiority ratio, equation 5) is varied, for instance, using the iterative solver itself (&p.y) or the direct Finite Difference (FD) solver (Bx). The high-fidelity reference solver (against which errors measured) also varied, comparing against either the direct FD solution (3g) or the fully analytical Fourier-spectral solution (å+). This highlights how the emulator's performance advantage shifts based on whether is compared to its direct training data or more accurate solver, and whether the ground truth' for error calculation includes discretization errors (FD) is error-free frre (analytical). The impact of the training data's 's convergence level (g) on achieving superiority is also demonstrated.
Trained on Analytic scheme Reference-Bseeline Analylic-FD2
Reference-Saseline Anallyic-lteratiive
2.0
10
100
Psi
0.003
0.5
0.03
0.3
0.005
0.01
0.05
0.005
0.01
0.05
phi
phi
Figure 12: Analysis of emulator superiority for the Poisson equation when the emulator is trained exclusively on data from the fully analytical Fourier-spectral solution (denoted d.p). (Training on the direct Finite Difference (FD) scheme, Bp. is not shown as it would lead to the emulator re-learning that scheme due identical functional forms, precluding superiority against it). The panels explore how superiority, defined by the ratio in equation 5, is observed under different evaluation configurations. The baseline for comparison (denominator of the ratio) is varied, potentially including less accurate iterative scheme (e.g., (g.q) the FD scheme itself (B²). The high-fidelity reference (against which errors of both emulator and baseline are measured) is always the true analytical solution (áæ). This reveals how the choice of baseline and ultimate reference truth influences the observed superiority,
even when the emulator is trained on high-qualit data.
Detailed expressions are provided in the supplemental HTML material. Illustrative examples of such superiority for combinations of training, testing and baseline solvers given in figure 11 and figure 12.
A.4 Additional Derivations
A.4.1 Derivation of the Closed-Form Optimizer for the Linear Case
We aim to show that for the linear simulator and the linear emulator, the optimization problem be solved in closed-form, and that for single-mode training, the solution independent the specific distribution over the mode's amplitude and phase.
Recall the one-step prediction error objective using the mean-squared error (MSE) loss, C(a, b)
N la blli:
33

---

1 Ilfo (u;) - Ph Pn(un))l
(( ~In
(94)
where Un In denotes an initial state sampled from the distribution In. Ph is the low-fidelity solver (our training reference), and fg is our linear emulator.
The key to simplifying this is to move to the Fourier domain. By Parseval's theorem, the L2-norm in state space is related the L2-norm Fourier space by a scaling factor: ||x 2 N [ where is the Discrete Fourier Transform (DFT) of x. Applying this, the loss function becomes:
[(0)) uu I<nn
N [o(up)
Pa (mp))
(95)
By the Convolution Theorem, the cross-correlation defining our linear emulator fo becomes an element-wise product in the Fourier domain. Let (0) be the Fourier multiplier for our emulator and F4 the multiplier for the reference solver Ph. The loss function then:
I(0) Eu, In
N2 (0) Q FA O mnll Eus, In
1 N2 (l(ân(0) 1x) o un
(96)
Now, introduce the crucial assumption from theoretical analysis: the training data distribution I, only produces states with energy in a single Fourier mode, M. This means the Fourier represen- tation ü₁ is sparse vector with only non-zero complex- valued entry, üM, at index M. The I 1.2-norm (which the square root of the sum of squared magnituees) collapses to single term:
1 N2
L(Ø) Ea..
(()r(())
FM)) uml?
(97)
Using the property lab|?= |/| for complex numbers, we can separate the terms:
(( East
N2 (0)
PMMl" | 2
(98)
The term (0) - AM 2 depends only on the parameters and the known properties of the solvers; it is constant with respect to the expectation over the data distribution. We can therefore factor it out of the expectation:
L(0) = N2 14M - Fasl
Eiss
[uAs (2]
(99)
This final form makes the key point clear. The term | is simply the expected energy of the training data at mode M. This is a positive constant determined by the distribution over the amplitude of the mode. Crucially, the phase of üss is entirely absent, and the overall distribution of amplitude only contributes a constant positive scaling factor to the loss. Hence, we can write
L(0) ~ (0) - FM
(100)
To minimize T(()), we only need to minimize the term 1Â (0) P. This minimum achieved and only if:
(0)(0) fm
(0) f.
(101)
Substituting the functional form for our linear ansatz at mode M (with corresponding primitive root of unity WM M/N arrive at the complex- valued equation presented in the main text:
34

---

0o + 8, fM
(102)
This equation can be solved for the two real-valued parameters 8u and 01. This derivation confirms that under the single-mode training assumption, the optimal parameters independent of the specific distribution over the amplitude and phase of that mode. worth noting that if the training data contained multiple modes, the loss would become weighted least-squares problem, where the optimal parameters would depend on the relative energy across the modes.
A.4.2 Proof that One-Step Supervised Training is equivalent to Discrete Residua
We show here that for the linear systems considered in paper, training with PI loss formulated as the discrete PDE residual leads to the same optimal emulator parameters as training with the supervised loss if only a single mode M is active. Therefore, our theoretical insights on superiority hold equally for this common PI training paradigm.
Let consider more general setting than in the main text of general linear, time-invariant [t] numerical solver that advances the state vector u to u
[t+1] according
Au |++]] Bu
(103)
where A and B are matrices representing the discretized spatial and temporal operators. This form covers explicit schemes (where A = I) and implicit schemes. The solver operator is thus Bu Our neural emulator, fe. aims approximate this operator.
(u ( A
We define two loss functions averaged over the training data distribution In
1. One-Step Supervised Loss (Loup): This measures the direct error between the emulator's output and the solver's output.
(0) Eus, ~i [N' (l.fo(u)) (Auu))) Eun-In
[o(()))) A- A- (104)
N
Discrete Residual Loss (Cres): This measures how well the emulator's output satisfies the underlying discrete PDE.
1 [IA.fo(um) Buall
Crs Eug In.
(105)
For linear PDEs on periodic domains considered in our paper, the matrices and are circulant. Circulant matrices diagonalized by the Discrete Fourier Transform (DFT), meaning the matrix operations become simple element-wise multiplications in the Fourier domain. Let and b be the Fourier multipliers (eigenvalues) of A and B, respectively. The Fourier multiplier of the solver Pn is then f = á C b. Let qn (0) be the multiplier for our linear emulator fo-
Using Parseval's theorem, becomes
express the losses in the Fourier domain. The supervised loss
Csup (8) = Ein
N2 [|qn ( un - fr @ 1l
Ean
1 [l()n()) - -.)) o az
(106)
.22
The residual loss can be rewritten by factoring out A
1
(rs(()) = Eun
fo(un) - - -A'Buu Buu un
Eun
A (fo(un) Pa Pa (ux))))
N
N
=iPr
(107)
35

---

where we identified the solver operator. In the Fourier domain, this becomes
Cres (0) Ein
N2 a
(0) un fn 1
Bag
á j (ān(0)
Un
(108)
Comparing the loss functions, see that Lres. is spectrally weighted version of Cup.. There three situations to consider:
1. The reference simulator an explicit scheme, i.e., A = I. As such à will be vector full of ones and the objectives are identical, Saup (0) = (0), leading to identical minimizers 6".
2. The initial condition only has a single mode active. Similarly to how we argued in sec- tion 1. the L2-norm collapse to single term and the relevant entry âM is simply a constant scaling factor in the optimization with respect to 0. Consequently, we have xx Les (0) and the optimizers l are identical.
Ssup (0)
3. In the most general case of an implicit scheme A + and 2 modes being present, supervised loss and residuum loss cannot easily be related each other. In alignment with section 1.4.1 the optimizers 0* of both objectives depend on the distribution of energy in the modes. It must be noted that the residuum-basee formulation additionally scales (and rotates) the energy in the modes of úU due to the element-wiss muliplication with à.
In conclusion, under the assumption of only a single mode being present, our superiority analysis holds equally for both supervised and physics-informa training.
B Experimental Details
Our experiments are implemented in JAX (Bradbury et al., 2018) and build upon components the APEBench suite (Koehler et al., 2024).
Architectures particular, we use the following neural architectures:
ConvNet: A feedforward convolutional network with 10 hidden layers of 34 channels each. Each layer transition except for the last uses the ReLU activation. The effective receptive field 10 + 11 and the network has 31757 learnable parameters.
UNet: A classical UNet with 2 hierarchical levels (i.e., three different spatial resolutions) and a base width of 12 hidden channels at the highest resolution. It uses double convolution blocks with ReLU activation (preceded by group normalization with one group) per level. The spatial resolution is halved at each of the levels while the channel count doubles. Skip connections exist between the encoder and decoder parts. This configuration has 27193 learnable parameters and an effective receptive field of 29 per direction.
Dilated ResNet: A Dilated Residual Network with 2 blocks, 32 hidden channels channels, and ReLU activation. Each block uses series of three convolutions with dilation rates 1, 2, and 1. Each convolution is followed by group normalization (with one group) and then the ReLU activation. This configuration has 31777 parameters and an effective receptive field of 20 per direction.
FNO: A Fourier Neural Operator with 4 blocks, 18 channels, and 12 active Fourier modes. Each block consists of spectral convolution and point-wise linear bypass, with the GELU activation applied to their sum. Lifting and projection layers are point-wise linear (Ixl) convolutions. This configuration has 32527 parameters and an infinite receptive field due to the global nature of Fourier transforms.
Vanilla Transf ormer: transformer architectue considereing each spatial degree of freedom a token and performing dense self-attention over them. We consider a setting with 31669 parameters using 28 hidden channels across 4 transformer blocks (each consisting of dense self-attention in space and two-layered multi-layer perceptron whose inner dimension twice as large as the number of hidden channels). Due to the global nature of the attention, the effective receptive field is infinite. However, in contrast to the convolutional
36

---

or Fourier-spectie based architectures, the Transformer does not have an inductive bias for the periodic boundary conditions.
All the architectures have trainable parameter numbers in the range of 27-33k to allow for a fair comparison.
Metries We train the emulators using the mean- squared error (MSE) as training metric and the normalized nott-eea-ssqaree error (nRMSE) test metric.
B.1 Nonlinear Emulators on Linear Advection
In this section, detail the experimental setup for the nonlinear emulators applied to the linear advection equation, as presented in section
Experimental Variables and Training Configuration The experiments systematically varied three key factors: the combined variable (which is proportional to the CFL number), the training (and corresponding baseline) simulator, and the initial condition distribution for testing. Specifically, we explored the following combinations:
Autoregressive Generalization against implicit scheme:
71 = -3.0
training simulator: implicit first-order upwind
training IC distribution: only first mode active
- test IC distribution: only first mode active
Autoregressive Generalization against explicit scheme:
71 = -0.9
training simulator: explicit first-order upwind
training IC distribution: only first mode active
test IC distribution: only first mode active
State-Space Space and Autoregressive Generalization against implicit scheme:
71 = -3.0
training simulator: implicit first-order upwind
training IC distribution: only first mode active
- test IC distribution: modes one to five active
These combinations were strategically chosen to investigate both autoregressive superiority (using the same initial condition distribution for training and testing) and state-space superiority (using different initial condition distributions). In all cases, the testing reference was the fully analytical solver. The domain is discretized into N 100 points.
Training and Evaluation Setup The training trajectories consist of 300 initial conditions rolled for I time step. Hence, the training dataset contains 300 samples with an input and output frame. The emulators were trained using a mean squared error (MSE) loss function to predict single step ahead. For evaluation, we used 10 initial conditions rolled out for the 50 time steps shown in the main text.
The optimization process utilized the Adam optimizer (Kingma and Ba, 2015) with batch size of 32. We used a warmup cosine decay learning rate schedule featuring warmup phase followed by cosine decay, with total of 10,000 training iterations. The learning rate ramped up from 0 to 1 x 10-3 during the first 1,000 steps (warmup phase), then gradually decayed to zero following a cosine schedule for the remaining steps.
Seed Study for Robustness We repeated each experiment 12 times with different random seeds. Results the main text show the median performance across these runs along with the 50% percentile intervals.
37

---

Computer Resources We conducted this experiment on a single NVIDIA RTX3060 GPU. Running for all the required seeds took less than 2 hours.
B.2
Burgers Superiority
B.2.1
Implicit Time Integrators for the Burgers Equation
The Burgers equation on the one-dimensional unit interval in non-conservative form with periodic boundary conditions reads
du Ou
8u dx2
+ u dx
u(t, u(t. 1).
(109)
Ōt
Similar the linear PDES in section A. we will discretize the domain into N intervals of equal length Ar = 1/N. Then, we consider the left end of each interval a nodal degree of freedom and collect the values of these points into the vector u In contrast the linear PDEs, the solver described below will exclusively work in state space.
The matrix associated with the discretized second derivative in one dimension reads
1
0
... 0
1
-2
1
0
0
1
0
-2
0
0
Lh
...
(110)
(Ar)²
0
-2
it is not exclusively tri-diagonal but also has entries in the top right and bottom left corners.
The convection term requires special treatment because of its nonlinearity and the advection charac- teristic of the first derivative. Let F, and B₁ represent the forward or backward approximation of the first derivative in one dimension on periodic boundaries, respectively, via
0
0
0
0
1
1
0
0
0
F2
0 -1
B₁ = Ar
-1
1
0
0
Aw
0 0
0
0
1
(111)
Again, note the element in the corner entries of the matrices. differentiation matrix based on the winds wh
Then, we can build an upwind
8 (=-(WW)) Wh
81 (1w)) o)
+ w/
(www)
diag max
BA + diag
max
Fh.
2
2
positive winds
negative
(112) Deducing the positive and negative winds from neighboring averages (using the periodic forward shift 8-1 and backward shift s1 operators) is necessary to have correct movement if the winds change sign over the domain. If we use the discrete state vector un winds wh, discretize the continuous equation via the method of lines as
dua
+ Pa(ux)us vLa.u
(113)
dt
Naturally, the spatial discretization of nonlinear PDE leads to system of nonlinear ODEs. If we treat it fully implicitly. here via a backward Euler in time for the time discretization, we get
[e+ I] 1
[c] To (u
(e- (e-"))yu, [++1]] "llyu(t vL u
[+1]]
(114)
At
In other words, in order to advance from u to u [[+1] (which is what the numerical simulators Ph, and eventually the neural emulator, will do), need to solve nonlinear algebraic system of
38

---

equations. The structure of the nonlinearity allows for Picard iteration (Turek, 1999) (which can be viewed quasi-Newton method) that introduces an additional iteration index [k] to the state vector at the next point in time
++ [+]]+++]] u + Tx (u)
(u] )
( [++1]]++1]]
vL / [++]][++]]
(115)
A
[t+1][A]] As such,
The key idea is that the upwiding matrix is linearized around the previous iterate ui one has to solve linear system for 2 [++]]]++]] at each Picard iteration step, which reads
I MTT (( ( [talkk) Atv Ln 2/
[[+]][k-]]
(116)
= Au
Note that across the Picard iterations, the right-hand side is constant, i.e., independent of the iteration index. This because it is always the previous state in time w/ The system matrix A₁, on the other hand, is dependent on the previous iteration index because the upwinding matrix Th (x ( to be reassembled at each iteration step. We use N = 60 for the spatial discretization. As such, the matrices are small and we employ direct solver via LU decomposition.
[++1][[8) has
While convergence of the exedd-poin iteration over [k] is desirable to mimize the nonlinear residual between two consecutive time steps, common approximation is to perform only single Picard iteration per time step. This leads to a linear system of the form
I + AtT (u)) Atv
Un
[++1]
U
[c]
(117)
=A, AA,=A,("!"))
We call this the one-step (P1) implicit method.
In our experiments, we use the one-step (PI) implicit method as the training reference (and the baseline during superiority evaluation) for the neural emulator. This method linearizes the nonlinear term by evaluating it at the previous time step and performs only a single matrix assembly and linear system solve per time step. As result, the nonlinear residual is not fully converged, which leads to errors in shock propagation, especially as the solution develops sharp gradients. This is qualitatively depicted in figure in the main text.
For evaluation, we use a fully converged nonlinear implicit method as the testing reference. This method iterates between assembling the system matrix (which depends the current iterate of the solution) and solving the resulting linear system, repeating until the nonlinear residual falls below a specified tolerance which set to 10- This approach yields the correct shock propagation and serves as a high-fidelln fidelity reference for comparison.
In both the training and testing references, the linear systems are solved using LU decomposition, ensuring solutions are obtained to machine precision. While we use direct solvers here, employing an unconverged iterative linear solver (e.g.. fixed number of Krylov or Jacobi iterations) could introduce additional sources of error and is a potential avenue for further investigations into emulator superiority.
As initial conditions for both training and testing, we use state with both the first Fourier mode and the mean/zero mode active. In figure 13, we show the effect of the unconverged Picard iteration on the solution.
B.2.2 Experimental Details
Experimental Setup We conducted our experiments using the UNet architecture described before. discretized version of the Burgers equation with viscosity coefficient of 1/2 = 0.2 and a convection difficulty parameter & = -8.0, representing a strongly nonlinear regime. The spatial domain was discretized using 60 grid points.
The model was trained
Training and Testing Simulators For training, we use the P1 method which truncates the fixed- point iteration after one iteration each time step. This represents our low-fidelity simulator. For testing and evaluation, we used the same scheme but with full convergence, ensuring full convergence of the nonlinear system and thus serving our high-fidelity reference.
39

---

Solution evolution
1.0
Space
0.4
0.2
0.0
One Step Relative Error P1 to converged
Picard Iterations until Convergence
50
10a
Picard Iterations
One-Step Error
10
20
10
10
10
10th
15
25
30
Time step
Figure 13: sample trajectory from the given initial condition distribution evolved with the converged Picard method (top figure). In the lower figure, the blue curve despicts the number of Picard iterations (consisting of linear system assembly and direct LU solve) it requires converge the fixed-point iteration to tolerance of 10-5, The data points are places between two consecutive snapshots to indicate that they relate to the state transition. Requiring > Picard iterations to converge while only performing one in the case of the PI method incurs an error in terms of an unconverged nonlinear residual. Here, we measure one-step error against the converged solution (depicted in red on logarithmic scale). Shaded areas around the blue and red curve indicate the 50 percentile interval around the median out of 300 different initial conditions. Both the number Picard iterations and the one-step error peak at the shock formation point and then go down due to the diffusive nature of problem.
Initial Conditions Both training and testing utilized Fourier-based initial conditions, focusing on the first Fourier mode with randomized phases and an additional offset (i.e., non-zero energy in the meen-moode/zer-moodle This approach generates smooth initial conditions that subsequently develop shocks through the nonlinear dynamics of the Burgers equation.
Training Protocol The network was trained using a one-step prediction framework (temporal horizon train 1) with mean squared error (MSE) loss function 500 training samples. For evaluation, we performed autoregressive rollouts over 20 time steps (temporal horizon test 20) on 100 test samples, measuring performance with normalized root mean squared error (nRMSE).
Optimization Strategy We employed the Adam optimizer (Kingma and Ba, 2015) with batch size of for 100,000 iterations. warmup cosine learning rate schedule was used, starting from 0.0 and reaching a maximum of 3e-4 after 5,000 warmup steps, followed by cosine decay for the remaining iterations.
Statistical Robustness To ensure statistical significance of results, conducted thorough seed study with 30 different network initialization seeds, while maintaining consistent training and testing data through fixed seeds for those processes. This approach allowed us quantify the variability in network performance while isolating the effects of initialization.
Computer Resources We conducted this experiment on a single NVIDIA RTX3060 GPU. Running it for all the required seeds took less than 5 hours.
40

---

C Ablational Experiments
C.1 Multi-Mode Case Study
If the training data contains multiple modes, the closed-form solution for the emulator parameters no longer a simple equation for single mode as shown in section A.4.1. I. Instead, it becomes a linear least-squares problem weighted by the concrete distribution of complex-values coefficients of ú. While an analytical solution still exists for this linear regression problem, the analysis becomes more complex because of this data dependency.
To provide an intuitionn we redid the theoretical experiments numerically (training on a combination of modes and only testing mode 5) with results displayed in table 1. We found that there is still superiority when more modes are present in the training initial condition as long as they are lower than what is being tested on (confirming again what called "forward superiority" in section 3.1).
Table The achieved state-space superiority in the linear scenario of section 3.1 (here measured on m 5) depends on how the energy is distributed in the data has seen during training. We still see "forward superiority" if the training modes are smaller than the testing mode.
Modes in Train IC Dist
Superiority {
1
0,57
1,2
0.61
1,2,3,4
0.72
1,2,3,4,5
0.81
1,2,3,4,5,6
0.87
1,2,3,4,5,677
1.00
1-10
1.10
1-20
2.30
4,6
1.10
3,4
0.79
C.2 Impact Receptive Field in CNNS and Active Fourier Modes in FNO
Our results in figure 5c already hint that local architectures (ConvNet) can achieve stronger state-space superiority than global ones (FNO, Transformer). We hypothesize this is because a constrained, local inductive bias acts powerful regularizer, preventing the model from overfitting to the global error patterns of the fow-fidelity solver.
To provide direct evidence and practical guidance, we performed new ablation study on the linear advection case, systematically varying the effective receptive field of the ConvNet. We present the resuls in table 2. The configuration with receptive field of 11 corresponds to the model used in the main paper. The key finding that superiority is maximized when the receptive field is appropriately matched to the physical characteristics of the problem (in this case, dictated by ? which is proportional to the CFL number). achieving superiority ratio g(i) of 0.73 after 10 time steps. Receptive fields that are too small fail to capture the necessary physics, while those that are too large provide excess capacity that diminishes the superiority effect by allowing the model to learn more of the solver's non-physical behavior.
receptive field 4 yields the best performance,
We also conducted an ablation study on the Fourier Neural Operator (FNO) by varying its number of active Fourier modes. The results of table 3 reveal clear trade-off between model capacity and the implicit regularization that drives superiority. When severely mode), the FNO fails to learn the problem. The optimal performance achieved active modes, which provides just enough capacity to model the target dynamics while the FNO's s inherent spectral truncation acts as a powerful regularizer, filtering out numerical artifacts from the solver and leading to strong superiority (ratio of 0.59). As the number of active modes increases further (3+), the FNO gains the capacity to partially overfit to the solver's structured errors across a wider frequency band, which, while still allowing for superiority, diminishes its magnitude. This demonstrates that superiority is maximized when the model is expressive enough to capture the core physics but constrained enough to regularize away the training data's flaws.
(I active
41

---

Table Influence of receptive field of the CNN on the achieved superiority ratio. This ablation is based on figure 5(c).
Superiority gld
at Time Step
Rec. Field/Time Step
1
2
3
5
6
7
8
9
10
4.0
4.0
3.9
3.9
3.9
3.8
3.8
3.7
3.7
3.6
2
1.0
1.0
1.0
1.0
1.0
1.0
1.0
1.0
1.0
1.0
3
1.0
0.95
0.91
0.88
0.86
0.84
0.83
0.81
0.80
0.80
4
1.0
0.94
0.90
0.86
0.83
0.80
0.78
0.76
0.75
0.73
5
1.0
0.94
0.90
0.86
0.83
0.81
0.79
0.77
0.76
0.74
7
1.0
0.96
0.93
0.91
0.89
0.87
0.86
0.85
0.84
0.83
11
1.0
0.96
0.93
0.91
0.89
0.87
0.86
0.84
0.83
0.82
Table Influence of active modes in FNO on the achieved superiority ratio. This ablation is based on figure 5(c).
Superiority glt)
at Time Step
Modes/Time Step
2
4
6
8
10
8.16
8.06
7.93
7.76
7.58
7.37
7.14
6.89
6.62
6.35
2
1.00
0.93
0.87
0.81
0.76
0.72
0.68
0.65
0.62
0.59
1.00
0.97
0.94
0.92
0.90
0.88
0.86
0.84
0.82
0.81
4
1.00
0.98
0.97
0.96
0.94
0.93
0.92
0.92
0.91
0.90
8
1.00
0.97
0.94
0.91
0.88
0.86
0.84
0.82
0.80
0.78
12
1.00
0.98
0.96
0.95
0.93
0.92
0.91
0.90
0.89
0.88
C.3 Effect of Coarse Solver Fidelity by ablating the number of Picard Iterations
In section 5, we argue that the coarse solver needs to be sufficiently inaccurate to enable superiority. Naturally, the question arises by how much. While a universal, quantifiable threshold is likely problem-dependent, we can demonstrate the relationship between solver fidelity and the potential for superiority. We have conducted an ablation study for the number of Picard iterations the Burgers example of section 4.2. We again used the UNet and display the superiority ratio rollout when trained and baselined against Picard solver truncated after certain number of iterations. single iteration represents our original low-fidelity setup, while a higher number of iterations produces more accurate, higher-fie higher-fidelity fidelity solver.
Table The effect of the number of Picard iterations on the autoregressive superiority achieved in the Burgers experiment.
Picard Iter/Time Step
1
3
5
7
(Lowest Fidelity)
1
0.91
0.85
0.88
0.92
0.96
1
1.05
2
0.97
0.98
1.01
1.09
1.16
1.25
3
1
1.05
1.05
1.11
1.19
1.27
1.37
1.47
1.36
1.34
1.46
1.53
1.58
1.67
1.84
These results clearly demonstrate our central thesis: the potential for superiority is directly related to the magnitude of the structured error in the training data. When the training solver is very coarse, there is significant room for the emulator to learn a more regularized operator, resulting in a strong superiority effect (gl4) ratio dropping to 0.85). As the fidelity the training solver increases, structured error in the training data decreases, and the potential for superiority vanishes. When trained on a reasonably converged solver, the emulator simply learns to replicate its already accurate behavior, and superiority is observed.
42

---

C.4 Additional Architectures for Burgers Study
To provide more complete picture, we extended our analysis the nonlinear Burgers' equation experiment of section 4.2 to include all architectures from our study and present the results in table (median superiority out of 30 seeds). The results robustly confirm that the superiority phenomenon is not limited to single architecture but is instead closely tied to the inductive biases of the models. The key finding is that architectures with spatial or spectral inductive bias (ConvNet, Dilated ResNet, FNO, UNet) all achieve significant superiority, outperforming the low-fidelity solver they were trained on. In stark contrast, the vanilla Transformer, which lacks any inherent spatial or spectral bias and treats the input sequence of tokens, fails to achieve superiority and becomes progressively worse than the baseline solver.
Table 5: Achieved superiority for Burgers emulation. This extends the experiment of figure 6.
Superiority g(4 at Time Step
Arch/Time Step
2
4
6
8
10
ConvNet
1.00
0.97
0.87
0.82
0.80
0.80
0.80
0.81
0.81
0.81
Dilated ResNet
1.00
0.91
0.79
0.73
0.69
0,66
0.65
0,67
0.66
0.65
FNO
1.00
1.01
0.91
0.87
0.86
0.89
0.93
1.00
1.08
1.12
Transformer
1.00
1.06
1.13
1.24
1.35
1.45
1.52
1.58
1.63
1.68
UNet
1.00
0.90
0.86
0.88
0.91
0.96
1.00
1.04
1.07
1.11
43


Reviews:
[{"review_id": "WryMAwbNLT", "reviewer": "Reviewer_s59z", "summary": "This paper introduces emulator superiority which defines the phenomenon that machine learning models trained on low-fidelity solver may outperform the data-generating solver in terms of several metrics that measure autoregressive generalization, state-space generalization as defined in the paper. This emulator superiority is brought about by two primary sources: smaller accumulation error from multi-step rollouts and the inductive bias of the specific neural network used as an operator solver.", "questions": "(1) Do any variations in the model contribute to higher superiority? Such as the kernel size or the number of modes in FNO. These experiments may provide guidance when solving a new problem not shown in this paper, such as Navier-Stokes.\n\n(2) Real-world observations often contain structured noise. Assuming we use such data rather than idealized solvers, what would this mean to this superiority?", "limitations": "Yes", "rating": 5, "confidence": 3, "soundness": null, "presentation": null, "contribution": null, "strengths": "", "weaknesses": "", "quality": 3, "clarity": 4, "significance": 3, "originality": 3, "strengths_and_weaknesses": "Strength: \n\n- This paper clearly illustrates experimental motivations, preliminaries and experimental settings, helping readers to follow. The granularity of discussions seems well spread over the paper --- where main descriptions are in the manuscript and in-depth mathematical descriptions are appropriately deferred to appendix.\n\n- Motivations are well described and strong – that I agree the conventional wisdoms goes that the model is as good as the data (but not better). However, recent empirical findings suggest strong extrapolating capacities of neural networks especially in predicting high frequency modes given low-fidelity data, which calls for the analyses over the underlying reasons of such.\n\n- This paper seems novel in that it differs from prior approaches which explicitly trains a model to perform superresolution or extrapolation.\n\nWeaknesses: \n\n- Limited nonlinear problems: would bear larger contribution if the superiority is met across a more complex and harder prediction tasks.\n\n- While several potential error source are enumerated, such as spatial discretization errors, linearization errors, etc, it seems impossible to dissect exactly why and where this emulator superiority arises seems a daunting task. E.g., what inductive bias of ConvNet / ResNET / FNO / Transformer constitute for what portion of the superiority is hard to disentangle. In such case, the practical value of this work seems questionable.", "was_revised": true, "final_justification": "My concerns and questions are well explained and also I reviewed other reviewers' points. So, I increased my score to accept.", "rebuttal": "We sincerely thank the reviewer for their thoughtful and constructive feedback.\nTheir recognition of our paper's strong motivation, clear presentation, and\nnovelty is very encouraging. It is especially gratifying to find agreement on\nthe importance of challenging the \"model is as good as the data\" paradigm. We\nappreciate the opportunity to address the reviewer's comments and believe the\nresulting discussion and additions will significantly strengthen our manuscript.\n\n---\n\n> Limited nonlinear problems: would bear larger contribution if the superiority\n> is met across a more complex and harder prediction tasks.\n\nWe agree with the reviewer that demonstrating superiority in more complex\nsettings is an important future direction. Our decision to focus on foundational\nlinear PDEs and the Burgers' equation was a deliberate methodological choice\naimed at achieving maximum clarity and rigor.\n\nComplex systems, such as the Navier-Stokes equations, involve multiple, coupled\nsources of numerical error (e.g., from discretization, linearization, and\npressure-correction steps). In such settings, isolating and analyzing the\nspecific interactions that give rise to emulator superiority would be\ndifficult. By starting with simpler, well-understood systems, we\nwere able to provide a rigorous theoretical proof (Section 3) for the existence\nof superiority. This controlled environment allowed us to precisely attribute\nthe phenomenon to the interplay between an emulator's inductive bias and the\nstructured, mode-dependent errors of the low-fidelity solver (as shown in Figures\n1 and 2), providing a clear and interpretable foundation for the concept.\n\nWe believe that this stepwise approach—beginning with canonical PDEs and\ngradually moving toward more complex, multi-physics systems—is essential for\nbuilding a robust understanding of emulator superiority. Importantly, many\ncomplex PDEs, such as the Navier-Stokes equations, are fundamentally composed of\nthe same building blocks we have investigated: advection, diffusion, and\nnonlinear interactions. By rigorously demonstrating the superiority phenomenon\nacross these foundational components, we provide strong evidence that the\nprinciple is not an isolated curiosity, but rather a fundamental aspect of\nlearning from numerical data that is likely to persist in more complicated,\nreal-world systems. While additional sources of error and coupling may introduce\nnew challenges, it is not unreasonable to hypothesize that our insights will\ngeneralize to these richer domains. We are committed to extending our analysis\nto such settings in future work, and see our current results as a necessary and\nrigorous foundation for that endeavor.\n\n\n> While several potential error source are enumerated, such as spatial\n> discretization errors, linearization errors, etc, it seems impossible to\n> dissect exactly why and where this emulator superiority arises seems a\n> daunting task. E.g., what inductive bias of ConvNet / ResNET / FNO /\n> Transformer constitute for what portion of the superiority is hard to\n> disentangle. In such case, the practical value of this work seems\n> questionable.\n\nThis is a crucial point that connects to our choice of problems. The primary\npractical value of our work lies in formally identifying and proving the\nexistence of the emulator superiority phenomenon itself. Our linear examples\nwere purposefully chosen precisely because they do allow us to disentangle the\nsources of superiority in a controlled manner. The practical implication is\nimmediate: our findings challenge the fundamental assumption that an emulator is\nstrictly bound by its training data's fidelity. This has direct consequences for\nhow the community benchmarks and evaluates neural emulators. As we argue (lines\n342-347), using flawed numerical references can unfairly penalize emulators that\nhave correctly learned more physically plausible dynamics. Our work provides the\nvocabulary and framework to understand this paradoxical effect. While\ndisentangling every factor in a highly complex system remains a future\nchallenge, establishing the core principle on a solid theoretical and empirical\nfoundation is, we believe, a significant and practical contribution.\n\n> Do any variations in the model contribute to higher superiority? Such as the\n> kernel size or the number of modes in FNO. These experiments may provide\n> guidance when solving a new problem not shown in this paper, such as\n> Navier-Stokes.\n\nThis is an excellent question that gets to the heart of how architectural\nchoices influence superiority. Our results in Figure 4c already hint that local\narchitectures (ConvNet) can achieve stronger state-space superiority than global\nones (FNO, Transformer). We hypothesize this is because a constrained, local\ninductive bias acts as a powerful regularizer, preventing the model from\noverfitting to the global error patterns of the low-fidelity solver.\n\nTo provide direct evidence and practical guidance, we performed a new ablation\nstudy on the linear advection case, systematically varying the effective\nreceptive field of the ConvNet. The key finding is that superiority is maximized\nwhen the receptive field is appropriately matched to the physical\ncharacteristics of the problem (in this case, dictated by the CFL number). A\nreceptive field of 4 yields the best performance, achieving a superiority ratio\n$\\xi^{[t]}$ of 0.73 after 10 time steps. Receptive fields that are too small\nfail to capture the necessary physics, while those that are too large provide\nexcess capacity that diminishes the superiority effect by allowing the model to\nlearn more of the solver's non-physical behavior.\n\nThe table below shows the superiority ratio rollout for ConvNets with different\nreceptive fields. The configuration with a receptive field of 11 corresponds to\nthe model used in the main paper.\n\n|Rec.Field/TimeStep|1|2|3|4|5|6|7|8|9|10|\n|-----------------:|--:|--:|--:|--:|--:|--:|--:|--:|--:|--:|\n|1|4.0|4.0|3.9|3.9|3.9|3.8|3.8|3.7|3.7|3.6|\n|2|1.0|1.0|1.0|1.0|1.0|1.0|1.0|1.0|1.0|1.0|\n|3|1.0|0.95|0.91|0.88|0.86|0.84|0.83|0.81|0.80|0.80|\n|4|1.0|0.94|0.90|0.86|0.83|0.80|0.78|0.76|0.75|0.73|\n|5|1.0|0.94|0.90|0.86|0.83|0.81|0.79|0.77|0.76|0.74|\n|7|1.0|0.96|0.93|0.91|0.89|0.87|0.86|0.85|0.84|0.83|\n|11|1.0|0.96|0.93|0.91|0.89|0.87|0.86|0.84|0.83|0.82|\n\nWe also conducted an ablation study on the Fourier Neural Operator (FNO) by\nvarying its number of active Fourier modes. The results reveal a clear trade-off\nbetween model capacity and the implicit regularization that drives superiority.\nWhen severely under-parameterized (1 active mode), the FNO fails to learn the\nproblem. The optimal performance is achieved at 2 active modes, which provides\njust enough capacity to model the target dynamics while the FNO’s inherent\nspectral truncation acts as a powerful regularizer, filtering out numerical\nartifacts from the solver and leading to strong superiority (ratio of 0.59). As\nthe number of active modes increases further (3+), the FNO gains the capacity to\npartially overfit to the solver's structured errors across a wider frequency\nband, which, while still allowing for superiority, diminishes its magnitude.\nThis demonstrates that superiority is maximized when the model is expressive\nenough to capture the core physics but constrained enough to regularize away the\ntraining data's flaws.\n\n|Act. Modes/Time Step|1|2|3|4|5|6|7|8|9|10|\n|---------------:|-----:|-----:|-----:|-----:|-----:|-----:|-----:|-----:|-----:|-----:|\n|1|8.16|8.06|7.93|7.76|7.58|7.37|7.14|6.89|6.62|6.35|\n|2|1|0.93|0.87|0.81|0.76|0.72|0.68|0.65|0.62|0.59|\n|3|1|0.97|0.94|0.92|0.9|0.88|0.86|0.84|0.82|0.81|\n|4|1|0.98|0.97|0.96|0.94|0.93|0.92|0.92|0.91|0.9|\n|8|1|0.97|0.94|0.91|0.88|0.86|0.84|0.82|0.8|0.78|\n|12|1|0.98|0.96|0.95|0.93|0.92|0.91|0.9|0.89|0.88|\n\n> Real-world observations often contain structured noise. Assuming we use such\n> data rather than idealized solvers, what would this mean to this superiority?\n\nThe reviewer raises an important point about the applicability of our framework\nto real-world data. Our formal definition of \"emulator superiority\" (Equation 5)\nfundamentally relies on the existence of three distinct components: a\nlow-fidelity solver (the training source), a high-fidelity solver (the reference\n\"truth\"), and the emulator.\n\nIn the case of real-world observations, there is no explicit \"low-fidelity\nsolver,\" only the data-generating process of nature, and the \"high-fidelity\"\nground truth is unknown. Therefore, our superiority analysis, as framed, does\nnot directly apply.\n\nHowever, we can address the spirit of the question by considering a modified\nscenario where perfect synthetic data is corrupted by structured noise. In this\ncase, the emulator's ability to \"outperform\" its noisy training data would\ndepend on its capacity for denoising. This relates to the standard machine\nlearning insight that models can learn robust representations from noisy\nexamples, but it is a phenomenon that is distinct from correcting the structured\nnumerical errors that are the focus of our paper.\n\nAnother scenario would be that only initial conditions are corrupted by\n(observational) noise, while the dynamics themselves are noiseless. In such\ncases, the impact on superiority depends critically on the structure and\nintensity of the noise. If the noise significantly alters the spectral content\nof the initial state, it may obscure the clear separation of modes or patterns\nrequired for the emulator to achieve state-space or autoregressive superiority.\nIn particular, strong or broadband noise can mask the underlying signal, making\nit difficult for the emulator to distinguish between true dynamics and\nnoise-induced artifacts, and thus potentially undermining the conditions\nnecessary for superiority to emerge. However, in such a setting we also believe\nthat learning a good emulator will be generally challenging.\n\nWe believe this is a fascinating avenue for future work. We will add a\ndiscussion to the Limitations & Outlook section."}, {"review_id": "6fvtxnBKU9", "reviewer": "Reviewer_k6sC", "summary": "This paper studies a surprising phenomenon in neural PDE modelling where a neural network solver trained only on low-fidelity numerical solver can outperform that solver in accuracy when evaluated against high-fidelity simulation results. The author conducts analysis on linear and non-linear PDEs with both linear and non-linear ansatz like UNet.", "questions": "Can the authors provide a general intuition for what kinds of non-linear PDEs are more amenable to neural networks with limited receptive fields (e.g., feedforward ConvNets), and which ones may require larger context (e.g., long-range dependencies in turbulent flows or stiff PDEs)?", "limitations": "The authors acknowledge one major limitation in the paper, the solver needs to be sufficiently “coarse” for neural emulator to exhibit superiority.", "rating": 4, "confidence": 3, "soundness": null, "presentation": null, "contribution": null, "strengths": "", "weaknesses": "", "quality": 3, "clarity": 3, "significance": 3, "originality": 3, "strengths_and_weaknesses": "Strengths:\n\n1. The study of whether a neural solver trained on low-fidelity data can outperform the source is novel and important to scientific ML.\n\n2. The analysis for linear PDEs is clear and provides intuition for more complex non-linear cases. They show how superiority arises by analyzing linearized spectral errors and optimal filters. \n\nWeakness:\n\n1. The experiments and analysis are limited to 1D PDEs, which are far from the complexity of most real-world applications. The analysis for linear PDEs relies heavily on mode-wise error decomposition, which works for 1D linear PDEs due to the orthogonality and separability of modes. It remains unclear whether the analysis will extend to more realistic scenarios involving nonlinear, multi-dimensional PDEs, such as those found in fluid dynamics, where mode interactions, geometry are significantly more complex and high-fidelity simulations are expensive or difficult to obtain.\n\n2. The analysis of nonlinear neural architectures (e.g., UNets, Transformers, FNOs) is relatively scarce. It could be interesting to see more in-depth discussion regarding the relationship between effective receptive field of model and its rollout error / superiority.", "was_revised": true, "final_justification": "The paper studies an interesting topic and provides clarification in the rebuttal stage. Despite certain limitations related to problem scope and the novelty of the theoretical analysis, I lean towards acceptance and keep my original rating.", "rebuttal": "We thank the reviewer for highlighting the novelty of our core question and for\nfinding our analysis of linear PDEs clear and intuitive. It was precisely our\ngoal to build a solid, understandable foundation that could inform research on\nmore complex systems. The feedback is highly constructive, and we appreciate the\nopportunity to clarify our contributions and strengthen the paper by addressing\nthe points raised.\n\n---\n\n> The experiments and analysis are limited to 1D PDEs, which are far from the\n> complexity of most real-world applications. The analysis for linear PDEs\n> relies heavily on mode-wise error decomposition, which works for 1D linear\n> PDEs due to the orthogonality and separability of modes. It remains unclear\n> whether the analysis will extend to more realistic scenarios involving\n> nonlinear, multi-dimensional PDEs, such as those found in fluid dynamics,\n> where mode interactions, geometry are significantly more complex and\n> high-fidelity simulations are expensive or difficult to obtain.\n\nThe reviewer raises a very important point regarding the generalization of our\nfindings. We agree that extending this analysis to multi-dimensional, nonlinear\nsystems is the critical next step. Our choice to focus on 1D PDEs was a\ndeliberate methodological decision to ensure rigor and clarity in establishing\nthis new phenomenon.\n\nThe 1D setting, especially for linear PDEs, allows for exact analytical\ntreatment via Fourier analysis. This enabled us to provide a closed-form proof\nfor the existence of superiority (Section 3), isolating the core mechanism\nwithout confounding factors from complex geometries or multi-dimensional mode\ninteractions. Establishing this principle on a solid theoretical footing was our\nprimary goal.\n\nWe believe our findings are more generalizable than they might appear at first\nglance. Many complex problems in fluid dynamics are governed by equations (e.g.,\nNavier-Stokes) that can be seen as a combination of fundamental physical\nprocesses: nonlinear convection (like the Burgers' equation), diffusion, and an\nincompressibility constraint (enforced via a Poisson equation). Our work\ndemonstrates the principle of superiority for each of these core building\nblocks. It is therefore reasonable to hypothesize that the phenomenon will also\nmanifest in more complex systems that combine these elements. Our paper provides\nthe foundational understanding necessary to begin investigating these more\nchallenging scenarios.\n\n> The analysis of nonlinear neural architectures (e.g., UNets, Transformers,\n> FNOs) is relatively scarce. It could be interesting to see more in-depth\n> discussion regarding the relationship between effective receptive field of\n> model and its rollout error / superiority.\n\n\nThis is an excellent suggestion. A deeper dive into how architectural properties\ninfluence superiority would indeed strengthen the paper. We will expand on this\nin the revised manuscript. Below, we present additional interpretation and new\nresults that provide more intuition.\n\n### Performance on Linear vs. Nonlinear Problems\n\nIn the linear advection case (Figure 4), we saw that only architectures with\nstrong, constraining inductive biases (like the local ConvNet) could achieve\nstate-space superiority. Higher-capacity models like UNet and FNO failed to do\nso, as the problem lacked the dynamic complexity to prevent them from simply\noverfitting to the solver's errors.\n\nHowever, the situation changes for the nonlinear Burgers' equation, where the\ndynamics naturally generate a rich spectrum of interacting modes. To provide a\nmore complete picture, we extended our Burgers' experiment to all architectures.\nThe results robustly confirm that superiority is not limited to a single model\nbut is instead closely tied to having a favorable inductive bias. We found that\nall architectures with a spatial or spectral bias (ConvNet, Dilated ResNet, FNO,\nUNet) achieve significant superiority by learning a better representation of the\nemergent shock dynamics. The Dilated ResNet was a standout performer, reaching a\nsuperiority ratio of 0.65. In contrast, the vanilla Transformer, which lacks an\ninherent structural bias, failed to achieve superiority.\n\nBelow are the median superiority ratios ($\\xi^{[t]}$) for the Burgers' equation\nexperiment.\n\n|Network/Time Step|1|2|3|4|5|6|7|8|9|10|\n|:------------------|--:|--:|--:|--:|--:|--:|--:|--:|--:|--:|\n|ConvNet|1|0.97|0.87|0.82|0.8|0.8|0.8|0.81|0.81|0.81|\n|Dilated ResNet|1|0.91|0.79|0.73|0.69|0.66|0.65|0.67|0.66|0.65|\n|FNO|1|1.01|0.91|0.87|0.86|0.89|0.93|1|1.08|1.12|\n|Transformer|1|1.06|1.13|1.24|1.35|1.45|1.52|1.58|1.63|1.68|\n|UNet|1|0.9|0.86|0.88|0.91|0.96|1|1.04|1.07|1.11|\n\nThis suggests that even high-capacity models can achieve superiority in\nphysically rich settings, provided their inductive bias aligns with the problem\nstructure (e.g., local convolutions or spectral filters). It should be noted\nthat modern Transformer variants like SWin resemble a convolution-like inductive\nbias and could thus again show superiority.\n\n### Ablation Study on Model Capacity and Inductive Bias\n\nTo further isolate the effect of architectural properties, we conducted new\nablation studies on both ConvNet and FNO architectures, this time on the linear\nadvection problem. The detailed results can be found in our response to Reviewer\ns59z, but the key findings strongly reinforce the importance of matching the\nmodel's inductive bias to the problem.\n\nFor the ConvNet, we systematically varied its receptive field. We observed that\nsuperiority was maximized when the receptive field size was closely aligned with\nthe physical scale of information propagation (dictated by the CFL number).\nReceptive fields that were too small failed to learn the dynamics, while those\nthat were too large provided excess capacity that allowed the model to begin\nreplicating the solver's non-physical errors, thus weakening the superiority\neffect.\n\nWe found a similar trend for the FNO when varying its number of active Fourier\nmodes. A configuration with just enough modes to represent the core physics\nachieved the best performance, as the FNO's spectral truncation acted as a\npowerful regularizer against the solver's high-frequency artifacts. As we added\nmore modes, the FNO's increased capacity led to it partially overfitting to the\nsolver's error profile, diminishing the magnitude of superiority.\n\nTogether, these studies provide strong evidence that achieving superiority is a\ndirect consequence of the interplay between model capacity and implicit\nregularization. The best results occur when an architecture is expressive enough\nto capture the necessary physics but constrained enough by its inductive bias to\nregularize away the structured flaws of its training data.\n\n> Can the authors provide a general intuition for what kinds of non-linear PDEs\n> are more amenable to neural networks with limited receptive fields (e.g.,\n> feedforward ConvNets), and which ones may require larger context (e.g.,\n> long-range dependencies in turbulent flows or stiff PDEs)?\n\nThis is a very insightful question that gets to the heart of designing effective\nneural emulators (independent of whether one wants to achieve superiority). Our\ngeneral intuition is that the ideal receptive field size is dictated by the\nphysical characteristics of information propagation within the PDE system over a\nsingle time step. In the case of the prototypical advection equation, this is\ndictated by the CFL number which we used to parameterize our experiments.\n\nFor more complicated problems, this depends on the dominant physical effect. If\nthe problem behaves primarily hyperbolic, e.g., compressible Navier-Stokes\nsimulations, then there might be a local region of influence and as long as the\ntime steps are not too large (i.e., information does not travel too far) local\nconvolutional architectures can be used. However, if the dominant effect\nincludes a Poisson-like component (like incompressibility constraints for\nincompressible Navier-Stokes simulations or volume conservation or other global\nconstraints) that could require global information to be propagated. There is\nreason to believe that local emulators can to some extent approximate those\nbehaviors but very likely will not achieve the same performance as global\nemulators. We refer the reviewer to benchmark works that attempt to provide\ncorresponding data points more holistically, e.g., PDEBench, PDEArena or\nAPEBench.\n\n> The authors acknowledge one major limitation in the paper, the solver needs to\n> be sufficiently “coarse” for neural emulator to exhibit superiority.\n\nThis is an important and valid limitation. The requirement that the low-fidelity\nsolver be \"sufficiently coarse\" for emulator superiority to emerge may seem\nrestrictive at first glance. However, we argue that this is not an unreasonable\nassumption in practice. In many real-world scientific and engineering\napplications, practitioners routinely make deliberate simplifications or\napproximations to reduce computational cost or to enable tractable simulations.\nFor example, truncating nonlinear solves is common practice in algorithms like\nPISO for incompressible flow, and reduced-order models are widely used in\nclimate and fluid dynamics. These practical choices often introduce systematic\nerrors or coarsen the solution, creating precisely the kind of setting where a\nneural emulator could, in principle, learn to correct or surpass the baseline\nsolver.\n\nThus, our focus on scenarios with intentionally limited or coarse solvers is\nwell-motivated by practical considerations. While our results suggest that\nneural emulators can sometimes recover or even exceed the accuracy of the\noriginal solver in these settings, we acknowledge that this is not something we\ncan guarantee in general at this stage. The phenomenon is promising, but its\nrobustness and generality across broader classes of problems remain important\ndirections for future research."}, {"review_id": "uun9GPQJ82", "reviewer": "Reviewer_8yEP", "summary": "This paper outlines the idea of emulator superiority, where a neural surrogate trained on low fidelity solver data is able to outperform said low-fidelity solver when evaluated on a high fidelity data. The paper also provides some theoretical analysis showing that a linear emulator in limited settings can have lower error than implicit and explicit numerical schemes in comparison to an analytical solution.", "questions": "1. In the nonlinear PDE experiments, isn't it possible that the reason for neural network superiority on certain timescales is due to the ability of neural networks to smoothly interpolate, while the static low fidelity solver cannot do this?\n\n2. In the case of the advection equation, the narrative says \"Our theoretical analysis reveals several crucial conditions for achieving superiority. First, we must train a low-capacity linear ansatz on a high-capacity (non-analytical) reference method.\" How does this assumption fit with neural network architectures when applied to real data? It is often the case that modern neural networks are high capacity, so doesn't this mean that in many practical scenarios, neural network superiority will not be achieved?\n\n3. In Figure 1, there are labels called \"Amplifying\" and \"Dampening\". Does \"Amplifying\" mean that the scale of the solution is increasing or diverging? And does \"Dampening\" mean the solution scale is decreasing?\n\n4. For the Burger's equation experiment, it is stated that \"for this experiment, we train a vanilla UNet (Ronneberger et al., 2015) on only the initial \"smooth\" frames from this low-fidelity solver.\" Does that mean that the UNet is trained with a limited time-series rollout from the low fidelity solver and then during the test phase, it is made to do a longer rollout than it was trained with?", "limitations": "yes", "rating": 5, "confidence": 4, "soundness": null, "presentation": null, "contribution": null, "strengths": "", "weaknesses": "", "quality": 3, "clarity": 3, "significance": 3, "originality": 3, "strengths_and_weaknesses": "# Strengths:\n1. The authors do a good job of formalizing their problem as well as describing the different types of generalization that are important for neural emulators applied to differential equations.\n2. The paper has useful theoretical results with closed form solutions for linear PDEs using a linear parametric ansatz.\n\n# Weaknesses:\n1. While the theoretical results for linear PDEs are good, the results for the nonlinear neural architectures are mixed, especially in the case of state-space superiority.\n\n2. The practical applicability of this paper is limited because of the problem setup. If a high fidelity simulation data is available, it is unlikely that a researcher would train a neural network with low fidelity data and then evaluate it on high fidelity data. More clarity on when exactly the proposed problem setting is valid is necessary.", "was_revised": true, "final_justification": "In light of some of the clarifications offered by the authors during the rebuttal phase, I believe the question answered by the paper regarding \"emulator superiority\" and corresponding theoretical insights will serve the scientific machine learning community well and hence I will raise my score to `Accept'.", "rebuttal": "We sincerely thank the reviewer for the valuable feedback. We were glad to hear\nthe problem formalization and theoretical results to be considered a key strength\nof the paper. It was our goal to establish a clear, formal foundation for\n\"emulator superiority\" and to ground it in rigorous theoretical analysis, and\nhappy that this was recognized. We appreciate the opportunity to address their\nquestions and concerns.\n\n---\n\n> While the theoretical results for linear PDEs are good, the results for the\n> nonlinear neural architectures are mixed, especially in the case of\n> state-space superiority.\n\nThis is a key observation that we believe reinforces, rather than contradicts,\nour central thesis. The \"mixed\" results are a direct consequence of the\nconditions for superiority we outline, and our aim was to provide a realistic\nassessment: Superiority arises from an interplay\nbetween the emulator's inductive biases and the solver's error characteristics,\nand as such cannot be taken for granted.\n\nOur results in Figure 4c show that high-capacity architectures with global\nreceptive fields (FNO, Transformer) struggle to achieve state-space superiority.\nWe hypothesize this is because their powerful representational capacity allows\nthem to \"overfit\" to the low-fidelity solver, perfectly replicating its\ndynamics—including its systematic errors.\n\nIn contrast, the ConvNet, whose local inductive bias is a poorer match for the\nglobal dynamics of the solver, is forced to learn a more regularized and\nphysically plausible local operator. This \"disadvantage\" in representational\npower becomes an advantage for generalization, leading to superiority. The mixed\nresults are therefore not a weakness of the experiment, but rather an empirical\nvalidation of the conditions required for superiority to emerge.\n\nPlease also note the additional experiments on the Burgers equation we present\nbelow as answer to your second question.\n\n> The practical applicability of this paper is limited because of the problem\n> setup. If a high fidelity simulation data is available, it is unlikely that a\n> researcher would train a neural network with low fidelity data and then\n> evaluate it on high fidelity data. More clarity on when exactly the proposed\n> problem setting is valid is necessary.\n\nWe agree that clarifying the practical context is crucial. Indeed, if enough\nhigh-fidelity (HF) data is available, it is preferable to train on it or to\ncleverly mix it with excessively available low-fidelity (LF) data.\n\nHowever, there are many fields of science and engineering in which high-fidelity\ndata is a scarce \"validation\" resource, e.g., climate science, and turbulence\nmodeling. A researcher may only be able to afford a handful of HF runs to serve\nas a \"gold standard\" for validation. In contrast, they can generate massive\ndatasets using computationally cheap, low-fidelity (LF) solvers. One can imagine\nemulators to be trained on this large LF dataset. Our work demonstrates the\nsurprising outcome that in certain cases such an emulator can learn dynamics that\nare actually more accurate than the training data suggests, potentially \napproaching the quality of an expensive HF reference.\n\nMore broadly, our work highlights that (almost) all numerical solvers are\nimperfect data generators (that being \"high-fidelity\" is not a binary attribute). Our\nframework provides a lens through which to understand how neural networks\ninteract with the inherent, structured errors present in any simulation data,\nwhich is a fundamental question for scientific machine learning.\n\nUltimately, the fairness of benchmarking an emulator depends critically on the\nquality of the reference data: If the reference trajectories used for evaluation\ncontain structured errors, and the neural emulator happens to learn an operator\nthat systematically corrects or improves upon these errors, the benchmarking\nprocess may underestimate the emulator’s true performance. In such cases, the\nemulator could be penalized for being more accurate than the reference, leading\nto an unfair assessment. This is something that is universally applicable and\nmight even appear in settings where the reference data is high-fidelity. \n\nWe will add a paragraph to the introduction to make these practical scenarios\nmore explicit.\n\n> In the nonlinear PDE experiments, isn't it possible that the reason for neural\n> network superiority on certain timescales is due to the ability of neural\n> networks to smoothly interpolate, while the static low fidelity solver cannot\n> do this?\n\nOur interpretation as to why we see autoregressive superiority together with\nstate-space superiority is because the emulators (by virtue of their inductive\nbiases) learn a better generalization to higher frequency modes than their\nlow-fidelity training data generator. Actually, correctly predicting shocks for\nthe Burgers would not require smooth interpolation. \n\n> In the case of the advection equation, the narrative says \"Our theoretical\n> analysis reveals several crucial conditions for achieving superiority. First,\n> we must train a low-capacity linear ansatz on a high-capacity (non-analytical)\n> reference method.\" How does this assumption fit with neural network\n> architectures when applied to real data? It is often the case that modern\n> neural networks are high capacity, so doesn't this mean that in many practical\n> scenarios, neural network superiority will not be achieved?\n\nThe reviewer's intuition aligns with results of the experiment in Figure 4c (for\nwhich nonlinear emulators are tested to generalize from training on only mode 1\nto testing on mode 5). Only the ConvNet with low capacity but strong inductive\nbias shows (autoregressive and) state-space superiority. However, this is only\nfor the advection equation which does not move any energy between the modes. The\nBurgers equation does this and we see that the fairly high-capacity UNet is able\nto achieve superiority on it.\n\nTo provide a more complete picture, we extended our analysis of the nonlinear\nBurgers' equation experiment (Section 4.2) to include all architectures from our\nstudy. The results robustly confirm that the superiority phenomenon is not\nlimited to a single architecture but is instead closely tied to the inductive\nbiases of the models. The key finding is that architectures with a spatial or\nspectral inductive bias (ConvNet, Dilated ResNet, FNO, UNet) all achieve\nsignificant superiority, outperforming the low-fidelity solver they were trained\non. The Dilated ResNet is a standout performer, achieving a remarkable\nsuperiority ratio of 0.65. In stark contrast, the vanilla Transformer, which\nlacks any inherent spatial or spectral bias and treats the input as an\nunstructured set of tokens, fails to achieve superiority and becomes\nprogressively worse than the baseline solver.\n\nBelow are the median superiority ratios ($\\xi^{[t]}$) over 30 random seeds for\neach architecture.\n\n|Network/Time Step|1|2|3|4|5|6|7|8|9|10|\n|:------------------|--:|--:|--:|--:|--:|--:|--:|--:|--:|--:|\n|ConvNet|1|0.97|0.87|0.82|0.8|0.8|0.8|0.81|0.81|0.81|\n|Dilated ResNet|1|0.91|0.79|0.73|0.69|0.66|0.65|0.67|0.66|0.65|\n|FNO|1|1.01|0.91|0.87|0.86|0.89|0.93|1|1.08|1.12|\n|Transformer|1|1.06|1.13|1.24|1.35|1.45|1.52|1.58|1.63|1.68|\n|UNet|1|0.9|0.86|0.88|0.91|0.96|1|1.04|1.07|1.11|\n\nIt should be noted that modern Transformer variants like SWin resemble a\nconvolution-like inductive bias and could thus again show superiority.\n\n>  In Figure 1, there are labels called \"Amplifying\" and \"Dampening\". Does\n>  \"Amplifying\" mean that the scale of the solution is increasing or diverging?\n>  And does \"Dampening\" mean the solution scale is decreasing?\n\nYes, the reviewer's interpretation is correct. These labels describe the effect\nof a numerical scheme on the magnitude of each Fourier mode's complex-valued\ncoefficient. \"Dampening\" means the magnitude of a mode is reduced\n($|\\hat{u}_k^{[t+1]}| < |\\hat{u}_k^{[t]}|$ for mode $k$), which corresponds to\ndiffusion. \"Amplifying\" means the magnitude is increased ($|\\hat{u}_k^{[t+1]}| >\n|\\hat{u}_k^{[t]}|$). For a purely hyperbolic problem like advection, any\nmodification of the magnitude is a numerical artifact. While unphysical,\nnumerical diffusion is often unavoidable and preferred over the opposite. If the\nmagnitude increases, it typically leads to instability (after multiple steps,\nthe solution diverges). This occurs for the explicit scheme beyond its stability\nlimit ($|\\gamma_1| > 1$). We will clarify this directly in the figure caption in\nthe revised manuscript.\n\n> For the Burger's equation experiment, it is stated that \"for this experiment,\n> we train a vanilla UNet ... on only the initial \"smooth\" frames from this\n> low-fidelity solver.\" Does that mean that the UNet is trained with a limited\n> time-series rollout from the low fidelity solver and then during the test\n> phase, it is made to do a longer rollout than it was trained with?\n\nThe reviewer's understanding is correct. The training dataset consists\nonly of one-step transitions $(u^{[0]}, u^{[1]})$ where $u^{[0]}$ is a smooth\ninitial state (low-frequency content). During testing, the emulator is\ninitialized with similar smooth states but must perform a longer autoregressive\nrollout (10 steps). During this longer rollout, the solution develops sharp\nshocks (high-frequency content). The superiority arises because the UNet,\ndespite only being trained on smooth-to-smooth transitions, learns a physical\noperator that correctly propagates these emergent shocks more accurately than\nthe truncated Picard solver it was trained on. \n\n> - Certain figure captions and text refer [...]\n> - Some equations in the appendix are too long [...]\n\nWe thank the reviewer for catching these formatting issues and will improve our\nmanuscript accordingly."}, {"review_id": "XtHbDnnLmk", "reviewer": "Reviewer_jfet", "summary": "Authors investigate when a learned PDE solver can be more accurate than the classical PDE solver used to generate training data.\n\nTheoretical analysis performed for three linear PDEs (convection, diffusion, stationary diffusion) and linear emulators revealed that learned solver can be more accurate if distributions of train and test data do not match.\n\nFor the nonlinear case the improved accuracy of learned solver was observed for nonlinear emulators (FNO, Dilated ResNet, Transformer, etc) on convection and Burgers equations for the case of identical distributions for train and test data.", "questions": "1. **Assumptions in the proof of theoretical results**\n\n   Authors assumed that the training data come from a very specific distribution that covers only a single mode. In this case it is clear (not explained explicitly) that the learned model will depend on the train set only from the mode number and both amplitude and phase are irrelevant. I have several questions on that:\n\n   a. Can the authors please supply missing proof that both amplitude and phase are irrelevant? I would suggest showing that equation (8) can be achieved when the usual loss function is considered, e.g., mean (with respect to distribution on phase and amplitude) $L_2$ error.\n\n   b. The whole formulation seems rather artificial with a weak dependence on training data. Can the authors please defend their choice of learning problem? Why is this simple scheme sufficient and interesting for operator learning?\n\n   c. Since authors work with linear models, it seems to be possible to obtain closed-form solutions for arbitrary data. Can the authors comment on the possible extension of their results on more realistic data distributions (e.g., more than one mode)?\n   \n2. **Theoretical results are trivial from the standpoint of backward error analysis**\n\n   Solution for PDE can be considered as an evaluation problem $y = \\phi(x)$, where $x$ is input data, $y$ is solution to PDE and $\\phi$ is an unknown map that brings input data to the solution. The approximate method that solves PDE is some other map $\\hat{\\phi}(x)$ that approximates $\\phi(x)$ sufficiently well. Standard forward error used by authors is $\\left\\|\\phi(x)-\\hat{\\phi}(x)\\right\\|$, where $\\left\\|\\cdot\\right\\|$ is a suitable norm.\n   \n   Backward error is based on an entirely different idea. Consider perturbation of input data $\\Delta x$ such that $\\hat{\\phi}(x) = \\phi(x + \\Delta x)$. In case this perturbation is possible to find $\\hat{\\phi}(x)$ can be considered as an *exact* solver for *wrong* input data. For solver $\\hat{\\phi}$ backward error is thus $\\left\\|\\Delta x\\right\\|$ (here $\\left\\|\\cdot\\right\\|$ is potentially different norm) and is related to forward error by condition number of the problem.\n   \n   If both $\\phi(x)$ and $\\hat{\\phi}(x)$ are continuous one may expect that one can find input $\\widetilde{x}$ near $x + \\Delta x$ such that $\\left\\|\\phi(\\widetilde{x}) - \\hat{\\phi}(\\widetilde{x})\\right\\| \\leq \\epsilon$ for small $\\epsilon$. Such inputs can be found by solving $\\arg\\min_{x}\\left\\|\\phi(x) - \\hat{\\phi}(x)\\right\\|$ (at least locally near $x$).\n   \n   Authors considered two different solution methods for the same problem $\\phi_1(x)$, $\\phi_2(x)$. These methods in general have different backward errors $\\Delta x_1(x)$ and $\\Delta x_2(x)$. Because of that these methods also have different solutions for the optimisation problem above $x_1 = \\arg\\min_{x}\\left\\|\\phi(x) - \\phi_1(x)\\right\\|\\neq x_2 = \\arg\\min_{x}\\left\\|\\phi(x) - \\phi_2(x)\\right\\|$.\n   \n   These observations do not look too interesting or surprising. And the main problem is we deliberately try to find inputs such that one method is superior to the other one. Several questions are in order:\n\n   a. Can the authors please comment on the considerations on backward error analysis?\n\n   b. On the higher level authors simply note that different methods are good on different inputs. This observation is well known in the literature on numerical methods. For example, trapezoidal rules have exponential convergence for periodic functions and are expected to give very small errors in this case, but for non-periodic functions the error will increase. Given that it is easy to find functions for which trapezoidal rule is superior to, say, Simpson's rule and vice versa. Can the authors explain what made their result conceptually different?\n\n   c. There is a clear parallel between adversarial examples for DL classification and inputs that lead to superiority. Did authors try to consider this parallel? What authors think in general on this relation?\n   \n3. **Empirical results for Burgers equation are not explained by theory**\n\n   In my view the results presented for Burgers equation is the most surprising one. It is surprising because superiority shows up in the same solution curve, still out of distribution but not completely unrelated to the training data.\n\n   a. It does not seem to me that theoretical results explain this behaviour. Can authors comment on whether a linear case is related to nonlinear one?\n\n   b. Can this result be replicated for other equations (e.g., KS, KdV) and architectures (e.g., FNO, ResNet)?\n\n   c. Do these results survive when one considers more diverse training sets?\n\nIn conclusion, I find questions raised by the author interesting and potentially fruitful for further study. Unfortunately, I do not find answers that authors provide to be useful or particularly interesting. Being said that I remain open to change my opinion after the rebuttal. I encourage authors to defend their position, supply more results and continue research on the topic.", "limitations": "yes", "rating": 4, "confidence": 4, "soundness": null, "presentation": null, "contribution": null, "strengths": "", "weaknesses": "", "quality": 2, "clarity": 3, "significance": 2, "originality": 4, "strengths_and_weaknesses": "**Strengths:**\n1. Authors ask genuinely interesting question\n2. Both theoretical and empirical considerations are available\n\n**Weaknesses:**\n1. Both theoretical and empirical results are not convincing\n2. The setup authors consider is far from establish practices\n\nI provide more details in the section below.", "was_revised": true, "final_justification": "To summarise, my opinion does not change after reading other reviews and rebuttal by authors:\n\n1. The question authors ask is genuinely interesting\n2. Available theoretical results are not well-developed\n3. Numerical results are not convincing", "rebuttal": "We sincerely thank the reviewer for their deep and critical engagement with our\nwork. \n\n---\n\n> - Both theoretical and empirical results are not convincing  \n> - The setup authors consider is far from established practices\n\nWe understand the reviewer's concerns and will address the specific points that\nled to this conclusion in the detailed sections below. We believe our setup,\nwhile simplified for analytical tractability, is designed to isolate a\nfundamental phenomenon relevant to the established practice of using imperfect\nsimulation data to train emulators, a common scenario in resource-constrained\nscientific and engineering domains.\n\n> Can the authors please supply missing proof that both amplitude and phase are\n> irrelevant? [...]\n\nThis is an excellent point. The proof was omitted for brevity, but we agree it\nis essential for the rigor of our argument. The key insight is that for a linear\nconvolutional ansatz and a Mean-Squared Error (MSE) loss, the optimization\nproblem diagonalizes in Fourier space. *If only a single mode is present*, the\nexpectation over the data becomes a constant factor (i.e., the distribution over\nphase and amplitude is irrelevant as long as it is nonzero) and the minimizer is\nfound by the first-order condition, which leads directly to Equation (8). We\nwill add the derivation to the appendix of the revised pdf and are happy to\nprovide it as a discussion reply if the reviewer is interested.\n\n> The whole formulation seems rather artificial with a weak dependence on\n> training data. Can the authors please defend their choice of learning problem?\n> Why is this simple scheme sufficient and interesting for operator learning?\n\nWe acknowledge that the single-mode training setup is a simplification. This\nchoice was deliberate and serves two purposes. This controlled setting is what\nallows us to derive a data-independent closed-form analytical solution for the\ntrained emulator. This provides a clear, interpretable, and undeniable\ndemonstration of the superiority principle without confounding numerical\nfactors. Moreover, probing a system's behavior with single Fourier modes is a\nfoundational technique in linear systems theory and numerical analysis. It\nallows for a precise characterization of a system's response across different\nfrequencies.\n\nOur \"simple scheme\" is therefore sufficient and interesting because it allows us\nto *isolate and prove* the existence of a counter-intuitive phenomenon. This\nprovides the foundational insight needed to understand why superiority might\noccur in more complex, realistic operator learning scenarios where analytical\ntreatment is impossible.\n\n> Since authors work with linear models, it seems to be possible to obtain\n> closed-form solutions for arbitrary data. Can the authors comment on the\n> possible extension of their results on more realistic data distributions\n> (e.g., more than one mode)?\n\nIf the training data contains multiple modes, the closed-form solution for the\nemulator parameters is no longer a simple equation for a single mode. Instead,\nit becomes a linear least-squares weighted by the concrete distribution of\ncomplex-values coefficients. While an analytical solution still exists for this\nlinear regression problem, the analysis becomes more complex because of this\ndata dependency.\n\nTo illustrate this, we redid the theoretical experiments numerically (training\non a combination of modes and only testing on mode 5). We found that there is\nstill superiority when more modes are present in the training initial condition\nas long as they are lower than what is being tested on (confirming again what we\ncalled \"forward superiority\" in line 202 of the original manuscript).\n\n|Modes in Train IC Dist|Superiority $\\xi^{[1]}$|\n|:---|:---|\n|1|0.57|\n|1,2|0.61|\n|1,2,3,4|0.72|\n|1,2,3,4,5|0.81|\n|1,2,3,4,5,6|0.87|\n|1,2,3,4,5,6,7|1.0|\n|1-10|1.1|\n|1-20|2.3|\n|4,6|1.1|\n|3,4|0.79|\n\n> Can the authors please comment on the considerations on backward error\n> analysis?\n\nWe thank the reviewer for bringing up this insightful perspective. We agree that\nfrom a backward error analysis standpoint, our low-fidelity solver and our\nlearned emulator can be viewed as exact solvers for two slightly different,\nperturbed problems. However, this perspective, while valid, misses the central\nand most novel aspect of our work: the learning dynamic itself.\n\nOur contribution is not merely to show that two different numerical methods\nexist and have different error properties. Our key finding is that we can learn\na new method (the emulator) by only observing an existing, inferior method (the\nlow-fidelity solver), and this new method can provably be closer to a\nhigh-fidelity reference. Standard backward error analysis compares pre-existing\nmethods; it does not account for the process of creating a new, superior method\nfrom an imperfect one, which is the core of our insights.\n\n>  ...authors simply note that different methods are good on different inputs.\nThis observation is well known... Given that it is easy to find functions for\nwhich trapezoidal rule is superior to, say, Simpson's rule and vice versa. Can\nthe authors explain what made their result conceptually different?\n\nThe conceptual difference lies in the origin of the superior method. The analogy\nof choosing between the Trapezoidal and Simpson's rules involves selecting from\na set of pre-existing, analytically derived methods. A user knowingly picks the\nright tool for the job.\n\nIn our work, the superior method (the emulator) is not analytically pre-defined.\nIt is discovered through a data-driven learning process that has access only to\ndata from an inferior method. The novelty is in demonstrating that this learning\nprocess can implicitly correct the known, structured flaws of its training data\nsource, a mechanism different from selecting an existing quadrature rule.\n\n> There is a clear parallel between adversarial examples for DL classification\n> and inputs that lead to superiority. Did authors try to consider this\n> parallel? What authors think in general on this relation?\n\nBoth phenomena reveal the \"fragility\" and specificity of what a neural network\nlearns. An adversarial example exploits a model's learned decision boundary to\ncause a misclassification, often via a perceptually minor perturbation. This is\ntypically viewed as a failure mode. A superiority-inducing input, as we've\nframed it, exploits the difference between the emulator's learned dynamics and\nthe solver's flawed dynamics. However, we see this as a constructive outcome.\nThe emulator is \"fragile\" to the specific errors of the solver, and instead of\nreplicating them, it learns a more physically regularized behavior. One could\nview our work as identifying a regime where a model's sensitivity to\nout-of-distribution data does **not** lead to a nonsensical output, but to a\nphysically more plausible one.\n\n> It does not seem to me that theoretical results explain this behaviour. Can\n> authors comment on whether a linear case is related to nonlinear one?\n\nThe theoretical results provide the conceptual blueprint. The key insight from\nour linear analysis is that an emulator with a favorable inductive bias can\nachieve superior spectral generalization; it can learn to handle high-frequency\nmodes more accurately even when only trained on low-frequency data.\n\nThis is also what we observe in the Burgers' experiment. The low-fidelity solver\n(truncated Picard iteration) makes significant errors in propagating the\nhigh-frequency components that constitute the sharp shock front. The UNet is\ntrained only on the first two frames, which are dominated by low-frequency\ncontent. To succeed at the multi-step rollout, the UNet must learn an operator\nthat correctly handles the emergent high-frequency modes of the shock. We\nbelieve that its ability to do so, outperforming its training data, is a direct\nnonlinear analogue of the spectral generalization proven in our linear theory.\n\n> Can this result be replicated for other equations (e.g., KS, KdV) and\n> architectures (e.g., FNO, ResNet)?\n\n**Architectures:** Yes, the superiority phenomenon on the Burgers' equation is\nnot unique to the UNet. As detailed in our response to Reviewer K6sC, we\nrepeated this experiment for all architectures from Section 4.1. The key finding\nfrom this new experiment is that architectures with favorable inductive\nbiases—namely ConvNet, Dilated ResNet, FNO, and the UNet—all achieve significant\nperiods of superiority (superiority ratio < 1). This demonstrates the generality\nof the effect. Conversely, the vanilla Transformer, which lacks a strong spatial\nor spectral prior, fails to achieve superiority. This strongly suggests that a\nsuitable inductive bias is a more critical ingredient for this phenomenon than\nraw model capacity. We hypothesize that modern Transformer variants\nincorporating convolutional-like structures (e.g., Swin Transformers) would\nlikely exhibit superiority as well.\n\n**Equations:** Our study was intentionally designed to analyze the fundamental\n\"building blocks\" of more complex physical systems. Many sophisticated problems,\nsuch as the Navier-Stokes equations in fluid dynamics, are governed by a\ncombination of simpler physical processes. Our paper provides the foundational\nunderstanding, and such cases represent very interesting directions for future\nwork.\n\n**Regarding the specific examples:** We hypothesize that the Korteweg-de Vries\n(KdV) equation, which develops rich, non-chaotic spectral content (solitons), is\na strong candidate for exhibiting similar superiority as in the Burgers\n(primarily adding a third-order derivative). Chaotic systems like the KS\nequation present more complex challenges where regions with superiority may be\nharder to identify.\n\n> Do these results survive when one considers more diverse training sets?\n\nThis is a key question for future work; the new results above give first\nindicators for generalization to diverse training sets."}, {"review_id": "xMUEJkU70I", "reviewer": "Reviewer_GMRy", "summary": "This paper introduces and investigates \"emulator superiority\". An interesting phenomenon where neural emulators trained solely on low-fidelity numerical solver data can outperform these solvers when evaluated against higher-fidelity references. Through theoretical analysis on linear PDEs, the authors demonstrate that certain neural architectures' inductive biases enable better error accumulation or generalization than their training solvers. Empirical experiments extend this concept to nonlinear PDEs, showing superiority across different architectures.", "questions": "- The paper notes that superiority may disappear when training data is already high-fidelity. Can we quantify a fidelity threshold below which emulators have room to improve?\n- If the core motivation of the work is to advocate for low-fidelity simulation-based emulators, it would be beneficial to include a discussion on the trade-offs in computational cost. Currently, there is limited analysis of the cost-benefit balance when using neural emulators, particularly considering that solving high-fidelity PDEs with neural models can involve significant training and inference costs, which may offset the claimed advantages.\n- How robust are the closed-form superiority proofs for linear PDEs to perturbations in discretization, boundary conditions, or higher-dimensional settings? Are these proofs extensible beyond 1D periodic domains and synthetic solvers?\n- Were baseline solvers used in training (e.g., Picard-truncated Burgers solver) intentionally made poor to create room for superiority? Could the superiority vanish if the low-fidelity solver were just marginally improved? \n- Can you try “non-toy” PDE systems (e.g., Navier-Stokes)?\n- To me, some of the training conditions (e.g., ICs from a single Fourier mode) are overly restrictive. How would the results change if training data sampled from more physically diverse regimes?", "limitations": "yes", "rating": 4, "confidence": 3, "soundness": null, "presentation": null, "contribution": null, "strengths": "", "weaknesses": "", "quality": 2, "clarity": 3, "significance": 2, "originality": 3, "strengths_and_weaknesses": "Strengths:\n- One of the strong aspects of this paper is its theoretical analysis. The authors provide clear and detailed proofs in the Fourier domain demonstrating conditions under which the \"emulator superiority\" can arise for linear PDEs.\n- The core concept of emulator superiority is interesting.\n\nWeaknesses:\n- Despite robust theoretical treatment for linear PDEs, the paper’s empirical experiments are somewhat limited to idealized, canonical PDE examples. While appropriate for proof-of-concept, the generalizability to real-world multi-physics simulations (e.g., climate, combustion, plasma physics) remains unexplored.\n- The paper sticks to one-step supervised training. It acknowledges but does not explore other paradigms like unrolled training or physics-informed losses, which might enhance or inhibit superiority. This makes me wonder to what extent is “superiority” an artifact of the chosen loss metric, test regime, or evaluation horizon? Would the observed superiority still exist if evaluated on alternative physics-informed metrics?\n- The superiority proofs rely on mode-isolated training and evaluation. For practical PDE problems, this may not be the case. Is the assumption of single-mode initial conditions meaningful in real scientific scenarios?", "was_revised": false, "final_justification": "", "rebuttal": "We are deeply grateful to the reviewer for their thorough and critical review.\n\n---\n\n> Despite robust theoretical treatment for linear PDEs, the paper’s empirical\n> experiments are somewhat limited to idealized, canonical PDE examples. \n\nWe concur with the reviewer that extending this work to complex multi-physics\nsimulations is a critical direction for future research. Our focus on canonical\nPDEs was a strategic choice, allowing us to isolate and rigorously analyze the\nfoundational principles of emulator superiority. We argue that this approach has\nbroader implications because many complex systems, like those governed by the\nNavier-Stokes equations, are fundamentally composed of the very dynamics we\nstudy: advection, diffusion, and Poisson-like solves. By demonstrating that\nsuperiority emerges within these core building blocks, we provide strong\nevidence that the phenomenon is a general principle of learning from numerical\ndata, not an artifact of simple systems. Understanding these dynamics in\nisolation is a necessary foundation for investigating more complex scenarios\nwhere multiple error sources are coupled and confounded. We will highlight this\nin the revised version.\n\n> [...] to what extent is “superiority” an artifact of the chosen\n> loss metric, test regime, or evaluation horizon? Would the observed\n> superiority still exist if evaluated on alternative physics-informed metrics?\n\nWe agree that the training paradigm is a critical factor. Our linear theoretical\ninvestigations hold equally for unrolled training and physics-informed training\n(if in the form of discrete residuals over two consecutive time steps). If the\nreviewer is interested, we'd be happy to provide the full proof, and we will add\nit to the appendix of the revised version.\n\nFor the nonlinear problems, we also investigated with unrolled training and\nfound that it oftentimes reduces the effect of superiority because the emulator\nis expected to replicate the error accumulation of the solver. This is desirable\nif the reference solver is good, but not if it is bad. Since there are\nadditional axes to unrolled training (like its length, weighting, additional\nreference frames etc.), we decided to focus on the one-step training paradigm\nfor this paper.\n\nThe reviewer also raises an interesting point about training and\nevaluating on physics-informed (PI) metrics. While the equivalence of one-step\nsupervised training and discrete residua cannot directly be shown for more\ncomplicated setups, we still believe that the discrete residuum consistent with\na coarse solver will give rise to the same superiority effect. Note that our\nformal definition of Equation (5) requires solvers and does not directly\ntranslate to PI metrics.\n\n> The superiority proofs rely on mode-isolated training and evaluation. For\n> practical PDE problems, this may not be the case. Is the assumption of\n> single-mode initial conditions meaningful in real scientific scenarios?\n\nFirst, we wish to clarify that while our theoretical proofs utilize single-mode\nanalysis for analytical tractability, a significant portion of our empirical\nresults already incorporate multi-mode dynamics. The experiment in Figure 4c\nexplicitly tests on five active modes, while the nonlinear dynamics in the\nBurgers' experiment (Figure 5) generate a rich, interactive spectrum from smooth\ninitial data. The theoretical focus on a single mode is a standard diagnostic\ntool from numerical analysis, which allowed us to isolate the core mechanism of\nsuperiority with a degree of rigor and clarity that is often unattainable in\nmore complex numerical settings.\n\nTo further bridge this gap and directly address the reviewer's concern, we have\nconducted a new numerical experiment based on our linear theoretical setup. As\ndetailed further in our response to Reviewer jfet, we trained a linear emulator\non initial conditions containing various combinations of Fourier modes and\ntested its state-space superiority on an unseen higher mode. This study confirms\nthat \"forward superiority\" is a robust phenomenon: as long as the emulator is\ntrained on modes with lower frequencies than the test mode, it consistently\nachieves superiority ( $\\xi<1$). The results also reveal an intuitive pattern\nwhere the magnitude of this effect diminishes as the training data becomes more\nspectrally diverse. By seeing the solver's error profile across a wider range of\nmodes, the emulator begins to replicate those systematic flaws more closely.\n\n> The paper notes that superiority may disappear when training data is already\n> high-fidelity. Can we quantify a fidelity threshold below which emulators have\n> room to improve?\n\nThis is a key question. While a universal, quantifiable threshold is likely\nproblem-dependent, we can demonstrate the relationship between solver fidelity\nand the potential for superiority. We have conducted an ablation study for the\nnumber of Picard iterations in the Burgers example of section 4.2. We again used\nthe UNet and display the superiority ratio $\\xi^{[t]}$ rollout when trained and\nbaselined against a Picard solver truncated after a certain number of\niterations. A single iteration represents our original low-fidelity setup, while\na higher number of iterations produces a more accurate, higher-fidelity solver.\n\n|Picard Iter/Time Step|1|2|3|4|5|6|7|8|\n|---:|---:|---:|---:|---:|---:|---:|---:|---:|\n|1 (Lowest Fidelity)|1|0.91|**0.85**|0.88|0.92|0.96|1|1.05|\n|2|1|1|0.97|0.98|1.01|1.09|1.16|1.25|\n|3|1|1.05|1.05|1.11|1.19|1.27|1.37|1.47|\n|4|1|1.36|1.34|1.46|1.53|1.58|1.67|1.84|\n\nThese results clearly demonstrate our central thesis: the potential for\nsuperiority is directly related to the magnitude of the structured error in the\ntraining data. When the training solver is very coarse, there is significant\nroom for the emulator to learn a more regularized operator, resulting in a\nstrong superiority effect ($\\xi$ ratio dropping to 0.85). As the fidelity of the\ntraining solver increases, the structured error in the training data decreases,\nand the potential for superiority vanishes. When trained on a reasonably\nconverged solver, the emulator simply learns to replicate its already accurate\nbehavior, and no superiority is observed. \n\n> If the core motivation of the work is to advocate for low-fidelity simulation-based emulators, it would be beneficial to include a discussion on the trade-offs in computational cost.\n\nWe thank the reviewer for this point, as it allows us to clarify our motivation.\nOur work addresses the common scientific workflow where emulators are trained on\ncomputationally feasible but imperfect solver data. Rather than a simple\ncost-benefit analysis, our core finding is twofold. First, we show that an\nemulator can learn a more regularized operator and surpass the fidelity of its\ntraining data, effectively bridging the gap towards a more accurate reference.\nSecond, this has critical implications for benchmarking, as a superior emulator\nmay be unfairly penalized by a flawed reference, a paradox that urges more\ncareful evaluation practices across the field. We will revise the introduction\nto make this practical context and these fundamental implications more explicit.\n\n> How robust are the closed-form superiority proofs for linear PDEs? \n> ... Are they extensible?\n\n- Higher Dimensions: Our proofs extend directly to n-dimensional periodic\n  domains (on equidistant rectangular grids), as the Fourier basis remains\n  separable.  \n- Boundary Conditions: For non-periodic boundaries, the Fourier basis is no\n  longer applicable. One would need to use other spectral bases (e.g.,\n  Chebyshev), which do not diagonalize the convolution operator in the same\n  simple way, making a closed-form proof significantly more complex.\n- Discretization: Our analysis is at a fixed discretization. Changing it changes\n  the parameters and thus the magnitude of superiority, but the underlying\n  mechanism remains.\n\nIn summary, we believe our theoretical analysis is already presented in the most\ngeneral setting that still permits a rigorous, closed-form proof. Beyond these\nconditions, numerical experiments are more appropriate.\n\n> Were baseline solvers used in training (e.g., Picard-truncated Burgers solver)\n> intentionally made poor to create room for superiority? Could the superiority\n> vanish if the low-fidelity solver were just marginally improved?\n\nThe solvers were chosen because they have clear, well-understood error\ncharacteristics (e.g., the known numerical diffusion of an implicit scheme or\nthe incorrect shock movement of truncated nonlinear solves). This is essential\nfor transparent analysis and not unrealistic; for example, truncating nonlinear\nsolves is common practice in algorithms like PISO for incompressible flow (see\nalso (Turek, 1999) for a further discussion). As noted in our response to your\nprevious question, superiority does indeed diminish as the solver is improved.\nThis is an expected and central part of our thesis: the potential for\nimprovement is directly related to the magnitude of the structured error in the\ntraining data.\n\n> Can you try “non-toy” PDE systems (e.g., Navier-Stokes)?\n\nThank you for this suggestion. Extending our analysis to more complex PDE\nsystems such as Navier-Stokes is an exciting direction for future work, but unfortunately beyond the scope of the rebuttal.\n\n> How would the results change if training data sampled\n> from more physically diverse regimes?\n\nThis is an important point we addressed in our reply to a previous point with\nthe experiment of mixed mode in the training dataset. We re-emphasize that our\nkey nonlinear experiments already use multi-mode data for testing and/or involve\nspectrally rich dynamics. A difference between the training and testing\ndistributions is a necessary precondition for state-space superiority, which is\none of the phenomena we aim to characterize."}]

Author rebuttal:

Meta-review: This paper formalizes and studies an interesting phenomenon where neural operators trained on data produced from low-fidelity solver (e.g., only apply 1 iteration of spatial Newton solver in each time step) can sometimes surpass the solver itself when compared to high-fidelity references (analytical for example). All reviewers found the experiments interesting, so did I. I would like to mention that there is some new development in interpreting operator learning as a statistical learning method (i.e., discretization invariance is really an implication of in-distribution inference), see e.g., Reinhardt et al. arXiv:2412.17582, so the coarse data still may correctly capture the correct distribution.

Here are my revision suggestions in addition to the reviewers': CFL condition ($\gamma$, or Courant number) is usually given as positive, please refer to any FDM textbook and fix this. Note for heat eq-like spatiotemporal problems (3.2), bounds on $\gamma_2$ are refered to as CFL condition as well (for explicit marching schemes). Section 3 title reads "Proof", however, after reading it I felt it was more like "A Case for Supriority". The letter $m$ and $M$ need to be more clearly defined in Section 3, I have to refer to the Appendix to find their meanings in Figure 1.
Decision: Accept (poster)