Paper deep dive
HERO: History-Enriched Rollout Training for Long-Horizon Autoregressive Neural Operators
Jiaquan Zhang, Shuxu Chen, Haifan Meng, Yi Lu, Zhihan Lyu, Fan Mo, Wei Dong, Yang Yang, Chaoning Zhang
Intelligence
Status: succeeded | Model: Gemma-4-26B-A4B | Prompt: intel-v1 | Confidence: 88%
Last extracted: 8/3/2026, 2:40:05 AM
Summary
The paper introduces HERO (History-Enriched Rollout Training), a framework for training autoregressive neural operators to solve time-dependent partial differential equations (PDEs). It addresses the issue of error accumulation in long-horizon rollouts by supplementing standard absolute trajectory supervision with relative supervision derived from the model's optimization history. HERO dynamically selects a 'failure trajectory' from a lagged operator, the current model, or a perturbed input to serve as a reference baseline. This reference is used in a margin-based objective to reweight the ground-truth rollout gradient, improving long-horizon accuracy, stability, and robustness without increasing inference-time cost.
Entities (7)
Relation Signals (6)
HERO → improves → long-horizon accuracy
confidence 95% · Experiments on nine PDE benchmarks ... show that HERO consistently improves long-horizon accuracy, stable rollout length, and out-of-distribution robustness
HERO → supplements → absolute trajectory supervision
confidence 90% · HERO, which augments conventional absolute trajectory supervision with relative supervision derived from the model's optimization history.
HERO → uses → Lagged Operator
confidence 90% · HERO constructs up to three candidate rollouts: ... a rollout from the lagged operator
Lagged Operator → generates → failure trajectory
confidence 85% · selects the strongest failure trajectory as reference ... from a periodically refreshed lagged operator
Margin-based Objective → induces → gradient reweighting
confidence 85% · inducing a bounded, sample-dependent reweighting of the ground-truth rollout gradient rather than an independent gradient direction
FNO-PF → isvariantof → FNO
confidence 80% · we compare the standard Fourier neural operator (FNO) with its push-forward variant (FNO-PF)
Cypher Suggestions (0)
No Cypher suggestions yet.
Abstract
Abstract:Neural operators provide fast surrogates for time-dependent partial differential equations (PDEs) by applying a learned evolution operator recursively to its own predictions, but this autoregressive rollout feeds every prediction error back as input, so local errors accumulate. Existing rollout-training strategies reduce the mismatch between training inputs and self-generated states, yet their supervision still measures only the absolute discrepancy from the ground-truth trajectory. Such supervision is therefore uninformative about whether the operator has overcome the long-horizon failure behaviors it exhibited earlier during optimization. We propose history-enriched rollout training (HERO), which augments conventional absolute trajectory supervision with relative supervision derived from the model's optimization history. HERO ranks detached candidate rollouts from a periodically refreshed lagged operator, the current model, and a perturbed input by rollout error, spectral discrepancy, energy drift, and error growth, and selects the strongest failure trajectory as reference. This reference enters a margin-based objective as a fixed comparison baseline, inducing a bounded, sample-dependent reweighting of the ground-truth rollout gradient rather than an independent gradient direction, which we further analyze theoretically. Experiments on nine PDE benchmarks with spectral and attention-based backbones show that HERO consistently improves long-horizon accuracy, stable rollout length, and out-of-distribution robustness at no inference-time cost. These results indicate that history-enriched relative supervision is effective for stabilizing long-horizon autoregressive prediction.
Tags
Links
- Source: https://arxiv.org/abs/2607.29135v1
- Canonical: https://arxiv.org/abs/2607.29135v1
Trouble viewing inline? Open PDF directly →
Full Text
131,054 characters extracted from source content.
Expand or collapse full text
HERO: History-Enriched Rollout Training for Long-Horizon Autoregressive Neural Operators Jiaquan Zhang School of Computer Science and Engineering, UESTCChengduChina , Shuxu Chen Electronics and Information Convergence Engineering, KHUYongin-siKorea , Haifan Meng School of Computer Science and Engineering, UESTCChengduChina , Yi Lu Department of Mathematical Sciences, UOLLiverpoolEngland , Zhihan Lyu School of Computer Science and Technology, XDUShanxiChina , Fan Mo School of Mechanical and Electrical Engineering, UESTCChengduChina , Wei Dong College of Computer and Information Engineering, XAUATXi’anChina , Yang Yang School of Computer Science and Engineering, UESTCChengduChina and Chaoning Zhang School of Computer Science and Engineering, UESTCChengduChina (2018) Abstract. Neural operators provide fast surrogates for time-dependent partial differential equations (PDEs) by applying a learned evolution operator recursively to its own predictions, but this autoregressive rollout feeds every prediction error back as input, so local errors accumulate. Existing rollout-training strategies reduce the mismatch between training inputs and self-generated states, yet their supervision still measures only the absolute discrepancy from the ground-truth trajectory. Such supervision is therefore uninformative about whether the operator has overcome the long-horizon failure behaviors it exhibited earlier during optimization. We propose history-enriched rollout training (HERO), which augments conventional absolute trajectory supervision with relative supervision derived from the model’s optimization history. HERO ranks detached candidate rollouts from a periodically refreshed lagged operator, the current model, and a perturbed input by rollout error, spectral discrepancy, energy drift, and error growth, and selects the strongest failure trajectory as reference. This reference enters a margin-based objective as a fixed comparison baseline, inducing a bounded, sample-dependent reweighting of the ground-truth rollout gradient rather than an independent gradient direction, which we further analyze theoretically. Experiments on nine PDE benchmarks with spectral and attention-based backbones show that HERO consistently improves long-horizon accuracy, stable rollout length, and out-of-distribution robustness at no inference-time cost. These results indicate that history-enriched relative supervision is effective for stabilizing long-horizon autoregressive prediction. †copyright: acmlicensed†journalyear: 2018†doi: X.X†conference: Make sure to enter the correct conference title from your rights confirmation email; June 03–05, 2018; Woodstock, NY†isbn: 978-1-4503-X-X/2018/06†ccs: Computing methodologies Neural networks 1. Introduction During autoregressive inference, each prediction is fed back as input, repeatedly propagating local errors and shifting the rollout away from the training distribution (Li et al., 2023a; McCabe et al., 2023). Nonlinear, multiscale, and repeated spectral dynamics can amplify these deviations, causing spectral distortion and eventual rollout instability (McCabe et al., 2023; Lippe et al., 2023; Li et al., 2026; Jiang et al., 2025). Low one-step error therefore does not guarantee accurate and stable long-horizon rollouts (McCabe et al., 2023; Lippe et al., 2023). For chaotic systems, long-term statistical fidelity must be evaluated separately from pointwise one-step accuracy (Li et al., 2022; Schiff et al., 2024). Existing approaches mitigate this mismatch by exposing the operator to self-generated states during training, including multi-step rollout training, push-forward strategies, differentiable-solver strategies, and recurrent operator training (Brandstetter et al., 2022; Um et al., 2020; Ye et al., 2025). These rollout-training strategies reduce the mismatch between training and inference inputs by modifying the states encountered during training, while their objectives remain grounded in discrepancies between predicted and reference trajectories rather than an explicitly history-enriched failure reference (Brandstetter et al., 2022; Um et al., 2020; Ye et al., 2025). Such supervision measures how far the current rollout deviates from the target evolution, but does not explicitly reveal whether its long-horizon error growth or previously observed failure patterns persist (McCabe et al., 2023; List et al., 2025; Lippe et al., 2023). To illustrate this limitation, we compare the standard Fourier neural operator (FNO) with its push-forward variant (FNO-PF), which propagates model-generated states into subsequent steps and updates the operator on the resulting rollout errors. As shown in Figure 1, FNO-PF delays the onset of degradation relative to FNO but retains a similar long-horizon growth trend. Exposure to self-generated states therefore improves long-horizon prediction, yet absolute supervision still cannot indicate whether the operator has overcome the failure behaviors it produced earlier during optimization. The model’s optimization history provides a natural source of comparative information for addressing this supervision gap. We treat parameter states saved at earlier optimization stages as historical checkpoints; each checkpoint defines a complete autoregressive operator whose long-horizon rollout is evaluated from the same initial conditions. Comparing these rollouts against the final operator (Figure 2) shows that a checkpoint’s relative quality depends on the rollout horizon and does not improve monotonically with training: some checkpoints stay competitive over the initial steps but fall behind at longer horizons, whereas others differ in their transition times and error-growth patterns. Historical checkpoints are therefore not uniformly weaker copies of the final operator; each retains a distinct long-horizon failure behavior. Moreover, as the current operator improves, a fixed historical trajectory becomes too easy to outperform, reducing the contribution of the relative objective. These observations motivate adapting the reference to both the rollout behavior and the current optimization stage. We propose a history-enriched rollout training framework (HERO) that supplements conventional ground-truth trajectory regression with relative long-horizon supervision. Instead of relying on a fixed historical checkpoint, HERO dynamically constructs up to three detached candidate trajectories from a periodically refreshed lagged operator, the current rollout, and an optional locally perturbed input. It evaluates their rollout error, spectral discrepancy, energy deviation, and error-growth behavior, and selects an informative failure trajectory as the comparison reference. The current operator is then optimized both to match the ground-truth evolution and to outperform the selected reference through a margin-based objective. The selected reference is used only as a detached comparison baseline and introduces no independent gradient direction. Instead, the relative objective induces a bounded, sample-dependent reweighting of the current ground-truth rollout gradient. We characterize this mechanism through a covariance decomposition of the aggregated update direction and establish conditional finite-step reachability of a prescribed long-horizon risk region under explicit regularity and alignment conditions. HERO is architecture agnostic, and all reference-construction components are used only during training, leaving the inference-time architecture and cost unchanged. The main contributions are summarized as follows: (i) We propose HERO, an architecture-agnostic history-enriched rollout training framework that converts dynamically selected long-horizon failure behaviors into relative supervision while retaining conventional ground-truth trajectory regression. (i) We characterize the optimization mechanism induced by stop-gradient reference trajectories. The reference applies a bounded, sample-dependent reweighting to the ground-truth rollout gradients rather than introducing an independent distillation direction. We further derive a covariance decomposition of the aggregated update direction and establish conditional finite-step reachability of a prescribed long-horizon risk region. (i) We conduct experiments on nine PDE benchmarks spanning one-, two-, and three-dimensional dynamics with spectral and attention-based backbones. The results show that HERO consistently improves long-horizon rollout accuracy, stability, and out-of-distribution robustness over strong autoregressive-training baselines without increasing inference-time cost. 2. Related Work Rollout-based training reduces the train–test mismatch of autoregressive neural operators (Kovachki et al., 2023; Li et al., 2021) by exposing the model to self-generated states, through push-forward and temporal bundling (Brandstetter et al., 2022), differentiable solver-in-the-loop correction (Um et al., 2020), recursive graph simulators (Sanchez-Gonzalez et al., 2020; Pfaff et al., 2021), unrolling taxonomies (Koehler et al., 2024; List et al., 2025), and recurrent operators (Ye et al., 2025); yet their supervision measures only the absolute discrepancy from the ground-truth trajectory. A complementary direction stabilizes long horizons through spectral and structural design (McCabe et al., 2023; Lippe et al., 2023; Li et al., 2026; Jiang et al., 2025; Huang and Greenberg, 2025) or long-term statistical constraints (Li et al., 2022; Jiang et al., 2023; Schiff et al., 2024), but modifies the operator or objective rather than supervising against the model’s own failures. HERO addresses this complementary gap, leaving the inference-time operator unchanged and using dynamically selected historical failure trajectories to supervise whether the current rollout improves over previously observed failure behavior. The complete introduction of related works is in Appendix A. 3. Motivation This section presents two observations motivating HERO: prediction errors accumulate substantially over long autoregressive rollouts, and historical checkpoints exhibit distinct horizon-dependent behaviors. These findings motivate complementing absolute trajectory supervision with history-enriched relative supervision. 3.1. Long-Horizon Error Accumulation Figure 1 compares FNO, FNO-PF, and HERO over 200-step free autoregressive rollouts from the same test initial conditions. The three methods remain relatively close during the initial interval, but their errors separate substantially as the rollout horizon increases. FNO exhibits rapid error accumulation. FNO-PF delays part of the degradation, yet retains a similar long-range growth trend. In contrast, HERO maintains lower error and slower amplification throughout the medium- and long-horizon rollout. Figure 1 shows the ground-truth and predicted fields at step 100 for a representative test trajectory. The FNO and FNO-PF predictions exhibit clear phase, amplitude, and local oscillatory errors, whereas HERO more closely preserves the spatial structure of the ground-truth field. These results show that short-horizon accuracy and exposure to self-generated states do not by themselves eliminate long-horizon degradation. Absolute supervision measures the discrepancy between the current rollout and the ground truth, but does not explicitly indicate whether the operator has overcome failure behaviors produced earlier during optimization. Figure 1. Short-horizon accuracy versus long-horizon rollout stability of FNO, FNO-PF, and HERO. (a) Mean normalized root-mean-square error (nRMSE) over 200-step free autoregressive rollouts from the same test initial conditions. (b) Ground-truth and predicted fields at step 100 for a representative test trajectory selected using a predefined criterion Figure 2. Relative rollout-error landscape across training checkpoints over 200-step free autoregressive rollouts. Each row shows the base-10 logarithmic ratio between the rollout error of a historical checkpoint and that of Final HERO. The bottom row corresponds to Final HERO and is identically zero. Dashed lines mark rollout steps 40, 100, and 160. 3.2. Distinct Rollout Behaviors across Training Checkpoints Each intermediate checkpoint defines a complete autoregressive operator and therefore induces its own long-horizon rollout behavior. Figure 2 compares checkpoints saved at training steps 4000, 6000, and 8000 with Final HERO. All checkpoints are evaluated from the same test initial conditions using 200-step free autoregressive rollouts. For checkpoint s, the mean relative error at rollout step t is (1) Es(t)=1Ntest∑i=1Ntest‖u^s,t(i)−ut(i)‖2‖ut(i)‖2+ε,E_s(t)= 1N_test _i=1^N_test \| u_s,t^(i)-u_t^(i) \|_2 \|u_t^(i) \|_2+ , where NtestN_test is the number of test trajectories, u^s,t(i) u_s,t^(i) and ut(i)u_t^(i) are the predicted and ground-truth states for trajectory i at rollout step t, respectively, and ε>0 >0 is a numerical-stability constant. We compare each checkpoint with Final HERO through the base-10 logarithmic ratio (2) Δs(t)=log10(Es(t)+εEFinal(t)+ε), _s(t)= _10 ( E_s(t)+ E_Final(t)+ ), where EFinal(t)E_Final(t) is the corresponding error of Final HERO. Positive values of Δs(t) _s(t) indicate a larger rollout error than Final HERO, negative values indicate a smaller error, and values near zero indicate comparable errors. Stacking Δs(t) _s(t) across checkpoints yields the relative error landscape in Figure 2. The horizontal axis represents the rollout step, while the vertical axis represents the training checkpoint. The color scale is centered at zero, and the Final-HERO row is identically zero. The landscape shows that the relative quality of a checkpoint depends on the rollout horizon and does not improve monotonically throughout training. Some historical checkpoints remain competitive during the initial rollout steps but become substantially less accurate as the rollout continues. Different checkpoints also exhibit distinct transition times and error-growth patterns. Historical checkpoints are therefore not merely uniformly weaker versions of the final model; they retain different long-horizon error-propagation behaviors. These observations motivate dynamic history-enriched supervision. Historical rollouts can provide informative failure references, but a fixed checkpoint can become stale, whereas an excessively frequent refresh provides little contrast. HERO therefore maintains a controlled training lag, evaluates candidate rollouts using long-horizon diagnostics, and uses the selected trajectory to complement conventional absolute supervision. Figure 3. Overview of the HERO training framework. A lagged operator, refreshed every R steps, supplies one of three detached candidate rollouts, ranked by rollout, spectral, energy, and error-growth diagnostics to select the failure reference τ−τ^-. The absolute rollout loss is combined with a margin objective driving the current trajectory error below the reference error. 4. Method Figure 3 gives an overview of HERO. From a shared input context, the current operator produces the rollout to be optimized while the lagged and perturbed branches produce detached candidates; the diagnostics then rank these candidates into a failure reference, which meets the current rollout in a margin comparison that is added to the absolute rollout loss. 4.1. Problem Formulation Consider a class of time-dependent PDEs defined on a spatial domain Ω⊆ℝd ^d. The state variable is denoted by u(,t)∈ℝcuu(x,t) ^c_u for ∈Ωx∈ , where cuc_u is the number of physical channels. At the discrete time point tn=nΔt_n=n t, we denote the spatial state by un=u(⋅,tn)u_n=u(·,t_n). Let H denote the number of most recent states provided to the operator at each prediction step. The input context at time step n is defined as zn=(un−H+1,…,un;μ)z_n=(u_n-H+1,…,u_n;μ), where μ denotes conditioning information that may influence the underlying dynamics, such as equation parameters, external forcing, boundary conditions, or geometric information. We consider a differentiable time-stepping operator GθG_θ, parameterized by trainable parameters θ, that maps the current context to a one-step prediction, namely Gθ:zn↦u^n+1G_θ:z_n u_n+1. To generate a length-K autoregressive rollout, the operator is recursively applied to its own predictions. For notational consistency, we initialize u^j=uj u_j=u_j for j=n−H+1,…,nj=n-H+1,…,n. For each rollout step k=1,…,Kk=1,…,K, the subsequent state is generated as (3) u^n+k=Gθ(u^n+k−H,…,u^n+k−1;μ). u_n+k=G_θ ( u_n+k-H,…, u_n+k-1;μ ). The predicted trajectory is denoted by τθ=(u^n+1,…,u^n+K) _θ=( u_n+1,…, u_n+K), and the corresponding ground-truth trajectory is denoted by τ⋆=(un+1,…,un+K) _ =(u_n+1,…,u_n+K). HERO is agnostic to the internal architecture of GθG_θ. It only requires the base operator to support differentiable one-step prediction and recursive rollout. HERO therefore modifies the training procedure without changing the inference-time architecture. 4.2. Absolute Rollout Supervision Given the predicted trajectory τθ _θ and the corresponding ground-truth trajectory τ⋆ _ , we first define a trajectory-level discrepancy. This discrepancy serves as both the conventional absolute supervision and a common measure for evaluating the current and candidate trajectories against the same ground-truth trajectory. Let X denote the discrete state space induced by the chosen spatial discretization. For two states v,w∈v,w , we define the normalized state error as (4) d(v,w)=‖v−w‖w‖+ε,d_X(v,w)= \|v-w\|_X\|w\|_X+ , where ε>0 >0 is a constant for numerical stability. The spatial norm is defined as (5) ‖v‖2=∑i=1Nxωi‖v(i)‖22.\|v\|_X^2= _i=1^N_x _i\|v(x_i)\|_2^2. Here, NxN_x is the number of spatial sampling points ix_i, ∥⋅∥2\|·\|_2 denotes the Euclidean norm over the cuc_u physical channels, and ωi≥0 _i≥ 0 is the discrete integration weight associated with the i-th spatial point. On a regular grid, the weights are uniform. On irregular grids, finite-element meshes, or point clouds, they may represent cell volumes, nodal areas, or quadrature weights. For any length-K trajectory τ=(v1,…,vK)τ=(v_1,…,v_K), we define its distance from the ground-truth trajectory as (6) D(τ,τ⋆)=1K∑k=1Kd(vk,un+k).D(τ, _ )= 1K _k=1^Kd_X (v_k,u_n+k ). Applying this distance to the current rollout gives the absolute rollout objective (7) ℒroll=(zn,τ⋆)[D(τθ,τ⋆)].L_roll=E_(z_n, _ ) [D( _θ, _ ) ]. The expectation is taken over sampled input contexts znz_n and their corresponding ground-truth trajectories τ⋆ _ . The objective ℒrollL_roll anchors the current rollout to the reference PDE evolution. However, it evaluates the current trajectory only through its absolute discrepancy from τ⋆ _ and does not explicitly measure improvement over structured long-horizon behaviors produced during optimization. HERO therefore supplements the absolute objective with history-enriched relative supervision. 4.3. History-Enriched Relative Supervision To provide the missing comparative signal, HERO compares the ground-truth error of the current rollout with that of a dynamically selected, model-generated failure trajectory. Let s denote the optimization iteration and let θs _s denote the current parameters before the parameter update at iteration s. Relative reference construction starts at iteration s0s_0. For a sampled context znz_n, the iteration-indexed current rollout is defined as (8) τθs=Rollout(Gθs,zn,K). _ _s=Rollout (G_ _s,z_n,K ). Here, Rollout(G,z,K)Rollout(G,z,K) denotes the length-K autoregressive rollout of operator G initialized from context z. Thus, τθs _ _s is the optimization-iteration-specific form of the trajectory τθ _θ defined above. The dependence of subsequent trajectory notation on s and n is omitted when no ambiguity arises. HERO constructs the relative supervision in five stages: 1) updating a lagged operator when the refresh condition is satisfied, 2) generating candidate rollouts, 3) computing trajectory diagnostics, 4) scoring and selecting a failure trajectory, and 5) constructing the relative objective. The lagged operator is updated periodically, whereas the remaining stages are performed at each relative-training iteration. These stages correspond to the main operations in Algorithm 1. 1) Lagged Operator: HERO maintains the current operator GθsG_ _s and a lagged operator Gθ¯sG_ θ_s. The current operator is updated by gradient descent, whereas the lagged operator is periodically copied from the current operator and remains frozen between refreshes. Let R be a positive integer denoting the refresh interval. For s≥s0s≥ s_0, the lagged parameters are updated as (9) θ¯s=sg(θs),s=s0,sg(θs),s>s0and(s−s0)modR=0,θ¯s−1,otherwise, θ_s= casessg( _s),&s=s_0,\\[2.0pt] sg( _s),&s>s_0\ and\ (s-s_0) R=0,\\[2.0pt] θ_s-1,&otherwise, cases where sg(⋅)sg(·) denotes the stop-gradient operator. After each refresh, Gθ¯sG_ θ_s remains frozen until the next refresh and receives no backpropagated gradient. Given the same input context znz_n, the lagged operator generates (10) τlag=Rollout(Gθ¯s,zn,K). _lag=Rollout (G_ θ_s,z_n,K ). The lagged rollout is used only as a comparison reference. HERO does not train the current operator to reproduce or imitate τlag _lag. Updating the lagged operator at every iteration would make the current and lagged rollouts nearly identical, whereas keeping it fixed throughout training would produce an increasingly stale reference. Periodic refresh maintains a controlled optimization delay between the two operators. 2) Candidate Rollouts: HERO constructs up to three candidate rollouts: a detached version of the current rollout, a rollout from the lagged operator, and a rollout generated by the current operator from a locally perturbed input context. (i) The self-reference is obtained by detaching the current rollout: (11) τself=sg(τθs). _self=sg ( _ _s ). For a deterministic base operator, τself _self and τθs _ _s are numerically identical, but gradients propagate only through τθs _ _s. The self-reference therefore provides a fixed representation of the current rollout for comparison. (i) The lagged candidate is the rollout τlag _lag generated by the frozen lagged operator. Since the lagged operator receives no backpropagated gradient, τlag _lag provides a detached reference from an earlier optimization state. (i) The perturbed candidate evaluates the current operator near the original input context. Let ζ∼pζ p_ζ denote a random perturbation variable. The operator ζP_ζ perturbs only the state-history component of znz_n and leaves the conditioning information μ unchanged. The perturbed rollout is (12) τpert=sg[Rollout(Gθs,ζ(zn),K)]. _pert=sg [Rollout (G_ _s,P_ζ(z_n),K ) ]. The perturbation preserves the basic scale and boundary structure of the input so that ζ(zn)P_ζ(z_n) remains in a local neighborhood of the original context. On regular grids, ζP_ζ may be implemented as a weak spectral perturbation. On irregular grids or point clouds, it may use local smoothing, neighborhood jitter, or a graph-spectral perturbation. When all three branches are enabled, the active candidate set is (13) =τlag,τpert,τself.C= \ _lag, _pert, _self \. If no suitable perturbation is available, τpert _pert is removed from C. The number of active candidates is M=||M=|C|, with M≥2M≥ 2 because the lagged and self-reference branches remain active. For subsequent evaluation, we fix a deterministic enumeration of the active candidates and denote them by τ(1),…,τ(M)τ^(1),…,τ^(M). The enumeration follows the priority τlag≻τpert≻τself _lag _pert _self, with unavailable branches omitted. All active candidates are detached from the computation graph and are used only for reference selection. 3) Trajectory Diagnostics: This stage evaluates the long-horizon behavior of each candidate and determines which trajectory provides the strongest measured failure reference. For each candidate j=1,…,Mj=1,…,M, its length-K trajectory is τ(j)=(v1(j),…,vK(j))τ^(j)=(v_1^(j),…,v_K^(j)). HERO computes scalar diagnostics that measure the overall rollout error, spectral-amplitude mismatch, quadratic-energy drift, and temporal error growth. These quantities are used only to rank the candidates. The same numerical-stability constant ε>0 >0 is used throughout. The overall rollout error is (14) Eroll(j)=D(τ(j),τ⋆).E_roll^(j)=D (τ^(j), _ ). This diagnostic measures the average state-space discrepancy between the candidate trajectory and the ground-truth trajectory over the rollout horizon. When a spatial spectral transform is available, the spectral-amplitude discrepancy is (15) Espec(j)=1K∑k=1K‖|ℱ[vk(j)]|−|ℱ[un+k]|‖1‖|ℱ[un+k]|‖1+ε.E_spec^(j)= 1K _k=1^K \| |F[v_k^(j)] |- |F[u_n+k] | \|_1 \| |F[u_n+k] | \|_1+ . Here, ℱF denotes the spatial spectral transform applied independently to each physical channel, |⋅||·| denotes the elementwise magnitude, and ∥⋅∥1\|·\|_1 is taken over all spectral coefficients and physical channels. This diagnostic measures whether the candidate assigns incorrect amplitudes to the spatial modes. If no suitable spectral transform is available, this diagnostic is omitted or replaced by one defined in an appropriate spectral basis. To measure global amplitude drift, HERO uses the quadratic state-energy proxy (16) ℰ(u)=∑i=1Nxωi‖u(i)‖22.E(u)= _i=1^N_x _i \|u(x_i) \|_2^2. This quantity equals the squared spatial norm defined above and is used as an energy proxy rather than a PDE-specific conserved energy. The corresponding energy drift is (17) Eenergy(j)=1K∑k=1K|ℰ(vk(j))−ℰ(un+k)|ℰ(un+k)+ε.E_energy^(j)= 1K _k=1^K |E(v_k^(j))-E(u_n+k) |E(u_n+k)+ . This diagnostic measures global amplitude loss, excessive dissipation, and abnormal energy growth. Let ek(j)=d(vk(j),un+k)e_k^(j)=d_X(v_k^(j),u_n+k) denote the state error of candidate j at rollout step k. For K≥2K≥ 2, the error-growth diagnostic is (18) Egrowth(j)=1K−1∑k=1K−1max(0,ek+1(j)−ek(j)).E_growth^(j)= 1K-1 _k=1^K-1 (0,e_k+1^(j)-e_k^(j) ). This diagnostic accumulates only positive error increments and therefore measures how strongly the prediction error grows along the rollout. All candidate trajectories are detached before these diagnostics are computed. Consequently, the diagnostic values are used only for candidate ranking and reference selection, and no gradient is propagated through them. 4) Candidate Scoring and Selection: HERO uses the trajectory diagnostics to score the active candidates and select the candidate that exhibits the strongest measured long-horizon failure. Let D denote the set of active diagnostics. When all four diagnostics are available, it is defined as (19) =roll,spec,energy,growth.D= \roll,spec,energy,growth \. Diagnostics that are not applicable to the current dataset are removed from D. Because the diagnostics have different numerical scales, HERO normalizes each diagnostic across the M active candidates within the same training sample. For each q∈q and j=1,…,Mj=1,…,M, the normalized diagnostic is (20) E~q(j)=Eq(j)−min1≤ℓ≤MEq(ℓ)max1≤ℓ≤MEq(ℓ)−min1≤ℓ≤MEq(ℓ)+ε. E_q^(j)= E_q^(j)- _1≤ ≤ ME_q^( ) _1≤ ≤ ME_q^( )- _1≤ ≤ ME_q^( )+ . This normalization maps each diagnostic to [0,1][0,1] across the active candidates and avoids introducing separate diagnostic-weighting hyperparameters. HERO computes the score of candidate j as the average normalized diagnostic, S(j)=1||∑q∈E~q(j)S^(j)= 1|D| _q E_q^(j) (21). A larger score indicates stronger failure behavior across the active diagnostics. HERO selects the highest-scoring candidate, indexed by j⋆=argmax1≤j≤MS(j)j = _1≤ j≤ MS^(j). If multiple candidates attain the same maximum, HERO selects the first one according to the fixed candidate order defined above. The selected candidate is used as the failure reference: (22) τ−=τ(j⋆).τ^-=τ^(j ). The selected τ−τ^- is therefore a model-generated trajectory that exhibits the strongest measured failure behavior among the active candidates. The diagnostics determine which candidate is selected, whereas the relative objective compares the selected trajectory with the current rollout using the common trajectory distance D. The diagnostic values and candidate scores are not backpropagated. 5) Relative Objective: HERO converts the selected failure trajectory τ−τ^- into a relative improvement signal. The current trajectory error is dcur=D(τθs,τ⋆)d_cur=D( _ _s, _ ), while the selected reference error is dfail=sg[D(τ−,τ⋆)]d_fail=sg[D(τ^-, _ )]. The stop-gradient operation treats the selected trajectory as a fixed comparison reference, so gradients propagate only through the current rollout. HERO encourages the current rollout to outperform the selected reference by a margin m>0m>0, corresponding to dcur+m≤dfaild_cur+m≤ d_fail. We implement this requirement using the smooth margin loss (23) ℓrel=1βlog(1+exp(β[dcur−dfail+m])). _rel= 1β (1+ (β [d_cur-d_fail+m ] ) ). where β>0β>0 controls the sharpness of the transition around the margin. The loss assigns a larger penalty when the current rollout has not outperformed the reference by the required margin and decreases as the margin is satisfied. Averaging the sample-level loss over the training distribution gives the relative supervision objective (24) ℒrel=(zn,τ⋆)[ℓrel].L_rel=E_(z_n, _ ) [ _rel ]. When the perturbed branch is enabled, the expectation also includes the randomness of ζ. 4.4. Training Objective HERO introduces relative supervision only after the operator has learned a basic time-stepping model. Reference construction is enabled at iteration s0s_0, while the relative-supervision weight increases from zero over a warmup period. Let W∈ℕ+W _+ denote the warmup length and let λmax≥0 _ ≥ 0 denote the maximum relative-supervision weight. At optimization iteration s, the weight is (25) λs=λmaxclip(s−s0W,0,1), _s= _ clip ( s-s_0W,0,1 ), where clip(x,a,b)=min(max(x,a),b)clip(x,a,b)= ( (x,a),b). Thus, λs=0 _s=0 for s≤s0s≤ s_0, increases from zero to λmax _ over the following W iterations, and remains equal to λmax _ afterward. The complete training objective is (26) ℒtotal=ℒroll+λsℒrel.L_total=L_roll+ _sL_rel. The absolute term learns the reference PDE evolution, while the relative term encourages the current rollout to improve upon the selected failure trajectory. 5. Experiments 5.1. Experimental Setup Backbones and baselines. We build every method on two backbones: FNO, a spectral neural operator, and Transolver, an attention-based neural operator with Physics-Attention (Wu et al., 2024). On each backbone we compare one-step supervised training, push-forward training (PF), PDE-Refiner (Refine), recurrent operator training (RNO), and HERO, denoting variants by Backbone-Strategy (e.g., FNO-PF) (Brandstetter et al., 2022; Lippe et al., 2023; Ye et al., 2025). PF and RNO are closest to HERO, since all three train on self-generated states yet supervise only the absolute discrepancy from the ground truth, whereas Refine corrects the output itself; HERO keeps ground-truth rollout regression and adds relative supervision against a dynamically history-enriched failure reference. Full baseline descriptions are given in Appendix C.1.1. Benchmarks. We evaluate on nine PDE benchmarks spanning 1D, 2D, and 3D dynamics (Table 4); APEBench supplies the task definitions, while PDEBench provides complementary scientific-machine-learning benchmark coverage (Koehler et al., 2024; Takamoto et al., 2022): Dispersion (Korteweg and de Vries, 1895), Burgers (Burgers, 1948), KdV (Korteweg and de Vries, 1895), and Kuramoto–Sivashinsky (KS) (Sivashinsky, 1977; Kuramoto, 1978) in 1D; Anisotropic Diffusion (Fick, 1855), Kolmogorov Flow (Meshalkin and Sinai, 1961), and Decaying Turbulence (Navier–Stokes) (Navier and others, 1823) in 2D; and Swift Hohenberg (Swift and Hohenberg, 1977) and Unbalanced Advection in 3D (Courant et al., 1928). Together they cover dispersion, nonlinear transport, dissipative chaos, forced and decaying turbulence, pattern formation, and multi-directional advection, stressing phase, spectral, energy, and stability behavior over long rollouts. All benchmarks use periodic boundaries and a single scalar channel, with 5050 training and 3030 test trajectories and 200200-step free autoregressive rollouts. Governing equations and coefficients are given in Appendix C.1.2. Metrics. For evaluation trajectory j and rollout step t, the per-step normalized RMSE is nRMSEt(j)=‖u^t(j)−ut(j)‖2/(‖ut(j)‖2+ε)nRMSE_t^(j)=\| u_t^(j)-u_t^(j)\|_2/(\|u_t^(j)\|_2+ ), with ε>0 >0 a numerical-stability constant. We report the set-averaged nRMSE@1nRMSE@1 and nRMSE@100nRMSE@100 (one-step and 100100-step accuracy; lower is better), the geometric-mean rollout error (GM100) (27) GM100(j)=exp(1100∑t=1100log(nRMSEt(j)+ε)),GM100^(j)= \! ( 1100 _t=1^100 (nRMSE_t^(j)+ ) ), reported as its validation- and test-set means (Val/Test mean GM100; lower is better), and the mean stable step 1M∑jS(j) 1M _jS^(j) with S(j)=maxT∈[1,100]∣nRMSEt(j)≤τ∀t≤TS^(j)= \T∈[1,100] _t^(j)≤τ\ \ ∀\,t≤ T\ and τ=0.1τ=0.1 (higher is better). Val mean GM100 is used for checkpoint selection. The detailed metrics descriptions are given in Appendix C.1.3 Implementation. Each task provides 50 training trajectories of length 51 and 30 test trajectories of length 201; at test time every model rolls out 200200 steps from the first frame. All methods share the backbone, data split, optimizer, and evaluation protocol, and use a Kpf=5K_pf=5 push-forward base; HERO adds the relative objective ℒrelL_rel, giving ℒtotal=ℒroll+λsℒrelL_total=L_roll+ _sL_rel. Across all tasks HERO uses K=5K=5, s0=4000s_0=4000, R=2000R=2000, W=2000W=2000, λmax=0.02 _ =0.02, m=0.02m=0.02, β=5.0β=5.0, and candidate set =τlag,τpert,τselfC=\ _lag, _pert, _self\; reference construction is confined to training and adds no inference-time cost. We optimize with AdamW (learning rate 10−310^-3, warm-up–cosine schedule, 20002000 warm-up steps, 10 00010\,000 total steps, batch size 2020), and report test results at the checkpoint with the best Val mean GM100. Backbone architectures per dimension and further details are given in Appendix C.1.4. We fix the start step s0=4000s_0=4000 and the refresh interval R=2000R=2000 of (9) for all tasks: this pair introduces relative supervision once the backbone has acquired a stable one-step evolution while leaving 60006000 steps of history-enriched reweighting, and it refreshes the lagged operator θ¯s θ_s often enough to avoid a stale reference without collapsing it onto the current operator θs _s. A sensitivity analysis of both hyperparameters is given in Appendix C.3. 5.2. Main Results Table 1 reports the mean performance on six PDE benchmarks spanning one-, two-, and three-dimensional systems. The complete results on all nine benchmarks, including the remaining three benchmarks and the corresponding standard deviations over five seeds, are provided in Tables 9 and 10 in Appendix C.2. Across both the main and additional benchmarks, HERO consistently improves FNO and Transolver. On the six main benchmarks, HERO achieves the best result within each backbone family for all reported metrics. On FNO, HERO reduces nRMSE@100 from 0.5410.541 to 0.1530.153 on Burgers, from 1.1311.131 to 0.6400.640 on KS, from 1.8681.868 to 0.5210.521 on anisotropic diffusion, and from 13.62913.629 to 4.1014.101 on NS. Compared with the strongest non-HERO strategy, it reduces nRMSE@100 by 33.7%33.7\%–35.8%35.8\% across the six benchmarks. Similarly, on Transolver, HERO reduces nRMSE@100 by 33.6%33.6\%–36.7%36.7\% relative to the strongest competing strategy. The improvements in nRMSE@100 and Stable step are substantially larger than those in nRMSE@1. For example, HERO increases the stable rollout length from 35.635.6 to 81.781.7 steps on Burgers and from 8.28.2 to 26.626.6 steps on anisotropic diffusion. These results indicate that HERO primarily suppresses long-horizon error accumulation rather than merely improving one-step fitting, while its consistent gains on both backbones demonstrate architectural generality. Because reference construction is confined to training, HERO leaves the deployed operator identical to its backbone; a measured comparison of parameter count, inference latency, memory, and training throughput is reported in Appendix C.4. Table 1. Main results on six PDE benchmarks. Method Val GM100 ↓ Test GM100 ↓ nRMSE @1 ↓ nRMSE @100 ↓ Stable step ↑ Burgers (1D) FNO 0.473 0.526 0.0429 0.541 35.6 FNO-PF 0.218 0.208 0.0180 0.236 68.0 FNO-Ref. 0.407 0.457 0.0368 0.495 40.4 FNO-RNO 0.293 0.311 0.0249 0.335 55.0 FNO-HERO 0.147 0.138 0.0168 0.153 81.7 Trans. 0.487 0.537 0.0478 0.618 31.2 Trans.-PF 0.449 0.505 0.0422 0.562 35.5 Trans.-Ref. 0.428 0.483 0.0401 0.535 38.7 Trans.-RNO 0.384 0.434 0.0362 0.494 41.4 Trans.-HERO 0.250 0.281 0.0341 0.327 49.3 Kuramoto–Sivashinsky (1D) FNO 0.511 0.578 0.0190 1.131 43.8 FNO-PF 0.401 0.451 0.0179 0.972 44.7 FNO-Ref. 0.469 0.528 0.0159 1.064 46.4 FNO-RNO 0.458 0.515 0.0181 1.009 43.6 FNO-HERO 0.257 0.283 0.0150 0.640 53.3 Trans. 0.611 0.695 0.0250 1.284 37.6 Trans.-PF 0.501 0.571 0.0211 1.163 41.8 Trans.-Ref. 0.498 0.572 0.0190 1.238 42.4 Trans.-RNO 0.477 0.538 0.0201 1.099 42.6 Trans.-HERO 0.321 0.361 0.0181 0.715 49.2 Anisotropic Diffusion (2D) FNO 0.713 0.825 0.0875 1.868 8.2 FNO-PF 0.655 0.754 0.0716 1.355 12.9 FNO-Ref. 0.567 0.652 0.0611 0.786 20.7 FNO-RNO 0.692 0.778 0.0846 0.925 9.7 FNO-HERO 0.363 0.415 0.0571 0.521 26.6 Trans. 1.168 1.358 0.0799 2.968 8.1 Trans.-PF 0.965 1.111 0.0696 2.172 13.0 Trans.-Ref. 0.803 0.943 0.0692 1.835 14.7 Trans.-RNO 0.892 1.048 0.0751 2.099 10.4 Trans.-HERO 0.537 0.626 0.0652 1.185 18.2 Navier–Stokes (2D) FNO 2.322 2.915 0.0449 13.629 17.6 FNO-PF 1.441 1.682 0.0349 6.271 22.3 FNO-Ref. 1.959 2.330 0.0402 9.888 18.1 FNO-RNO 1.788 2.141 0.0388 7.829 20.7 FNO-HERO 0.931 1.058 0.0325 4.101 26.5 Trans. 3.075 3.825 0.0301 18.084 20.3 Trans.-PF 2.400 2.974 0.0278 12.438 24.3 Trans.-Ref. 2.169 2.591 0.0280 10.974 23.1 Trans.-RNO 2.224 2.704 0.0291 11.438 22.5 Trans.-HERO 1.423 1.687 0.0262 7.012 25.8 Swift–Hohenberg (3D) FNO 0.603 0.678 0.0788 0.941 13.4 FNO-PF 0.514 0.581 0.0607 0.843 20.4 FNO-Ref. 0.552 0.619 0.0642 0.871 19.6 FNO-RNO 0.520 0.596 0.0662 0.850 19.5 FNO-HERO 0.341 0.381 0.0570 0.541 27.7 Trans. 0.821 0.912 0.1079 1.286 5.3 Trans.-PF 0.686 0.782 0.0914 1.126 5.9 Trans.-Ref. 0.633 0.713 0.0823 1.004 11.3 Trans.-RNO 0.652 0.739 0.0860 1.064 7.6 Trans.-HERO 0.404 0.450 0.0780 0.667 15.0 Unbalanced Advection (3D) FNO 2.311 2.647 0.0869 7.204 6.4 FNO-PF 1.833 2.121 0.0748 5.546 10.0 FNO-Ref. 2.017 2.317 0.0776 6.086 7.8 FNO-RNO 1.697 1.931 0.0791 4.833 9.0 FNO-HERO 1.101 1.244 0.0702 3.185 12.7 Trans. 2.723 3.106 0.0731 7.777 9.8 Trans.-PF 2.096 2.426 0.0736 6.122 10.3 Trans.-Ref. 1.901 2.201 0.0713 5.319 9.6 Trans.-RNO 1.784 2.079 0.0742 5.015 10.6 Trans.-HERO 1.197 1.389 0.0673 3.175 12.1 Stable step denotes the mean number of stable rollout steps. 5.3. Ablation Study To examine which part of the history-enriched design produces the long-horizon gain, we conduct component-level ablations on the 2D Navier–Stokes benchmark. All variants share the FNO backbone, data split, optimizer, evaluation protocol, and the Kpf=5K_pf=5 push-forward absolute objective; only two factors are modified: whether the relative objective ℒrelL_rel is active, and which branches remain in the candidate set C in (13). FNO and FNO-PF are the one-step and push-forward baselines (λmax=0 _ =0). Self-relative and Lag-only HERO restrict C to τself\ _self\ and τlag\ _lag\, so the reference is fixed rather than selected. HERO w/o Lag, w/o Pert, and w/o Self each drop one branch from the default set, whereas HERO w/o Relative keeps candidate construction and diagnostics but sets λs=0 _s=0 in (26). Table 2. Component-level ablation of HERO on the 2D Navier–Stokes benchmark Variant Val GM100 ↓ Test GM100 ↓ nRMSE @1 ↓ nRMSE @100 ↓ Stable step ↑ FNO 2.3222.322 2.9152.915 0.04490.0449 13.62913.629 17.617.6 FNO-PF 1.4411.441 1.6821.682 0.03490.0349 6.2716.271 22.322.3 Self-relative 1.3271.327 1.5741.574 0.03430.0343 5.7145.714 23.623.6 Lag-only HERO 1.1681.168 1.3461.346 0.03360.0336 4.8734.873 24.824.8 HERO w/o Lag 1.2051.205 1.3971.397 0.03350.0335 5.0265.026 24.524.5 HERO w/o Pert 1.0611.061 1.2211.221 0.03290.0329 4.5274.527 25.425.4 HERO w/o Self 1.0241.024 1.1781.178 0.03310.0331 4.3814.381 25.725.7 HERO w/o Relative 1.4251.425 1.6641.664 0.03480.0348 6.1986.198 22.522.5 Full HERO 0.9310.931 1.0581.058 0.03250.0325 4.1014.101 26.526.5 Table 2 reports the ablation results. HERO w/o Relative performs at the level of FNO-PF (Test GM100: 1.6641.664 vs. 1.6821.682), confirming that candidate construction and diagnostics contribute nothing by themselves and that the gain comes from the relative gradient. Activating it reduces Test GM100 from 1.6641.664 to 1.0581.058 and nRMSE@100 from 6.1986.198 to 4.1014.101, approximately 36.4%36.4\% and 33.8%33.8\%, while nRMSE@1 changes only from 0.03480.0348 to 0.03250.0325, showing that relative supervision acts on error accumulation rather than on one-step fitting. Self-relative improves over FNO-PF but remains far from Full HERO (1.5741.574 vs. 1.0581.058), since a reference drawn from the current parameters provides little contrast, and Lag-only HERO recovers only part of the gap (1.3461.346) with a fixed-lag reference. Among the reduced candidate sets, HERO w/o Lag degrades most (1.3971.397), whereas w/o Pert (1.2211.221) and w/o Self (1.1781.178) lose less; although w/o Self comes closest overall, its shorter stable horizon (25.725.7 vs. 26.526.5) localizes the residual difference in late-rollout stability. Overall, these results show that the long-horizon gain of HERO arises from relative supervision against a dynamically history-enriched failure reference, with the lagged branch and diagnostic-based selection both necessary for the full effect. 5.4. Out-of-Distribution Generalization To test whether HERO remains stable outside the training regime, we conduct zero-shot evaluation on the 2D Kolmogorov flow. All models are trained at ν=1.00×10−2ν=1.00× 10^-2, Af=1.0A_f=1.0, and N=642N=64^2, then evaluated without fine-tuning under i) parameter shifts to ν=6.67×10−3ν=6.67× 10^-3 or Af=1.50A_f=1.50 and i) resolution shifts to 96296^2 or 1282128^2. All methods use the same initial conditions and rollout protocol. We report GM100, nRMSE@200, and the mean stability horizon (MSH). Table 3 reports that HERO achieves the lowest error and longest stability horizon under all four shifts. Under parameter extrapolation, HERO reduces GM100 from 1.4871.487 to 0.7420.742 for lower viscosity and from 1.6031.603 to 0.8090.809 for stronger forcing relative to FNO. Compared with FNO-PF, the strongest baseline, HERO reduces GM100 by 37.9%37.9\% and 38.6%38.6\%, respectively, while approximately doubling MSH. HERO also maintains its advantage under resolution extrapolation. At 96296^2 and 1282128^2, it reduces GM100 from 1.2741.274 to 0.6540.654 and from 1.3611.361 to 0.7040.704 relative to FNO, respectively. Compared with FNO-PF, these values correspond to reductions of 37.2%37.2\% and 37.5%37.5\%, while MSH increases from 1.91.9 to 3.53.5 and from 1.61.6 to 3.03.0. These results indicate that HERO suppresses long-horizon error accumulation under unseen physical parameters and spatial resolutions. Table 3. Zero-shot out-of-distribution generalization on the 2D Kolmogorov flow under parameter and resolution shift Method GM100 ↓ nRMSE@200 ↓ MSH ↑ Parameter: ν=6.67×10−3ν=6.67× 10^-3, Af=1.0A_f=1.0, N=642N=64^2 FNO 1.4871.487 2.2862.286 1.21.2 FNO-PF 1.1941.194 1.9361.936 1.71.7 FNO-RNO 1.2661.266 2.0442.044 1.51.5 FNO-HERO 0.7420.742 1.2111.211 3.13.1 Parameter: ν=1.00×10−2ν=1.00× 10^-2, Af=1.50A_f=1.50, N=642N=64^2 FNO 1.6031.603 2.4722.472 1.01.0 FNO-PF 1.3181.318 2.1272.127 1.41.4 FNO-RNO 1.3741.374 2.2162.216 1.31.3 FNO-HERO 0.8090.809 1.3361.336 2.72.7 Resolution: ν=1.00×10−2ν=1.00× 10^-2, Af=1.0A_f=1.0, N=962N=96^2 FNO 1.2741.274 1.7541.754 1.31.3 FNO-PF 1.0411.041 1.4931.493 1.91.9 FNO-RNO 1.0961.096 1.5721.572 1.71.7 FNO-HERO 0.6540.654 0.9460.946 3.53.5 Resolution: ν=1.00×10−2ν=1.00× 10^-2, Af=1.0A_f=1.0, N=1282N=128^2 FNO 1.3611.361 1.9671.967 1.11.1 FNO-PF 1.1261.126 1.6811.681 1.61.6 FNO-RNO 1.1871.187 1.7741.774 1.51.5 FNO-HERO 0.7040.704 1.0581.058 3.03.0 6. Conclusion HERO introduced history-enriched rollout training for long-horizon autoregressive neural operators. Instead of relying only on absolute trajectory error, HERO compares the current rollout with dynamically selected failure trajectories from earlier optimization states and converts this comparison into a margin-based relative objective. Experiments across nine PDE benchmarks with spectral and attention-based backbones showed consistent improvements in long-horizon accuracy, stable rollout length, and out-of-distribution robustness, without changing the inference-time architecture or cost. Ablation studies further confirmed that the gains arise from relative supervision and that the lagged branch and diagnostic-based reference selection are necessary for the full improvement. These results indicate that optimization history provides an effective training signal for suppressing autoregressive error accumulation in neural PDE surrogates. 7. Limitations and Ethical Considerations This work uses only simulated PDE data and does not involve human subjects, personal information, or privacy-sensitive content. The current evaluation focuses on the PDE benchmarks and backbone architectures. HERO also introduces additional training-time computation for candidate rollout construction and reference selection, while leaving the inference-time architecture and cost unchanged. 8. Generative AI Usage Generative AI was used only for language editing. References S. Bengio, O. Vinyals, N. Jaitly, and N. Shazeer (2015) Scheduled sampling for sequence prediction with recurrent neural networks. In Advances in Neural Information Processing Systems, Vol. 28. Cited by: §A.1. Y. Bengio, J. Louradour, R. Collobert, and J. Weston (2009) Curriculum learning. In Proceedings of the 26th Annual International Conference on Machine Learning, ICML ’09, New York, NY, USA, p. 41–48. External Links: ISBN 9781605585161, Link, Document Cited by: §A.1. B. Bonev, T. Kurth, C. Hundt, J. Pathak, M. Baust, K. Kashinath, and A. Anandkumar (2023) Spherical Fourier neural operators: learning stable dynamics on the sphere. In Proceedings of the 40th International Conference on Machine Learning, A. Krause, E. Brunskill, K. Cho, B. Engelhardt, S. Sabato, and J. Scarlett (Eds.), Proceedings of Machine Learning Research, Vol. 202, p. 2806–2823. External Links: Link Cited by: §A.1. J. Brandstetter, D. E. Worrall, and M. Welling (2022) Message passing neural PDE solvers. In The Tenth International Conference on Learning Representations, ICLR 2022, Virtual Event, April 25-29, 2022, Cited by: §A.1, §1, §2, §5.1. S. L. Brunton, B. R. Noack, and P. Koumoutsakos (2020) Machine learning for fluid mechanics. Annual review of fluid mechanics 52 (1), p. 477–508. Cited by: §A.2. J.M. Burgers (1948) A mathematical model illustrating the theory of turbulence. R. Von Mises and T. Von Kármán (Eds.), Advances in Applied Mechanics, Vol. 1, p. 171–199. External Links: ISSN 0065-2156, Document, Link Cited by: §5.1. Q. Cao, S. Goswami, and G. E. Karniadakis (2024) Laplace neural operator for solving differential equations. Nature Machine Intelligence 6 (6), p. 631–640. External Links: Document, Link, ISBN 2522-5839 Cited by: §A.1. S. Cao (2021) Choose a transformer: fourier or galerkin. In Advances in Neural Information Processing Systems, Vol. 34. Cited by: §A.1. W. Chen, J. Song, P. Ren, S. Subramanian, D. Morozov, and M. W. Mahoney (2024) Data-efficient operator learning via unsupervised pretraining and in-context learning. In Advances in Neural Information Processing Systems, Vol. 37. Cited by: §A.1. R. Courant, K. Friedrichs, and H. Lewy (1928) Über die partiellen differenzengleichungen der mathematischen physik. Mathematische Annalen 100 (1), p. 32–74. External Links: Document, Link, ISBN 1432-1807 Cited by: §5.1. A. Fick (1855) Ueber diffusion. Annalen der physik 170 (1), p. 59–86. Cited by: §5.1. G. Gupta, X. Xiao, and P. Bogdan (2021) Multiwavelet-based operator learning for differential equations. In Advances in Neural Information Processing Systems, Vol. 34, p. 24048–24062. External Links: Link Cited by: §A.1. Z. Hao, C. Su, S. Liu, J. Berner, C. Ying, H. Su, A. Anandkumar, J. Song, and J. Zhu (2024) DPOT: auto-regressive denoising operator transformer for large-scale PDE pre-training. In Proceedings of the 41st International Conference on Machine Learning, R. Salakhutdinov, Z. Kolter, K. Heller, A. Weller, N. Oliver, J. Scarlett, and F. Berkenkamp (Eds.), Proceedings of Machine Learning Research, Vol. 235, p. 17616–17635. External Links: Link Cited by: §A.1. Z. Hao, Z. Wang, H. Su, C. Ying, Y. Dong, S. Liu, Z. Cheng, J. Song, and J. Zhu (2023) GNOT: a general neural operator transformer for operator learning. In Proceedings of the 40th International Conference on Machine Learning, A. Krause, E. Brunskill, K. Cho, B. Engelhardt, S. Sabato, and J. Scarlett (Eds.), Proceedings of Machine Learning Research, Vol. 202, p. 12556–12569. External Links: Link Cited by: §A.1. M. Herde, B. Raonić, T. Rohner, R. Käppeli, R. Molinaro, E. de Bézenac, and S. Mishra (2024) Poseidon: efficient foundation models for pdes. In Advances in Neural Information Processing Systems, Vol. 37. Cited by: §A.1. Y. Huang and D. S. Greenberg (2025) Geometric and physical constraints synergistically enhance neural PDE surrogates. In Forty-second International Conference on Machine Learning, ICML 2025, Vancouver, BC, Canada, July 13-19, 2025, Vol. 267. Cited by: §A.2, §2. R. Jiang, P. Y. Lu, E. Orlova, and R. Willett (2023) Training neural operators to preserve invariant measures of chaotic attractors. In Advances in Neural Information Processing Systems, Vol. 36, p. 27645–27669. Cited by: §A.2, §2. R. Jiang, X. Zhang, K. Jakhar, P. Y. Lu, P. Hassanzadeh, M. Maire, and R. Willett (2025) Hierarchical implicit neural emulators. In Advances in Neural Information Processing Systems, Vol. 38, p. 73718–73751. Cited by: §A.2, §1, §2. G. E. Karniadakis, I. G. Kevrekidis, L. Lu, P. Perdikaris, S. Wang, and L. Yang (2021) Physics-informed machine learning. Nature Reviews Physics 3, p. 422–440. External Links: Document, Link Cited by: §A.1, §A.2. D. Kochkov, J. A. Smith, A. Alieva, Q. Wang, M. P. Brenner, and S. Hoyer (2021) Machine learning–accelerated computational fluid dynamics. Proceedings of the National Academy of Sciences 118 (21), p. e2101784118. External Links: Document, Link Cited by: §A.2. F. Koehler, S. Niedermayr, R. Westermann, and N. Thuerey (2024) APEBench: a benchmark for autoregressive neural emulators of pdes. In Advances in Neural Information Processing Systems, Vol. 37, p. 120252–120310. Cited by: §A.1, §2, §5.1. D. J. Korteweg and G. de Vries (1895) XLI. on the change of form of long waves advancing in a rectangular canal, and on a new type of long stationary waves. The London, Edinburgh, and Dublin Philosophical Magazine and Journal of Science 39 (240), p. 422–443. External Links: Document, Link, https://doi.org/10.1080/14786449508620739 Cited by: §5.1. N. Kovachki, Z. Li, B. Liu, K. Azizzadenesheli, K. Bhattacharya, A. Stuart, and A. Anandkumar (2023) Neural operator: learning maps between function spaces with applications to pdes. Journal of Machine Learning Research 24 (89), p. 1–97. Cited by: §A.1, §2. Y. Kuramoto (1978) Diffusion-induced chaos in reaction systems. Progress of Theoretical Physics Supplement 64, p. 346–367. External Links: ISSN 0375-9687, Document, Link, https://academic.oup.com/ptps/article-pdf/doi/10.1143/PTPS.64.346/5293041/64-346.pdf Cited by: §5.1. A. Lamb, A. Goyal, Y. Zhang, S. Zhang, A. Courville, and Y. Bengio (2016) Professor forcing: a new algorithm for training recurrent networks. In Advances in Neural Information Processing Systems, Vol. 29. Cited by: §A.1. J. Li, P. Jiang, H. Saleem, Z. Wang, P. Koniusz, and F. D. Salim (2026) SGNO: spectral generator neural operators for stable long horizon pde rollouts. arXiv preprint arXiv:2602.18801. Cited by: §A.2, §1, §2. Z. Li, D. Shu, and A. Barati Farimani (2023a) Scalable transformer for pde surrogate modeling. In Advances in Neural Information Processing Systems, Vol. 36, p. 28010–28039. Cited by: §1. Z. Li, D. Z. Huang, B. Liu, and A. Anandkumar (2023b) Fourier neural operator with learned deformations for pdes on general geometries. Journal of Machine Learning Research 24 (388), p. 1–26. External Links: Link Cited by: §A.1. Z. Li, N. Kovachki, K. Azizzadenesheli, B. Liu, K. Bhattacharya, A. Stuart, and A. Anandkumar (2021) Fourier neural operator for parametric partial differential equations. In International Conference on Learning Representations, External Links: Link Cited by: §A.1, §2. Z. Li, M. Liu-Schiaffini, N. Kovachki, K. Azizzadenesheli, B. Liu, K. Bhattacharya, A. Stuart, and A. Anandkumar (2022) Learning chaotic dynamics in dissipative systems. In Advances in Neural Information Processing Systems, Vol. 35, p. 16768–16781. Cited by: §A.2, §1, §2. Z. Li, H. Zheng, N. Kovachki, D. Jin, H. Chen, B. Liu, K. Azizzadenesheli, and A. Anandkumar (2024) Physics-informed neural operator for learning partial differential equations. ACM/IMS Journal of Data Science 1 (3). External Links: Link, Document Cited by: §A.1. P. Lippe, B. Veeling, P. Perdikaris, R. Turner, and J. Brandstetter (2023) Pde-refiner: achieving accurate long rollouts with neural pde solvers. In Advances in Neural Information Processing Systems, Vol. 36, p. 67398–67433. Cited by: §A.2, §1, §2, §5.1. B. List, L. Chen, K. Bali, and N. Thuerey (2025) Differentiability in unrolled training of neural physics simulators on transient dynamics. Computer Methods in Applied Mechanics and Engineering 433, p. 117441. Cited by: §A.1, §1, §2. M. McCabe, P. Harrington, S. Subramanian, and J. Brown (2023) Towards stability of autoregressive neural operators. Trans. Mach. Learn. Res. 2023. External Links: Link Cited by: §A.2, §1, §2. L.D. Meshalkin and Ia.G. Sinai (1961) Investigation of the stability of a stationary solution of a system of equations for the plane movement of an incompressible viscous liquid. Journal of Applied Mathematics and Mechanics 25 (6), p. 1700–1705. External Links: ISSN 0021-8928, Document, Link Cited by: §5.1. C. Navier et al. (1823) Mémoire sur les lois du mouvement des fluides. Mémoires de l’Académie Royale des Sciences de l’Institut de France 6 (1823), p. 389–440. Cited by: §5.1. T. Pfaff, M. Fortunato, A. Sanchez-Gonzalez, and P. W. Battaglia (2021) Learning mesh-based simulation with graph networks. In 9th International Conference on Learning Representations, ICLR 2021, Virtual Event, Austria, May 3-7, 2021, Cited by: §A.1, §2. M. Raissi, P. Perdikaris, and G.E. Karniadakis (2019) Physics-informed neural networks: a deep learning framework for solving forward and inverse problems involving nonlinear partial differential equations. Journal of Computational Physics 378, p. 686–707. External Links: ISSN 0021-9991, Document, Link Cited by: §A.1, §A.2. B. Raonić, R. Molinaro, T. D. Ryck, T. Rohner, F. Bartolucci, R. Alaifari, S. Mishra, and E. de Bézenac (2023) Convolutional neural operators for robust and accurate learning of pdes. In Advances in Neural Information Processing Systems, Vol. 36. Cited by: §A.1. S. Rasp, P. D. Dueben, S. Scher, J. A. Weyn, S. Mouatadid, and N. Thuerey (2020) WeatherBench: a benchmark data set for data-driven weather forecasting. Journal of Advances in Modeling Earth Systems 12, p. e2020MS002203. External Links: Document, Link Cited by: §A.2. M. Ren, W. Zeng, B. Yang, and R. Urtasun (2018) Learning to reweight examples for robust deep learning. In Proceedings of the 35th International Conference on Machine Learning, J. Dy and A. Krause (Eds.), Proceedings of Machine Learning Research, Vol. 80, p. 4334–4343. External Links: Link Cited by: §A.1. S. Ross, G. Gordon, and D. Bagnell (2011) A reduction of imitation learning and structured prediction to no-regret online learning. In Proceedings of the Fourteenth International Conference on Artificial Intelligence and Statistics, G. Gordon, D. Dunson, and M. Dudík (Eds.), Proceedings of Machine Learning Research, Vol. 15, Fort Lauderdale, FL, USA, p. 627–635. External Links: Link Cited by: §A.1. A. Sanchez-Gonzalez, J. Godwin, T. Pfaff, R. Ying, J. Leskovec, and P. Battaglia (2020) Learning to simulate complex physics with graph networks. In International conference on machine learning, p. 8459–8468. Cited by: §A.1, §2. Y. Schiff, Z. Y. Wan, J. B. Parker, S. Hoyer, V. Kuleshov, F. Sha, and L. Zepeda-Núñez (2024) DySLIM: dynamics stable learning by invariant measure for chaotic systems. In Forty-first International Conference on Machine Learning, ICML 2024, Vienna, Austria, July 21-27, 2024, Vol. 235, p. 43649–43684. Cited by: §A.2, §1, §2. G.I. Sivashinsky (1977) Nonlinear analysis of hydrodynamic instability in laminar flames—i. derivation of basic equations. Acta Astronautica 4 (11), p. 1177–1206. External Links: ISSN 0094-5765, Document, Link Cited by: §5.1. J. Swift and P. C. Hohenberg (1977) Hydrodynamic fluctuations at the convective instability. Phys. Rev. A 15, p. 319–328. External Links: Document, Link Cited by: §5.1. M. Takamoto, T. Praditia, R. Leiteritz, D. MacKinlay, F. Alesiani, D. Pflüger, and M. Niepert (2022) PDEBENCH: an extensive benchmark for scientific machine learning. In Advances in Neural Information Processing Systems, Vol. 35. Cited by: §5.1. A. Tran, A. Mathews, L. Xie, and C. S. Ong (2023) Factorized fourier neural operators. In The Eleventh International Conference on Learning Representations, External Links: Link Cited by: §A.1. K. Um, R. Brand, Y. R. Fei, P. Holl, and N. Thuerey (2020) Solver-in-the-loop: learning from differentiable physics to interact with iterative pde-solvers. In Advances in Neural Information Processing Systems, Vol. 33, p. 6111–6122. Cited by: §A.1, §1, §2. H. Wu, H. Luo, H. Wang, J. Wang, and M. Long (2024) Transolver: a fast transformer solver for pdes on general geometries. In International Conference on Machine Learning, Cited by: §5.1. Z. Ye, C. Zhang, and W. Wang (2025) Recurrent neural operators: stable long-term pde prediction. arXiv preprint arXiv:2505.20721. Cited by: §A.1, §1, §2, §5.1. Y. D. Zhong, B. Dey, and A. Chakraborty (2020) Symplectic ode-net: learning hamiltonian dynamics with control. In International Conference on Learning Representations, External Links: Link Cited by: §A.2. Appendix A Related Works A.1. Rollout Training with Self-Generated States Neural operators learn mappings between function spaces and provide efficient surrogate models for parameterized PDEs (Kovachki et al., 2023). The FNO parameterizes global integral kernels in the Fourier domain and is widely used as a backbone for time-dependent PDE modeling (Li et al., 2021). Other operator families use multiwavelet bases, Fourier- or Galerkin-style attention, and factorized Fourier layers (Gupta et al., 2021; Cao, 2021; Tran et al., 2023). Learned spatial deformations, spherical representations, and Laplace-domain modeling adapt operators to different geometries or dynamics (Li et al., 2023b; Bonev et al., 2023; Cao et al., 2024). Continuous convolutions and attention-based neural operators provide further architectural alternatives (Raonić et al., 2023; Hao et al., 2023). Physics-informed methods encode governing-equation constraints in neural surrogates (Raissi et al., 2019; Li et al., 2024; Karniadakis et al., 2021). Pretraining and foundation-model approaches instead emphasize transfer across PDEs and data efficiency (Hao et al., 2024; Herde et al., 2024; Chen et al., 2024). When a local evolution operator is recursively applied to its own predictions, local errors enter subsequent inputs and accumulate throughout the rollout. Rollout-based training reduces this train–test mismatch by exposing the model to self-generated states. Message Passing Neural PDE Solvers combine temporal bundling with push-forward training (Brandstetter et al., 2022), whereas Solver-in-the-Loop couples learned corrections with differentiable numerical solvers during multi-step optimization (Um et al., 2020). Graph-based simulators similarly improve autoregressive robustness through recursive training and perturbations of training states (Sanchez-Gonzalez et al., 2020; Pfaff et al., 2021). APEBench organizes related strategies into one-step learning, supervised unrolling, diverted-chain training, and differentiable solver interaction (Koehler et al., 2024). List et al. distinguish the effect of exposure to model-induced states from that of long-term gradients (List et al., 2025), while recurrent neural operators directly train on recursively generated predictions (Ye et al., 2025). More generally, scheduled sampling and Professor Forcing reduce the teacher-forcing/free-running mismatch in sequence modeling, while DAgger aggregates on-policy data to address covariate shift in imitation learning (Bengio et al., 2015; Lamb et al., 2016; Ross et al., 2011). Curriculum learning and example reweighting also adapt training signals, but do not construct a relative target from an explicitly selected historical rollout (Bengio et al., 2009; Ren et al., 2018). These methods align training more closely with autoregressive inference, but their supervision predominantly measures the absolute discrepancy between the current rollout and the ground-truth trajectory. Such objectives penalize the current prediction error without comparing it against long-horizon failure behaviors produced earlier during optimization. HERO complements this absolute supervision with a relative objective constructed from dynamically selected failure trajectories in the model’s optimization history. A.2. Long-Horizon Structural and Statistical Stability A complementary line of research controls long-horizon degradation through spectral modeling, architectural constraints, and structure-preserving objectives. McCabe et al. connect autoregressive instability to aliasing and spectral sensitivity (McCabe et al., 2023), while PDE-Refiner improves the representation of low-amplitude frequency components through iterative refinement (Lippe et al., 2023). The Spectral Generator Neural Operator constrains repeatedly propagated spectral dynamics (Li et al., 2026), and hierarchical implicit neural emulators introduce multiscale representations of future states to improve temporal coherence (Jiang et al., 2025). Geometric and physical constraints further encode symmetry and conservation structure into neural PDE surrogates (Huang and Greenberg, 2025). Physics-informed formulations impose governing-equation residuals, whereas structure-preserving approaches can encode geometric invariants in learned dynamical systems (Raissi et al., 2019; Karniadakis et al., 2021; Zhong et al., 2020). For dissipative chaotic systems, where long-term pointwise tracking is limited by sensitivity to initial conditions, prior work instead targets dissipative behavior and long-term statistics. Existing methods constrain learned dynamics on the attractor (Li et al., 2022), preserve invariant measures through optimal-transport and contrastive objectives (Jiang et al., 2023), or introduce scalable measure-matching regularization (Schiff et al., 2024). At the broader application scale, fluid simulation and data-driven weather forecasting motivate evaluating multi-step trajectories and forecast lead times alongside one-step error (Brunton et al., 2020; Kochkov et al., 2021; Rasp et al., 2020). These approaches improve long-horizon behavior through operator design, iterative refinement, auxiliary representations, or predefined physical and statistical constraints. HERO addresses a complementary supervision problem: it leaves the inference-time operator unchanged and uses training-only references from the model’s optimization history to measure whether the current rollout has improved over previously observed failure behavior. Appendix B Method B.1. Gradient Flow Analysis HERO separates reference construction from parameter optimization. Within the relative branch, the current rollout τθs _ _s is the only differentiable trajectory. All active candidate trajectories, together with the diagnostic values, candidate scores S(j)S^(j), selected failure reference τ−τ^-, and reference error dfaild_fail, are detached from the computation graph. Consequently, gradients do not propagate through the lagged operator, perturbed branch, diagnostic computation, or candidate-selection procedure. For the sample-level relative loss, the detached reference satisfies ∇θsdfail=0 _ _sd_fail=0. Its gradient is therefore (28) ∇θsℓrel=σ(β[dcur−dfail+m])∇θsdcur, _ _s _rel=σ (β [d_cur-d_fail+m ] ) _ _sd_cur, where σ(x)=1/(1+exp(−x))σ(x)=1/(1+ (-x)) is the logistic sigmoid. The relative objective therefore does not optimize the selected failure trajectory or any candidate branch. It scales the gradient that moves the current rollout toward the ground-truth trajectory by a factor in (0,1)(0,1). This factor approaches one when the margin is strongly violated and decreases as the current rollout outperforms the selected reference. The trajectory diagnostics determine which candidate is used for comparison, but do not introduce separate gradients into the parameter update. Algorithm 1 History-Enriched Relative Supervision in HERO 1:Current operator GθsG_ _s, lagged parameters θ¯s−1 θ_s-1, current rollout τθs _ _s, input context znz_n, ground-truth trajectory τ⋆ _ , and active diagnostics D 2:Relative loss ℓrel _rel and updated lagged parameters θ¯s θ_s 3:Stage 1: Update the lagged operator 4:Refresh or retain θ¯s θ_s according to Eq. (9). 5:Stage 2: Construct candidate rollouts 6:Construct the detached self-reference and the lagged rollout according to Eqs. (11) and (10). 7:Add the perturbed rollout defined in Eq. (12) when the perturbation branch is enabled. 8:Form the active candidate set C and enumerate its members as τ(1),…,τ(M)τ^(1),…,τ^(M). 9:Stage 3: Evaluate candidate trajectories 10:for j=1,…,Mj=1,…,M do 11: Compute the active rollout, spectral, energy, and error-growth diagnostics defined in Eqs. (14)–(18). 12:end for 13:Stage 4: Score candidates and select a failure reference 14:Normalize each diagnostic across the active candidates using Eq. (20). 15:Compute the aggregate candidate scores using Eq. (20). 16:Select the highest-scoring candidate and set it as the failure reference τ−τ^-. 17:Stage 5: Construct the relative objective 18:Compute the current and reference errors relative to τ⋆ _ . 19:Compute the sample-level relative loss using Eq. (23). 20:return ℓrel _rel and θ¯s θ_s B.2. Theoretical Analysis Every historical component of HERO is detached from the computational graph: candidate trajectories, diagnostic scores, and reference-selection operations are all evaluated under stop-gradient, so model parameters receive gradients only through the current rollout. The historical reference therefore supplies no backpropagation direction of its own, and its influence on training enters through a different channel. That channel is a task-dependent activation weight, which rescales each task gradient by a bounded factor determined by how the current rollout compares with its selected reference. HERO thus modifies the relative contributions of different tasks to the aggregated update direction rather than the direction of any individual task gradient. This reweighting view organizes the analysis into three steps: i) we establish the task-level gradient reweighting form of HERO and bound its modulation factor, i) we decompose the first-order effect of the reweighting on the long-horizon risk into a uniform scaling term and a covariance correction, and i) we convert a uniform long-horizon descent margin into a finite-step entry guarantee for a prescribed long-horizon sublevel set. B.2.1. Long-Horizon Trajectory Risk Let S denote the physical task space, equipped with a probability measure ρ. A task q∈q may encode an initial condition, equation parameter, external forcing, boundary condition, or geometric configuration. Let θ∈ℝpθ ^p denote the model parameters, and let X be a real Hilbert space representing the physical state space. For a task q, the model-generated trajectory over the time interval [0,T][0,T] is denoted by u^θq∈C([0,T];X) u_θ^q∈ C ([0,T];X ), and the corresponding reference trajectory is denoted by uq∈C([0,T];X)u^q∈ C ([0,T];X ). Given a continuous nonnegative state discrepancy d:X×X→ℝ≥0d:X× X _≥ 0, the trajectory loss for task q over the temporal horizon T is defined as (29) ℓT(q,θ)=1T∫0Td(u^θq(t),uq(t))dt. _T(q,θ)= 1T _0^Td\! ( u_θ^q(t),u^q(t) )\,dt. The discrepancy d is not required to satisfy all axioms of a mathematical metric. For instance, when d is a normalized relative error, ℓT _T should be interpreted as a trajectory discrepancy functional. Let TtrT_tr denote the rollout horizon used during training and let TLT_L denote the long-horizon evaluation window, with Ttr≤TLT_tr≤ T_L. The task-level training loss is ℓtr(q,θ):=ℓTtr(q,θ) _tr(q,θ):= _T_tr(q,θ), whereas the population-level long-horizon risk is (30) JL(θ)=∫ℓTL(q,θ)ρ(dq).J_L(θ)= _S _T_L(q,θ)\,ρ(dq). For a prescribed tolerance ϵ>0ε>0, define the long-horizon sublevel set (31) ϵ=θ∈ℝp:JL(θ)≤ϵ.G_ε= \θ ^p:J_L(θ)≤ε \. The set ϵG_ε contains all model parameters satisfying the prescribed long-horizon performance criterion. It is not required to be convex, connected, bounded, or associated with a unique minimizer. Throughout the analysis, JLJ_L and ℓtr(q,⋅) _tr(q,·) are differentiable, and all task-level gradients appearing under the integral with respect to ρ are integrable. B.2.2. Historical-Reference-Induced Relative Activation At iteration k, HERO constructs a collection of candidate trajectories from the current operator, a periodically lagged operator, and locally perturbed rollouts. A diagnostic and selection procedure identifies a task-dependent reference trajectory, denoted by u~kq u_k^q. The selected reference trajectory is treated as a constant during the current parameter update. Its trajectory loss is (32) bk(q)=1Ttr∫0Ttrd(u~kq(t),uq(t))dt,b_k(q)= 1T_tr _0^T_trd\! ( u_k^q(t),u^q(t) )\,dt, with (33) ∇θbk(q)=0. _θb_k(q)=0. For a prescribed relative margin m>0m>0, define (34) hk(q,θ)=ℓtr(q,θ)−bk(q)+m.h_k(q,θ)= _tr(q,θ)-b_k(q)+m. The quantity hk(q,θ)h_k(q,θ) is positive when the current rollout has not yet achieved an improvement of at least m relative to the selected reference. The smooth relative trajectory loss is defined by (35) ℓrel,k(q,θ)=1βlog(1+exp(βhk(q,θ))), _rel,k(q,θ)= 1β \! (1+ \! (β h_k(q,θ) ) ), where β>0β>0 controls the sharpness of the transition around the relative-margin boundary. Using (33), differentiation of (35) gives (36) ∇θℓrel,k(q,θ)=αk(q,θ)∇θℓtr(q,θ), _θ _rel,k(q,θ)= _k(q,θ) _θ _tr(q,θ), where (37) αk(q,θ)=σ(βhk(q,θ))∈(0,1), _k(q,θ)=σ\! (β h_k(q,θ) )∈(0,1), and σ(z)=(1+exp(−z))−1σ(z)= (1+ (-z) )^-1 is the sigmoid function. The coefficient αk(q,θ) _k(q,θ) represents the relative activation level of task q. When the current trajectory fails to satisfy the required relative margin, αk _k approaches one. When the current trajectory already substantially outperforms the selected reference, αk _k approaches zero. Near the relative-margin boundary, αk _k provides a smooth interpolation between these regimes. Let λk≥0 _k≥ 0 denote the relative-loss coefficient at iteration k. The task-level HERO objective is (38) ℓH,k(q,θ)=ℓtr(q,θ)+λkℓrel,k(q,θ). _H,k(q,θ)= _tr(q,θ)+ _k _rel,k(q,θ). Define the task-level training gradient as g(q,θ)=∇θℓtr(q,θ)g(q,θ)= _θ _tr(q,θ). Combining (36) and (38) yields (39) ∇θℓH,k(q,θ)=wk(q,θ)g(q,θ), _θ _H,k(q,θ)=w_k(q,θ)g(q,θ), where (40) wk(q,θ)=1+λkαk(q,θ).w_k(q,θ)=1+ _k _k(q,θ). Since αk(q,θ)∈(0,1) _k(q,θ)∈(0,1) and λk≥0 _k≥ 0, the task weight is bounded as 1≤wk(q,θ)≤1+λk1≤ w_k(q,θ)≤ 1+ _k, where the inequalities are strict for λk>0 _k>0. Equation (39) shows that the historical reference does not directly introduce an additional task-level backpropagation direction. Instead, it induces a bounded task-dependent modulation of the current rollout gradient. Because different tasks generally produce different activation levels, HERO changes their relative contributions to the aggregated update. The exact scalar-reweighting form in (39) applies when the base rollout objective and the current-trajectory term in the relative objective are both given by ℓtr _tr. If a distinct base objective ℒbaseL_base is used, the task-level gradient takes the form ∇θℒbase+λkαk∇θℓtr _θL_base+ _k _k _θ _tr, which is not, in general, a scalar multiple of a single task gradient. B.2.3. Aggregated HERO Update Direction The population negative-gradient direction associated with standard trajectory training is (41) d0(θ)=−∫g(q,θ)ρ(dq).d^0(θ)=- _Sg(q,θ)\,ρ(dq). The corresponding HERO direction at iteration k is (42) dkH(θ)=−∫wk(q,θ)g(q,θ)ρ(dq).d_k^H(θ)=- _Sw_k(q,θ)g(q,θ)\,ρ(dq). The parameter update is (43) θk+1=θk+ηdkH(θk), _k+1= _k+η d_k^H( _k), where η>0η>0 is the learning rate. To characterize whether a task gradient contributes to reducing the long-horizon risk, define the long-horizon usefulness of task q as (44) a(q,θ)=⟨∇JL(θ),g(q,θ)⟩.a(q,θ)= ∇ J_L(θ),g(q,θ) . If a(q,θ)>0a(q,θ)>0, then the task-level negative-gradient direction −g(q,θ)-g(q,θ) is a first-order descent direction for JLJ_L. Define the mean activation (45) α¯k(θ)=∫αk(q,θ)ρ(dq), α_k(θ)= _S _k(q,θ)\,ρ(dq), the mean task usefulness (46) a¯(θ)=∫a(q,θ)ρ(dq), a(θ)= _Sa(q,θ)\,ρ(dq), and the covariance between relative activation and long-horizon usefulness: Covρ(αk,a)= _ρ( _k,a)= ∫(αk(q,θ)−α¯k(θ)) _S ( _k(q,θ)- α_k(θ) ) (47) ×(a(q,θ)−a¯(θ))ρ(dq). × (a(q,θ)- a(θ) )\,ρ(dq). This covariance quantifies whether HERO assigns larger activation weights to tasks whose training gradients provide stronger first-order descent contributions to the long-horizon risk. Proposition B.1 (Historical-reference-induced correction of the long-horizon descent component). For any parameter vector θ and iteration k, the HERO direction satisfies (48) ⟨∇JL(θ),dkH(θ)⟩ ∇ J_L(θ),d_k^H(θ) =(1+λkα¯k(θ))⟨∇JL(θ),d0(θ)⟩ = (1+ _k α_k(θ) ) ∇ J_L(θ),d^0(θ) −λkCovρ(αk,a). - _kCov_ρ( _k,a). Proof. Using (42), (40), and (44), (49) ⟨∇JL,dkH⟩ ∇ J_L,d_k^H =−∫(1+λkαk)aρ(dq) =- _S (1+ _k _k )a\,ρ(dq) =−a¯−λk∫αkaρ(dq). =- a- _k _S _ka\,ρ(dq). The covariance identity gives (50) ∫αkaρ(dq)=α¯ka¯+Covρ(αk,a). _S _ka\,ρ(dq)= α_k a+Cov_ρ( _k,a). Moreover, from (41), (51) ⟨∇JL,d0⟩=−a¯. ∇ J_L,d^0 =- a. Substituting (50) and (51) into (49) yields (48). ∎ Proposition B.1 separates the first-order long-horizon effect of HERO into two components. The first term, (1+λkα¯k)⟨∇JL,d0⟩ (1+ _k α_k ) ∇ J_L,d^0 , corresponds to a uniform scaling of the standard training direction. The second term, −λkCovρ(αk,a)- _kCov_ρ( _k,a), is the correction induced by nonuniform task reweighting. When Covρ(αk,a)>0Cov_ρ( _k,a)>0, tasks with larger long-horizon usefulness receive, on average, larger activation weights. The covariance term then contributes an additional negative component to ⟨∇JL,dkH⟩ ∇ J_L,d_k^H , thereby increasing the first-order descent contribution of the aggregated HERO update relative to uniform gradient scaling. To separate the nonuniform reweighting effect from the global magnitude change, define the scale-normalized HERO direction (52) d~kH(θ)=dkH(θ)1+λkα¯k(θ). d_k^H(θ)= d_k^H(θ)1+ _k α_k(θ). Since the denominator is strictly positive, d~kH d_k^H and dkHd_k^H have the same geometric direction. Proposition B.1 implies (53) ⟨∇JL,d~kH⟩ ∇ J_L, d_k^H =⟨∇JL,d0⟩ = ∇ J_L,d^0 −λk1+λkα¯kCovρ(αk,a). - _k1+ _k α_kCov_ρ( _k,a). Equation (53) isolates the first-order contribution produced exclusively by task-dependent reweighting. A positive covariance is not guaranteed solely by the form of the relative objective. The lagged reference, perturbed candidate construction, and diagnostic selection mechanism determine whether the induced activation weights are positively associated with task-level long-horizon usefulness. Corollary B.2 (Preservation and enhancement of a standard first-order descent trend). Suppose that, in a parameter region of interest, the standard training direction satisfies (54) ⟨∇JL(θ),d0(θ)⟩≤−γ0,γ0>0, ∇ J_L(θ),d^0(θ) ≤- _0, _0>0, and (55) Covρ(αk,a)≥ck,ck≥0.Cov_ρ( _k,a)≥ c_k, c_k≥ 0. Then (56) ⟨∇JL,dkH⟩ ∇ J_L,d_k^H ≤−(1+λkα¯k)γ0−λkck ≤- (1+ _k α_k ) _0- _kc_k ≤−γ0−λkck. ≤- _0- _kc_k. Proof. The result follows directly by substituting (54) and (55) into (48). ∎ Corollary B.2 shows that a nonnegative covariance preserves an existing first-order long-horizon descent trend. A strictly positive covariance provides an additional first-order descent margin. Corollary B.3 (Correction of a locally adverse standard direction). Suppose that, in a parameter region of interest, (57) ⟨∇JL(θ),d0(θ)⟩≤Bk,Bk≥0, ∇ J_L(θ),d^0(θ) ≤ B_k, B_k≥ 0, and Covρ(αk,a)≥ckCov_ρ( _k,a)≥ c_k. Then (58) ⟨∇JL,dkH⟩≤(1+λkα¯k)Bk−λkck. ∇ J_L,d_k^H ≤ (1+ _k α_k )B_k- _kc_k. If λk>0 _k>0 and (59) λkck>(1+λkα¯k)Bk, _kc_k> (1+ _k α_k )B_k, then ⟨∇JL,dkH⟩<0 ∇ J_L,d_k^H <0. Proof. Equation (58) follows from Proposition B.1. Condition (59) makes the right-hand side of (58) strictly negative. ∎ Since 0<α¯k<10< α_k<1, a more conservative sufficient condition is (60) λkck>(1+λk)Bk. _kc_k>(1+ _k)B_k. Corollary B.3 shows that even when the standard short-horizon training direction has a locally adverse first-order effect on the long-horizon risk, sufficiently strong positive covariance can make the aggregated HERO direction a long-horizon descent direction. B.2.4. Mini-Batch Form Consider a mini-batch ℬ=q1,…,qBB= \q_1,…,q_B \. For each task, define (61) gi=g(qi,θ),αi=αk(qi,θ),ai=⟨∇JL(θ),gi⟩.g_i=g(q_i,θ), _i= _k(q_i,θ), a_i= ∇ J_L(θ),g_i . The empirical standard and HERO directions are (62) d^0=−1B∑i=1Bgi,d^kH=−1B∑i=1B(1+λkαi)gi. d^0=- 1B _i=1^Bg_i, d_k^H=- 1B _i=1^B (1+ _k _i )g_i. Define the empirical means (63) α¯ℬ=1B∑i=1Bαi,a¯ℬ=1B∑i=1Bai, α_B= 1B _i=1^B _i, a_B= 1B _i=1^Ba_i, and the empirical covariance (64) Cov^ℬ(α,a)=1B∑i=1B(αi−α¯ℬ)(ai−a¯ℬ). Cov_B(α,a)= 1B _i=1^B ( _i- α_B ) (a_i- a_B ). The same algebraic decomposition gives (65) ⟨∇JL,d^kH⟩ ∇ J_L, d_k^H =(1+λkα¯ℬ)⟨∇JL,d^0⟩ = (1+ _k α_B ) ∇ J_L, d^0 −λkCov^ℬ(α,a). - _k Cov_B(α,a). Thus, the population-level mechanism in Proposition B.1 also applies directly to finite mini-batches: nonuniform relative activation changes the aggregated gradient through an empirical covariance correction. B.2.5. Finite-Step Entry into the Long-Horizon Sublevel Set The preceding analysis characterizes how historical-reference-induced reweighting modifies the first-order long-horizon descent component. The following result establishes finite-step entry into ϵG_ε when this descent contribution remains uniformly nondegenerate outside the target set. Theorem B.4 (Finite-step entry under uniform long-horizon descent). Assume that ϵ≠∅G_ε≠ , and consider the HERO iteration (66) θk+1=θk+ηdkH(θk). _k+1= _k+η d_k^H( _k). Before the first entry into ϵG_ε, suppose that the following conditions hold: (1) The gradient of JLJ_L is L-Lipschitz continuous along every update segment θk+tηdkH(θk):t∈[0,1] \ _k+tη d_k^H( _k):t∈[0,1] \. (2) There exists γϵ>0 _ε>0 such that (67) ⟨∇JL(θk),dkH(θk)⟩≤−γϵ. ∇ J_L( _k),d_k^H( _k) ≤- _ε. (3) There exists Dϵ<∞D_ε<∞ such that ‖dkH(θk)‖2≤Dϵ \|d_k^H( _k) \|_2≤ D_ε. If the learning rate satisfies (68) 0<η≤γϵLDϵ2,0<η≤ _εLD_ε^2, then the parameter sequence enters ϵG_ε after finitely many iterations. Define the first-entry time τH(ϵ)=infk≥0:θk∈ϵ _H(ε)= \k≥ 0: _k _ε \. Then (69) τH(ϵ)≤⌈2[JL(θ0)−ϵ]+ηγϵ⌉, _H(ε)≤ 2 [J_L( _0)-ε ]_+η _ε , where [x]+=maxx,0[x]_+= \x,0 \. Proof. Suppose that θk∉ϵ _k _ε. By the L-smoothness of JLJ_L along the update segment, (70) JL(θk+1) J_L( _k+1) ≤JL(θk)+⟨∇JL(θk),θk+1−θk⟩ ≤ J_L( _k)+ ∇ J_L( _k), _k+1- _k +L2‖θk+1−θk‖22. + L2 \| _k+1- _k \|_2^2. Using (66), θk+1−θk=ηdkH(θk) _k+1- _k=η d_k^H( _k). Therefore, (71) JL(θk+1) J_L( _k+1) ≤JL(θk)+η⟨∇JL(θk),dkH(θk)⟩ ≤ J_L( _k)+η ∇ J_L( _k),d_k^H( _k) +Lη22‖dkH(θk)‖22. + Lη^22 \|d_k^H( _k) \|_2^2. Applying (67) and the bound ‖dkH(θk)‖2≤Dϵ \|d_k^H( _k) \|_2≤ D_ε yields (72) JL(θk+1)≤JL(θk)−ηγϵ+Lη2Dϵ22.J_L( _k+1)≤ J_L( _k)-η _ε+ Lη^2D_ε^22. By (68), Lη2Dϵ2/2≤ηγϵ/2Lη^2D_ε^2/2≤η _ε/2. Hence, as long as θk∉ϵ _k _ε, (73) JL(θk+1)≤JL(θk)−ηγϵ2.J_L( _k+1)≤ J_L( _k)- η _ε2. Assume that the parameter sequence remains outside ϵG_ε for the first N iterations. Summing (73) gives (74) JL(θN)≤JL(θ0)−Nηγϵ2.J_L( _N)≤ J_L( _0)- Nη _ε2. For (75) N≥2[JL(θ0)−ϵ]+ηγϵ,N≥ 2 [J_L( _0)-ε ]_+η _ε, equation (74) implies JL(θN)≤ϵJ_L( _N)≤ε, and therefore θN∈ϵ _N _ε. This contradicts the assumption that the sequence remains outside ϵG_ε. Consequently, the first-entry time is finite and satisfies (69). ∎ B.2.6. Scope of the Theoretical Results Theorem B.4 is a conditional finite-step entry result rather than an unconditional convergence guarantee. First, the theorem establishes the first entry of the parameter sequence into ϵG_ε. It does not establish that ϵG_ε is invariant under subsequent updates, that the sequence remains permanently inside the set, or that the parameters converge to a unique point. Such conclusions require additional conditions inside the target region. Second, Proposition B.1 characterizes the method-specific effect of historical-reference-induced task reweighting. Theorem B.4, by contrast, applies a standard smoothness-based descent argument once a uniform long-horizon descent margin has been established. Third, a positive covariance alone does not imply the existence of a uniform constant γϵ>0 _ε>0. Finite-step entry additionally requires that the long-horizon descent contribution remain nondegenerate before reaching ϵG_ε. Fourth, the analysis is expressed in terms of population gradients or empirical mini-batch gradients. Practical optimization may additionally involve random candidate perturbations, finite-batch noise, momentum, and adaptive preconditioning. The theoretical results therefore characterize the idealized gradient-level mechanism induced by historical rollout references. The covariance correction is nonetheless an empirically accessible quantity rather than a purely formal device, since ∇JL∇ J_L can be approximated by the mean long-horizon task gradient on a finite validation task set, and the association between relative activation and long-horizon usefulness can then be measured directly through Corr(αi,ai)Corr( _i,a_i). Taken together, the analysis identifies bounded task reweighting, rather than any additional backpropagation direction, as the channel through which historical rollout references act on long-horizon training. Appendix C Experiments C.1. Experimental Setup C.1.1. Baselines and Backbones Backbones. To assess whether HERO is architecture agnostic, we instantiate every method on two structurally different neural operators. FNO is the standard spectral operator for PDE surrogate learning: it parameterizes the integral kernel in the Fourier domain and applies spectral convolutions that retain only the low-frequency Fourier modes, modeling global spatial dependence directly in frequency space. Transolver is a Transformer-based neural solver whose Physics-Attention adaptively partitions the discrete field into learnable slices and models interactions among these latent physical states rather than among all grid points, providing a global but non-spectral inductive bias. The two backbones thus probe complementary regimes—frequency-domain spectral convolution versus slice-level attention—and let us test whether HERO transfers from spectral neural operators to attention-based ones. Unless a suffix is appended, FNO and Transolver denote the base backbones trained with one-step supervision and no rollout-stabilization mechanism. Comparison methods. All competing training strategies are applied on both backbones, yielding matched backbone–strategy pairs so that any difference is attributable to the training strategy rather than to the architecture. Beyond the one-step base models, we compare against three rollout-stabilization strategies. Push-forward training (PF) unrolls the model on its own predictions during training to close the gap between one-step teacher forcing and test-time free rollout: for a backbone fθf_θ initialized at u^t=ut u_t=u_t, it forms a KpfK_pf-step self-rollout u^t+s=fθ(u^t+s−1) u_t+s=f_θ( u_t+s-1) for s=1,…,Kpfs=1,…,K_pf and minimizes (76) ℒpf=1Kpf∑s=1Kpf‖u^t+s−ut+s‖22,L_pf= 1K_pf _s=1^K_pf \|\, u_t+s-u_t+s\, \|_2^2, where u^t+s u_t+s and ut+su_t+s are the predicted and ground-truth states at rollout offset s. PDE-Refiner (Refine) instead targets the prediction itself, applying a multi-step refinement/denoising process that restores non-dominant spatial-frequency content lost during rollout, and thus tests whether output-side iterative correction alone stabilizes long-horizon prediction. RNO is a recurrent-training framework that recursively applies the operator to its own predictions so that training matches the autoregressive test-time dynamics, reducing exposure bias and long-horizon error amplification. We denote the resulting variants by Backbone-Strategy (e.g., FNO-PF, Transolver-Refine, FNO-RNO), and the HERO variants analogously by FNO-HERO and Transolver-HERO. Relation to HERO. PF and RNO are closest to HERO, as all three expose the operator to self-generated states during training, yet both supervise only the absolute discrepancy from the ground truth. Refine differs along another axis, correcting the current output rather than using optimization history. In contrast, HERO keeps ground-truth rollout regression but adds relative supervision against a dynamically history-enriched failure reference, rewarding improvement over earlier failure behavior rather than mere closeness to the ground truth. C.1.2. Benchmarks We evaluate HERO on nine time-dependent PDE benchmarks, spanning 1D, 2D, and 3D dynamics and covering dispersion, nonlinear transport, dissipative chaos, forced and decaying turbulence, pattern formation, and multi-directional advection. All benchmarks use periodic boundary conditions and a single scalar state channel, and each dataset provides 5050 training and 3030 test trajectories; long-horizon behavior is assessed with 200200-step free autoregressive rollouts. The spatial dimension, domain, resolution, and time step of every benchmark are summarized in Table 4. The 1D benchmarks isolate phase- and spectrum-dominated failure modes. The Dispersion equation governs waves whose wavenumber components propagate at different phase speeds, governed by ∂tu=α3∂xxxu _tu= _3\, _xu with α3=2.5×10−4 _3=2.5× 10^-4, where α3 _3 is the dispersion coefficient; this makes it a direct test of long-horizon phase consistency and frequency-domain structure. The Burgers equation couples nonlinear convection with viscous diffusion, governed by ∂tu=a2∂xxu+b1u∂xu _tu=a_2\, _xu+b_1\,u\, _xu with a2=3.0×10−4a_2=3.0× 10^-4 and b1=−0.125b_1=-0.125, where a2a_2 is the diffusion coefficient and b1b_1 the convection coefficient; it probes long-horizon error accumulation under shock formation and dissipation. The KdV equation adds high-order stabilization to nonlinear dispersive transport, governed by ∂tu=δu∂xu+γ3∂xxxu+γ4∂xxxxu _tu=δ\,u\, _xu+ _3\, _xu+ _4\, _xu with δ=−2.0δ=-2.0, γ3=−14.0 _3=-14.0, and γ4=−9.0 _4=-9.0, where δ is the nonlinear coefficient and γ3,γ4 _3, _4 the third- and fourth-order coefficients; it stresses long-horizon phase error, amplitude drift, and solitary-wave structure preservation. The Kuramoto–Sivashinsky (KS) equation is a canonical dissipative-chaotic system, (77) ∂tu=a2∂xxu+a4∂xxxxu+b2(∂xu)2, _tu=a_2\, _xu+a_4\, _xu+b_2\,( _xu)^2, where a2=−1.0a_2=-1.0 drives long-wavelength instability, a4=−1.0a_4=-1.0 provides high-order dissipation, and b2=−1.0b_2=-1.0 scales the gradient nonlinearity; its chaotic dynamics make it a stringent test of long-horizon error amplification and rollout stability. The 2D benchmarks emphasize energy- and structure-dominated behavior on diffusion and vorticity fields. Anisotropic Diffusion evolves a scalar field with direction-dependent and cross-directional diffusion strengths, (78) ∂tu=∇⋅(∇u),=(1.0×10−35.0×10−45.0×10−42.0×10−3), _tu=∇\!·(D\,∇ u), = pmatrix1.0× 10^-3&5.0× 10^-4\\[2.0pt] 5.0× 10^-4&2.0× 10^-3 pmatrix, where D is the symmetric diffusion tensor; it tests whether the operator preserves directionally coupled smoothing and energy decay. Kolmogorov Flow models forced 2D incompressible vorticity dynamics, (79) ∂tω+⋅∇ω=νΔω−αω+f(), _tω+v\!·\!∇ω=ν\, ω-α\,ω+f(x), where ω is the vorticity, v the associated velocity field, ν=1.0×10−2ν=1.0× 10^-2 the viscosity, α the linear drag, and f a fixed periodic forcing (injection mode 44, scale 1.01.0); the sustained energy injection and dissipation make it a test of long-horizon spectral balance. Decaying Turbulence is the unforced incompressible Navier–Stokes system in vorticity form, governed by ∂tω+⋅∇ω=νΔω _tω+v\!·\!∇ω=ν\, ω with ∇⋅=0∇\!·\!v=0 and ν=1.0×10−4ν=1.0× 10^-4, where ν is the viscosity; it evaluates vortex transport, energy dissipation, and the stability of a freely decaying field. The 3D benchmarks stress spatial expressivity and stability in higher dimensions. Swift–Hohenberg is a reaction–diffusion pattern-formation model, (80) ∂tu=ru−(qc2+Δ)2u+p(u),p(u)=u2−u3, _tu=r\,u-(q_c^2+ )^2u+p(u), p(u)=u^2-u^3, where r=0.7r=0.7 is the reactivity, qc=1.0q_c=1.0 the critical wavenumber, and p(u)p(u) the local reaction polynomial; the emergent 3D patterns demand both spatial expressivity and long-horizon stability. Unbalanced Advection transports a scalar field with distinct per-axis velocities, governed by ∂tu+⋅∇u=0 _tu+c\!·\!∇ u=0 with =(0.01,−0.04,0.005)c=(0.01,-0.04,0.005), where c is the constant advection velocity; the differing magnitudes and signs across axes induce direction-dependent phase drift over long rollouts. Table 4. Configuration of the nine PDE rollout benchmarks Benchmark Dim Domain Resolution Δt t Dispersion 1D [0,1)[0,1) 160160 1.0×10−31.0× 10^-3 Burgers 1D [0,1)[0,1) 160160 1.0×10−11.0× 10^-1 KdV 1D [0,2π)[0,2π) 160160 1.0×10−21.0× 10^-2 KS 1D [0,60)[0,60) 160160 1.0×10−11.0× 10^-1 Anisotropic Diff. 2D [0,1)2[0,1)^2 64264^2 1.0×10−11.0× 10^-1 Kolmogorov Flow 2D [0,2π)2[0,2π)^2 64264^2 1.0×10−11.0× 10^-1 Decaying Turb. (NS) 2D [0,1)2[0,1)^2 64264^2 1.0×10−11.0× 10^-1 Swift–Hohenberg 3D [0,10π)3[0,10π)^3 32332^3 1.0×10−11.0× 10^-1 Unbalanced Adv. 3D [0,1)3[0,1)^3 32332^3 1.0×10−11.0× 10^-1 • Δt t denotes the time step. Resolution denotes NxN_x (1D), Nx×NyN_x\!×\!N_y (2D), or Nx×Ny×NzN_x\!×\!N_y\!×\!N_z (3D). C.1.3. Evaluation Metrics We evaluate all models under a common rollout protocol and report metrics that jointly capture one-step fidelity, long-horizon accuracy, and rollout stability over the first 100100 autoregressive steps. All metrics build on the same per-step relative state error: for evaluation trajectory j and rollout step t, the per-step normalized RMSE is nRMSEt(j)=‖u^t(j)−ut(j)‖2/(‖ut(j)‖2+ε)nRMSE_t^(j)=\| u_t^(j)-u_t^(j)\|_2/(\|u_t^(j)\|_2+ ), the per-trajectory form of the mean relative error used in the motivation study (Eq. (1)), where u^t(j) u_t^(j) and ut(j)u_t^(j) are the predicted and ground-truth states and ε>0 >0 is a numerical-stability constant. Field-level accuracy is summarized by the set-averaged per-step error nRMSE@T=1M∑j=1MnRMSET(j)nRMSE@T= 1M _j=1^MnRMSE_T^(j), which we report at two steps. nRMSE@1nRMSE@1 isolates one-step prediction fidelity, whereas nRMSE@100nRMSE@100 measures accuracy after 100100 autoregressive steps; a small nRMSE@1nRMSE@1 together with a large nRMSE@100nRMSE@100 therefore directly exposes autoregressive error accumulation, phase drift, and high-frequency error growth, and lower values are better for both. A single late-step value does not reflect the whole rollout, and autoregressive errors tend to grow multiplicatively. We therefore summarize the entire 100100-step trajectory with a geometric-mean error. For trajectory j, (81) GM100(j)=exp(1100∑t=1100log(nRMSEt(j)+ε)),GM100^(j)= \! ( 1100 _t=1^100 \! (nRMSE_t^(j)+ ) ), where the logarithmic averaging measures relative error growth uniformly across the horizon rather than letting the largest late-step errors dominate. We report its split-level mean, ValmeanGM100=1Mval∑jGM100(j)Val\ mean\ GM100= 1M_val _jGM100^(j) and TestmeanGM100=1Mtest∑jGM100(j)Test\ mean\ GM100= 1M_test _jGM100^(j), where the validation value drives model selection and early stopping and the test value is the reported long-horizon accuracy; lower is better, and a low TestmeanGM100Test\ mean\ GM100 indicates that the model stays accurate throughout the rollout rather than only at isolated steps. Error magnitude alone does not indicate when a rollout destabilizes, so we additionally report the horizon over which a prediction remains usable. With a relative-error threshold τ=0.1τ=0.1, the stable horizon of trajectory j is (82) S(j)=maxT∈[1,100]|nRMSEt(j)≤τ∀t≤T,S^(j)= \! \\,T∈[1,100]\ |\ nRMSE_t^(j)≤τ\ \ ∀\,t≤ T\, \, where τ=0.1τ=0.1 marks a 10%10\% relative-error tolerance and S(j)S^(j) is the last step up to which the error stays within this tolerance without interruption. The Meanstablestep=1M∑j=1MS(j)Mean\ stable\ step= 1M _j=1^MS^(j) averages this horizon over the evaluation set; higher is better, and a value approaching 100100 means that most trajectories remain stable across the full evaluation horizon. C.1.4. Implementation Details Protocol. We follow the APEBench task definitions and data splits. Each task provides 5050 training trajectories of length Ttrain=51T_train=51 and 3030 test trajectories of length Ttest=201T_test=201; at test time every model is initialized from the first frame of a test trajectory and rolled out for Teval=200T_eval=200 autoregressive steps (Table 4). To keep the comparison fair, all compared methods use the same backbone, data split, optimizer, and evaluation protocol, and differ only in the training objective. The shared absolute supervision is a Kpf=5K_pf=5 push-forward rollout (Eq. (76)), so that every model is exposed to its own predicted intermediate states during training; HERO augments this with the relative objective ℒrelL_rel (Eq. (24)), giving the total objective ℒtotal=ℒroll+λsℒrelL_total=L_roll+ _sL_rel (Eq. (26)). HERO configuration. Unless stated otherwise, HERO uses a single set of defaults across all tasks: rollout horizon K=Kpf=5K=K_pf=5, start step s0=4000s_0=4000, refresh interval R=2000R=2000, warm-up length W=2000W=2000, maximum relative weight λmax=0.02 _ =0.02, margin m=0.02m=0.02, and sharpness β=5.0β=5.0, together with the default candidate set =τlag,τpert,τselfC=\ _lag, _pert, _self\. Relative supervision is enabled after s0s_0 steps, and the relative weight λs _s follows a linear warm-up from 0 to λmax _ over the following W steps. The candidate score is the within-candidate mean of the normalized diagnostics, without any additional diagnostic weighting; when a diagnostic does not apply to a given data form, it is removed from the active diagnostic set D and the remaining diagnostics are averaged. The lagged operator is initialized at step s0s_0 and refreshed every R steps, and all candidate trajectories and diagnostic scores are computed without gradients. Because reference construction, candidate generation, and diagnostics are confined to training, HERO leaves the architecture, parameter count, and inference procedure of its backbone unchanged and introduces no additional inference-time modules or parameters. Optimization. We use the AdamW optimizer with an initial learning rate of 10−310^-3 and a warm-up–cosine schedule with 20002000 warm-up steps, for a total of 10 00010\,000 steps at a batch size of 2020. All methods share the same random seed, train/validation/test split, and evaluation scripts. For each task we select the checkpoint with the best long-horizon validation metric (Val mean GM100) and report its performance on the test set. Backbone configurations. We evaluate two backbones with per-dimension configurations summarized in Table 5. The Config column reports each backbone’s core setting: for FNO, the number of retained low-frequency Fourier modes (2828 in 1D, 8×88× 8 in 2D, and 8×8×88× 8× 8 in 3D); for Transolver, its Physics-Attention token mixing. Width is the channel width of FNO and the hidden dimension of Transolver, and Depth is the number of spectral layers of FNO and attention blocks of Transolver. The Heads and Slices columns describe Transolver’s Physics-Attention only and are marked “–” for FNO. All tasks of the same dimension share one configuration, and every training strategy—one-step, PF, Refine, RNO, and HERO—uses the identical backbone configuration, so that any performance difference is attributable to the training objective rather than to model capacity. Table 5. Backbone architecture configurations. Each per-dimension configuration is shared across all tasks of that dimension and across all training strategies (one-step, PF, Refine, RNO, HERO); only the training objective differs. Backbone Dim Applied tasks Config Width Depth Heads Slices FNO 1D Disp, Burgers, KdV, KS 2828 2828 11 – – FNO 2D AnisoDiff, Kolm, NS 8×88× 8 2020 11 – – FNO 3D SwiftH, UnbalAdv 8×8×88× 8× 8 44 55 – – Transolver 1D Disp, Burgers, KdV, KS Physics-Attn. 224224 22 44 6464 Transolver 2D AnisoDiff, Kolm, NS Physics-Attn. 256256 44 88 128128 Transolver 3D SwiftH, UnbalAdv Physics-Attn. 512512 33 88 256256 C.2. Complete Benchmark Results Tables 9 and 10 report the complete results on all nine PDE benchmarks (mean ± standard deviation over five seeds), grouped by spatial dimension, including the three benchmarks (Dispersion, KdV, and Kolmogorov) omitted from the main-text Table 1. HERO retains the same long-horizon advantage on all tasks. On the FNO backbone, it reduces nRMSE@100 from 0.4250.425 to 0.2460.246 on Dispersion, from 3.2173.217 to 0.9190.919 on KdV, and from 1.3791.379 to 0.6480.648 on Kolmogorov. Transolver-HERO further attains the lowest nRMSE@100 within the Transolver family on KdV and Kolmogorov, showing that the rollout improvement extends across both backbones. C.3. Hyperparameter Sensitivity We analyze the two hyperparameters that govern the history-enriched mechanism on the 2D Navier–Stokes benchmark: the relative-training start step s0s_0 and the refresh interval R of the lagged operator in (9). In each study, only the analyzed hyperparameter varies and all remaining HERO settings, the backbone, and the evaluation protocol are held fixed, with 10 00010\,000 total optimization steps. FNO-PF is included as the absolute-supervision reference. C.3.1. Relative-Training Start Step Table 6 varies s0s_0, so that the relative objective is active over the interval [s0,10 000][s_0,10\,000], as listed in the relative-training interval column. Performance follows a clear interior optimum. At s0=2000s_0=2000 the backbone has not yet formed a stable local time-stepping operator, the lagged rollout carries largely unstructured prediction error, and enabling the relative objective this early perturbs optimization: nRMSE@1 rises to 0.03390.0339 and Test GM100 stays at 1.3381.338, better than FNO-PF but clearly behind the default. At s0=4000s_0=4000 the reference is already informative while 60006000 steps of relative reweighting remain, giving the lowest Test GM100 (1.0671.067) and nRMSE@100 (4.1244.124) and the longest stable horizon (25.925.9). Increasing s0s_0 to 60006000 improves the reference further and yields the best one-step error (0.03230.0323), but the shortened relative phase leaves long-horizon behavior less corrected (1.1421.142 Test GM100). At s0=8000s_0=8000 only 20002000 steps remain, the one-step error approaches FNO-PF, and both the long-horizon error and the stable horizon degrade toward push-forward training alone. Table 6. Sensitivity to the relative-training start step s0s_0 on the 2D Navier–Stokes benchmark. The relative objective is active over [s0, 10 000][s_0,\,10\,000]; FNO-PF is the absolute-supervision reference. Setting Relative-training interval Val GM100 ↓ Test GM100 ↓ nRMSE @1 ↓ nRMSE @100 ↓ Stable step ↑ FNO-PF — 1.4501.450 1.6971.697 0.03500.0350 6.3066.306 22.722.7 HERO (s0=2000s_0=2000) 20002000–10 00010\,000 1.1631.163 1.3381.338 0.03390.0339 4.9174.917 24.324.3 HERO (s0=4000s_0=4000) 40004000–10 00010\,000 0.9370.937 1.0671.067 0.03260.0326 4.1244.124 25.925.9 HERO (s0=6000s_0=6000) 60006000–10 00010\,000 0.9980.998 1.1421.142 0.03230.0323 4.3864.386 25.225.2 HERO (s0=8000s_0=8000) 80008000–10 00010\,000 1.2871.287 1.5061.506 0.03420.0342 5.7145.714 23.223.2 C.3.2. Lagged-Operator Refresh Interval Table 7 varies R with s0=4000s_0=4000 fixed, so the lagged operator is refreshed at steps 4000,4000+R,…4000,4000+R,… until training ends, as listed in the effective refresh steps column. The same interior optimum appears. A short interval (R=500R=500, ten refreshes) keeps the lagged operator close to the current one, so the reference degenerates toward the self-reference and provides little contrast, giving the weakest HERO result (1.3521.352 Test GM100). A long interval (R=4000R=4000, two refreshes) leaves the reference increasingly stale relative to the improving operator, weakening the effective activation of the margin objective (1.2471.247). Intermediate intervals perform best, with R=1000R=1000 and R=2000R=2000 within 7%7\% of each other on Test GM100 and R=2000R=2000 attaining the lowest long-horizon error (4.1244.124) and the longest stable horizon (25.925.9). Every configuration in both tables outperforms FNO-PF on all long-horizon metrics, so the effect of HERO does not depend on a narrow hyperparameter range. The two studies degrade for the same reason: the periodic-lag mechanism requires a reference that is neither too close to the current operator nor too far behind it, which is exactly the balance between reference relevance and historical contrast that motivates the design. Table 7. Sensitivity to the refresh interval R at s0=4000s_0=4000 on the 2D Navier–Stokes benchmark. FNO-PF is the absolute-supervision reference. Setting Effective refresh steps Val GM100 ↓ Test GM100 ↓ nRMSE @1 ↓ nRMSE @100 ↓ Stable step ↑ FNO-PF — 1.4501.450 1.6971.697 0.03500.0350 6.3066.306 22.722.7 HERO (R=500R=500) 4000,4500,…,95004000,4500,…,9500 1.1761.176 1.3521.352 0.03340.0334 4.9824.982 24.124.1 HERO (R=1000R=1000) 4000,5000,…,90004000,5000,…,9000 0.9910.991 1.1391.139 0.03240.0324 4.3764.376 25.325.3 HERO (R=2000R=2000) 4000,6000,80004000,6000,8000 0.9370.937 1.0671.067 0.03260.0326 4.1244.124 25.925.9 HERO (R=4000R=4000) 4000,80004000,8000 1.0841.084 1.2471.247 0.03310.0331 4.6414.641 24.824.8 Table 8. Parameter count, inference latency, peak memory, and training throughput of the FNO variants under a shared measurement protocol Method #Params Infer (ms/step) ↓ Peak mem (MB) ↓ Train (steps/s) ↑ FNO 30,54130,541 0.5030.503 66.1366.13 374.79374.79 FNO-PF 30,54130,541 0.5050.505 66.2866.28 352.61352.61 FNO-Ref. 30,91730,917 0.5530.553 69.8469.84 328.47328.47 FNO-RNO 30,54130,541 0.5070.507 66.9466.94 344.26344.26 FNO-HERO 30,54130,541 0.5090.509 67.1867.18 336.92336.92 C.4. Computational Cost HERO adds candidate generation, diagnostics, and reference selection to the training loop, so we quantify what this costs in practice. Table 8 compares all FNO variants under one protocol: the same 2D backbone, batch size, resolution, and hardware, with inference latency and peak memory measured over 200200-step free autoregressive rollouts and training throughput averaged over 10001000 optimization steps after warm-up. HERO leaves the deployed model untouched. It shares the exact parameter count of the plain backbone (30,54130,541), whereas Refine adds 376376 parameters for its refinement head, and its inference latency matches the backbone to within 1.2%1.2\% (0.5090.509 vs. 0.5030.503 ms/step), a gap of the same order as that of PF (0.5050.505) and RNO (0.5070.507) and therefore within measurement noise rather than an architectural cost. Refine is the only variant with a measurable inference overhead (0.5530.553 ms/step, 9.9%9.9\%), since its iterative correction runs at test time. Peak memory follows the same pattern: HERO uses 67.1867.18 MB, 1.6%1.6\% above the backbone and below Refine (69.8469.84 MB). The cost of HERO is therefore confined to training. Its throughput of 336.92336.92 steps/s is 4.4%4.4\% below PF, the closest rollout baseline, and 10.1%10.1\% below one-step training, reflecting the additional lagged and perturbed rollouts and the diagnostic evaluation; the lagged operator is refreshed only every R steps, so its contribution is amortized. HERO is nonetheless faster to train than Refine (328.47328.47 steps/s), which pays its overhead at both training and inference time. Read together with Table 2, a one-time training slowdown of a few percent buys a 33.8%33.8\% reduction in nRMSE@100, and none of it is charged to deployment. Table 9. Complete results on the four 1D PDE benchmarks (mean ± standard deviation over five seeds). GM100 denotes the geometric-mean error over the 100-step rollout; Stable step denotes the mean number of stable rollout steps. Within each benchmark, the best value per column is shown in bold separately for the FNO and Transolver families. Method Val GM100 ↓ Test GM100 ↓ nRMSE@1 ↓ nRMSE@100 ↓ Stable step ↑ Dispersion (1D) FNO 0.152±0.0050.152± 0.005 0.166±0.0060.166± 0.006 0.0048±0.00010.0048± 0.0001 0.425±0.0170.425± 0.017 69.7±2.669.7± 2.6 FNO-PF 0.144±0.0040.144± 0.004 0.160±0.0050.160± 0.005 0.0045±0.00010.0045± 0.0001 0.402±0.0140.402± 0.014 71.3±2.571.3± 2.5 FNO-Ref. 0.136±0.0040.136± 0.004 0.148±0.0050.148± 0.005 0.0040±0.00010.0040± 0.0001 0.381±0.0140.381± 0.014 73.5±2.873.5± 2.8 FNO-RNO 0.141±0.0040.141± 0.004 0.154±0.0050.154± 0.005 0.0044±0.00010.0044± 0.0001 0.395±0.0150.395± 0.015 72.3±2.472.3± 2.4 FNO-HERO 0.091±0.0020.091± 0.002 0.097±0.0030.097± 0.003 0.0038±0.00010.0038± 0.0001 0.246±0.0080.246± 0.008 80.8±2.580.8± 2.5 Trans. 0.140±0.0040.140± 0.004 0.150±0.0050.150± 0.005 0.0041±0.00010.0041± 0.0001 0.394±0.0160.394± 0.016 72.1±2.672.1± 2.6 Trans.-PF 0.134±0.0040.134± 0.004 0.146±0.0050.146± 0.005 0.0039±0.00010.0039± 0.0001 0.380±0.0150.380± 0.015 72.4±2.772.4± 2.7 Trans.-Refine 0.139±0.0040.139± 0.004 0.152±0.0050.152± 0.005 0.0040±0.00010.0040± 0.0001 0.393±0.0150.393± 0.015 71.9±2.771.9± 2.7 Trans.-RNO 0.136±0.0040.136± 0.004 0.149±0.0050.149± 0.005 0.0041±0.00010.0041± 0.0001 0.382±0.0140.382± 0.014 71.7±2.671.7± 2.6 Trans.-HERO 0.087±0.0020.087± 0.002 0.093±0.0030.093± 0.003 0.0037±0.00010.0037± 0.0001 0.248±0.0090.248± 0.009 81.1±2.281.1± 2.2 Burgers (1D) FNO 0.473±0.0250.473± 0.025 0.526±0.0280.526± 0.028 0.0429±0.00230.0429± 0.0023 0.541±0.0300.541± 0.030 35.6±4.435.6± 4.4 FNO-PF 0.218±0.0090.218± 0.009 0.208±0.0090.208± 0.009 0.0180±0.00080.0180± 0.0008 0.236±0.0120.236± 0.012 68.0±4.468.0± 4.4 FNO-Ref. 0.407±0.0170.407± 0.017 0.457±0.0230.457± 0.023 0.0368±0.00150.0368± 0.0015 0.495±0.0310.495± 0.031 40.4±4.240.4± 4.2 FNO-RNO 0.293±0.0140.293± 0.014 0.311±0.0150.311± 0.015 0.0249±0.00110.0249± 0.0011 0.335±0.0180.335± 0.018 55.0±4.155.0± 4.1 FNO-HERO 0.147±0.0060.147± 0.006 0.138±0.0070.138± 0.007 0.0168±0.00080.0168± 0.0008 0.153±0.0070.153± 0.007 81.7±4.081.7± 4.0 Trans. 0.487±0.0240.487± 0.024 0.537±0.0320.537± 0.032 0.0478±0.00270.0478± 0.0027 0.618±0.0420.618± 0.042 31.2±5.131.2± 5.1 Trans.-PF 0.449±0.0190.449± 0.019 0.505±0.0270.505± 0.027 0.0422±0.00180.0422± 0.0018 0.562±0.0290.562± 0.029 35.5±4.935.5± 4.9 Trans.-Refine 0.428±0.0210.428± 0.021 0.483±0.0270.483± 0.027 0.0401±0.00190.0401± 0.0019 0.535±0.0320.535± 0.032 38.7±4.138.7± 4.1 Trans.-RNO 0.384±0.0160.384± 0.016 0.434±0.0240.434± 0.024 0.0362±0.00160.0362± 0.0016 0.494±0.0290.494± 0.029 41.4±4.241.4± 4.2 Trans.-HERO 0.250±0.0110.250± 0.011 0.281±0.0130.281± 0.013 0.0341±0.00160.0341± 0.0016 0.327±0.0160.327± 0.016 49.3±4.049.3± 4.0 KdV (1D) FNO 0.088±0.0070.088± 0.007 0.074±0.0060.074± 0.006 0.0060±0.00040.0060± 0.0004 3.217±0.2813.217± 0.281 46.4±6.446.4± 6.4 FNO-PF 0.055±0.0040.055± 0.004 0.051±0.0040.051± 0.004 0.0052±0.00030.0052± 0.0003 2.663±0.2022.663± 0.202 48.8±5.648.8± 5.6 FNO-Ref. 0.060±0.0050.060± 0.005 0.056±0.0040.056± 0.004 0.0085±0.00050.0085± 0.0005 2.030±0.1742.030± 0.174 46.1±5.946.1± 5.9 FNO-RNO 0.074±0.0060.074± 0.006 0.069±0.0050.069± 0.005 0.0109±0.00070.0109± 0.0007 1.438±0.1391.438± 0.139 47.5±5.447.5± 5.4 FNO-HERO 0.036±0.0020.036± 0.002 0.034±0.0020.034± 0.002 0.0049±0.00020.0049± 0.0002 0.919±0.0710.919± 0.071 60.6±5.260.6± 5.2 Trans. 0.659±0.0540.659± 0.054 0.816±0.0800.816± 0.080 0.0683±0.00430.0683± 0.0043 1.078±0.1011.078± 0.101 16.4±6.116.4± 6.1 Trans.-PF 0.525±0.0390.525± 0.039 0.681±0.0600.681± 0.060 0.0552±0.00330.0552± 0.0033 0.906±0.0880.906± 0.088 23.4±6.323.4± 6.3 Trans.-Refine 0.562±0.0400.562± 0.040 0.735±0.0540.735± 0.054 0.0599±0.00330.0599± 0.0033 0.977±0.0770.977± 0.077 22.1±5.722.1± 5.7 Trans.-RNO 0.482±0.0370.482± 0.037 0.613±0.0550.613± 0.055 0.0582±0.00400.0582± 0.0040 0.819±0.0780.819± 0.078 23.6±6.223.6± 6.2 Trans.-HERO 0.319±0.0200.319± 0.020 0.401±0.0280.401± 0.028 0.0516±0.00260.0516± 0.0026 0.528±0.0450.528± 0.045 31.8±5.531.8± 5.5 Kuramoto–Sivashinsky (1D) FNO 0.511±0.0290.511± 0.029 0.578±0.0370.578± 0.037 0.0190±0.00110.0190± 0.0011 1.131±0.0771.131± 0.077 43.8±5.743.8± 5.7 FNO-PF 0.401±0.0220.401± 0.022 0.451±0.0260.451± 0.026 0.0179±0.00090.0179± 0.0009 0.972±0.0630.972± 0.063 44.7±4.744.7± 4.7 FNO-Ref. 0.469±0.0240.469± 0.024 0.528±0.0350.528± 0.035 0.0159±0.00090.0159± 0.0009 1.064±0.0731.064± 0.073 46.4±4.646.4± 4.6 FNO-RNO 0.458±0.0270.458± 0.027 0.515±0.0360.515± 0.036 0.0181±0.00090.0181± 0.0009 1.009±0.0701.009± 0.070 43.6±4.643.6± 4.6 FNO-HERO 0.257±0.0120.257± 0.012 0.283±0.0160.283± 0.016 0.0150±0.00070.0150± 0.0007 0.640±0.0370.640± 0.037 53.3±4.353.3± 4.3 Trans. 0.611±0.0360.611± 0.036 0.695±0.0530.695± 0.053 0.0250±0.00140.0250± 0.0014 1.284±0.0971.284± 0.097 37.6±4.937.6± 4.9 Trans.-PF 0.501±0.0290.501± 0.029 0.571±0.0420.571± 0.042 0.0211±0.00120.0211± 0.0012 1.163±0.0861.163± 0.086 41.8±4.641.8± 4.6 Trans.-Refine 0.498±0.0310.498± 0.031 0.572±0.0330.572± 0.033 0.0190±0.00090.0190± 0.0009 1.238±0.0801.238± 0.080 42.4±4.542.4± 4.5 Trans.-RNO 0.477±0.0250.477± 0.025 0.538±0.0320.538± 0.032 0.0201±0.00110.0201± 0.0011 1.099±0.0811.099± 0.081 42.6±5.342.6± 5.3 Trans.-HERO 0.321±0.0180.321± 0.018 0.361±0.0230.361± 0.023 0.0181±0.00090.0181± 0.0009 0.715±0.0400.715± 0.040 49.2±4.549.2± 4.5 Table 10. Complete results on the 2D and 3D PDE benchmarks (mean ± standard deviation over five seeds). GM100 denotes the geometric-mean error over the 100-step rollout; Stable step denotes the mean number of stable rollout steps. Within each benchmark, the best value per column is shown in bold separately for the FNO and Transolver families. Method Val GM100 ↓ Test GM100 ↓ nRMSE@1 ↓ nRMSE@100 ↓ Stable step ↑ Anisotropic Diffusion (2D) FNO 0.713±0.0370.713± 0.037 0.825±0.0570.825± 0.057 0.0875±0.00480.0875± 0.0048 1.868±0.1201.868± 0.120 8.2±3.28.2± 3.2 FNO-PF 0.655±0.0360.655± 0.036 0.754±0.0420.754± 0.042 0.0716±0.00350.0716± 0.0035 1.355±0.0771.355± 0.077 12.9±3.212.9± 3.2 FNO-Ref. 0.567±0.0300.567± 0.030 0.652±0.0400.652± 0.040 0.0611±0.00330.0611± 0.0033 0.786±0.0530.786± 0.053 20.7±3.320.7± 3.3 FNO-RNO 0.692±0.0340.692± 0.034 0.778±0.0490.778± 0.049 0.0846±0.00400.0846± 0.0040 0.925±0.0540.925± 0.054 9.7±3.39.7± 3.3 FNO-HERO 0.363±0.0180.363± 0.018 0.415±0.0190.415± 0.019 0.0571±0.00260.0571± 0.0026 0.521±0.0340.521± 0.034 26.6±2.926.6± 2.9 Trans. 1.168±0.0681.168± 0.068 1.358±0.0801.358± 0.080 0.0799±0.00420.0799± 0.0042 2.968±0.2072.968± 0.207 8.1±3.58.1± 3.5 Trans.-PF 0.965±0.0530.965± 0.053 1.111±0.0681.111± 0.068 0.0696±0.00340.0696± 0.0034 2.172±0.1572.172± 0.157 13.0±3.613.0± 3.6 Trans.-Refine 0.803±0.0390.803± 0.039 0.943±0.0580.943± 0.058 0.0692±0.00370.0692± 0.0037 1.835±0.1191.835± 0.119 14.7±2.914.7± 2.9 Trans.-RNO 0.892±0.0510.892± 0.051 1.048±0.0641.048± 0.064 0.0751±0.00390.0751± 0.0039 2.099±0.1502.099± 0.150 10.4±3.410.4± 3.4 Trans.-HERO 0.537±0.0270.537± 0.027 0.626±0.0340.626± 0.034 0.0652±0.00310.0652± 0.0031 1.185±0.0661.185± 0.066 18.2±2.618.2± 2.6 Kolmogorov Flow (2D) FNO 0.979±0.0720.979± 0.072 1.220±0.1001.220± 0.100 0.4678±0.02690.4678± 0.0269 1.379±0.1061.379± 0.106 1.2±0.81.2± 0.8 FNO-PF 0.731±0.0530.731± 0.053 1.002±0.0661.002± 0.066 0.3402±0.01660.3402± 0.0166 1.255±0.1001.255± 0.100 1.7±0.81.7± 0.8 FNO-Ref. 0.743±0.0460.743± 0.046 0.872±0.0660.872± 0.066 0.3085±0.01660.3085± 0.0166 1.024±0.0901.024± 0.090 2.6±0.72.6± 0.7 FNO-RNO 0.804±0.0520.804± 0.052 1.071±0.0871.071± 0.087 0.3605±0.01840.3605± 0.0184 1.304±0.1001.304± 0.100 1.5±0.61.5± 0.6 FNO-HERO 0.477±0.0250.477± 0.025 0.565±0.0350.565± 0.035 0.2860±0.01400.2860± 0.0140 0.648±0.0480.648± 0.048 3.4±0.63.4± 0.6 Trans. 1.113±0.0741.113± 0.074 1.452±0.1071.452± 0.107 0.3107±0.01860.3107± 0.0186 2.078±0.1812.078± 0.181 1.1±0.71.1± 0.7 Trans.-PF 0.950±0.0630.950± 0.063 1.183±0.0791.183± 0.079 0.2951±0.01890.2951± 0.0189 1.566±0.1141.566± 0.114 1.4±0.71.4± 0.7 Trans.-Refine 0.648±0.0470.648± 0.047 0.790±0.0530.790± 0.053 0.2880±0.01720.2880± 0.0172 0.938±0.0720.938± 0.072 2.9±0.72.9± 0.7 Trans.-RNO 0.755±0.0530.755± 0.053 0.952±0.0770.952± 0.077 0.2987±0.01510.2987± 0.0151 1.131±0.0951.131± 0.095 2.0±0.72.0± 0.7 Trans.-HERO 0.417±0.0270.417± 0.027 0.503±0.0310.503± 0.031 0.2704±0.01260.2704± 0.0126 0.617±0.0500.617± 0.050 3.7±0.63.7± 0.6 Decaying Turbulence / Navier–Stokes (2D) FNO 2.322±0.1572.322± 0.157 2.915±0.2292.915± 0.229 0.0449±0.00280.0449± 0.0028 13.629±1.07713.629± 1.077 17.6±3.817.6± 3.8 FNO-PF 1.441±0.1001.441± 0.100 1.682±0.1281.682± 0.128 0.0349±0.00190.0349± 0.0019 6.271±0.4986.271± 0.498 22.3±3.522.3± 3.5 FNO-Ref. 1.959±0.1331.959± 0.133 2.330±0.1662.330± 0.166 0.0402±0.00240.0402± 0.0024 9.888±0.9109.888± 0.910 18.1±4.018.1± 4.0 FNO-RNO 1.788±0.1071.788± 0.107 2.141±0.1452.141± 0.145 0.0388±0.00200.0388± 0.0020 7.829±0.6477.829± 0.647 20.7±3.420.7± 3.4 FNO-HERO 0.931±0.0590.931± 0.059 1.058±0.0621.058± 0.062 0.0325±0.00150.0325± 0.0015 4.101±0.3234.101± 0.323 26.5±3.426.5± 3.4 Trans. 3.075±0.2463.075± 0.246 3.825±0.3433.825± 0.343 0.0301±0.00200.0301± 0.0020 18.084±1.65218.084± 1.652 20.3±3.820.3± 3.8 Trans.-PF 2.400±0.1782.400± 0.178 2.974±0.2212.974± 0.221 0.0278±0.00160.0278± 0.0016 12.438±1.06712.438± 1.067 24.3±3.524.3± 3.5 Trans.-Refine 2.169±0.1552.169± 0.155 2.591±0.2012.591± 0.201 0.0280±0.00140.0280± 0.0014 10.974±0.86410.974± 0.864 23.1±3.423.1± 3.4 Trans.-RNO 2.224±0.1352.224± 0.135 2.704±0.2132.704± 0.213 0.0291±0.00150.0291± 0.0015 11.438±1.04311.438± 1.043 22.5±4.022.5± 4.0 Trans.-HERO 1.423±0.0911.423± 0.091 1.687±0.1031.687± 0.103 0.0262±0.00120.0262± 0.0012 7.012±0.4977.012± 0.497 25.8±2.925.8± 2.9 Swift–Hohenberg (3D) FNO 0.603±0.0360.603± 0.036 0.678±0.0420.678± 0.042 0.0788±0.00430.0788± 0.0043 0.941±0.0590.941± 0.059 13.4±3.413.4± 3.4 FNO-PF 0.514±0.0260.514± 0.026 0.581±0.0330.581± 0.033 0.0607±0.00270.0607± 0.0027 0.843±0.0530.843± 0.053 20.4±3.620.4± 3.6 FNO-Ref. 0.552±0.0310.552± 0.031 0.619±0.0360.619± 0.036 0.0642±0.00310.0642± 0.0031 0.871±0.0560.871± 0.056 19.6±3.119.6± 3.1 FNO-RNO 0.520±0.0290.520± 0.029 0.596±0.0310.596± 0.031 0.0662±0.00340.0662± 0.0034 0.850±0.0530.850± 0.053 19.5±3.119.5± 3.1 FNO-HERO 0.341±0.0140.341± 0.014 0.381±0.0200.381± 0.020 0.0570±0.00250.0570± 0.0025 0.541±0.0330.541± 0.033 27.7±3.327.7± 3.3 Trans. 0.821±0.0430.821± 0.043 0.912±0.0560.912± 0.056 0.1079±0.00540.1079± 0.0054 1.286±0.0911.286± 0.091 5.3±4.15.3± 4.1 Trans.-PF 0.686±0.0330.686± 0.033 0.782±0.0460.782± 0.046 0.0914±0.00420.0914± 0.0042 1.126±0.0731.126± 0.073 5.9±3.45.9± 3.4 Trans.-Refine 0.633±0.0290.633± 0.029 0.713±0.0400.713± 0.040 0.0823±0.00380.0823± 0.0038 1.004±0.0611.004± 0.061 11.3±3.611.3± 3.6 Trans.-RNO 0.652±0.0320.652± 0.032 0.739±0.0460.739± 0.046 0.0860±0.00380.0860± 0.0038 1.064±0.0671.064± 0.067 7.6±3.37.6± 3.3 Trans.-HERO 0.404±0.0170.404± 0.017 0.450±0.0240.450± 0.024 0.0780±0.00320.0780± 0.0032 0.667±0.0360.667± 0.036 15.0±3.315.0± 3.3 Unbalanced Advection (3D) FNO 2.311±0.1632.311± 0.163 2.647±0.1752.647± 0.175 0.0869±0.00460.0869± 0.0046 7.204±0.6067.204± 0.606 6.4±3.26.4± 3.2 FNO-PF 1.833±0.1231.833± 0.123 2.121±0.1422.121± 0.142 0.0748±0.00400.0748± 0.0040 5.546±0.3945.546± 0.394 10.0±2.510.0± 2.5 FNO-Ref. 2.017±0.1152.017± 0.115 2.317±0.1752.317± 0.175 0.0776±0.00430.0776± 0.0043 6.086±0.4346.086± 0.434 7.8±2.67.8± 2.6 FNO-RNO 1.697±0.1151.697± 0.115 1.931±0.1451.931± 0.145 0.0791±0.00380.0791± 0.0038 4.833±0.3864.833± 0.386 9.0±2.59.0± 2.5 FNO-HERO 1.101±0.0551.101± 0.055 1.244±0.0681.244± 0.068 0.0702±0.00290.0702± 0.0029 3.185±0.1973.185± 0.197 12.7±2.312.7± 2.3 Trans. 2.723±0.1922.723± 0.192 3.106±0.2473.106± 0.247 0.0731±0.00410.0731± 0.0041 7.777±0.6947.777± 0.694 9.8±3.19.8± 3.1 Trans.-PF 2.096±0.1452.096± 0.145 2.426±0.1892.426± 0.189 0.0736±0.00410.0736± 0.0041 6.122±0.4556.122± 0.455 10.3±2.610.3± 2.6 Trans.-Refine 1.901±0.1161.901± 0.116 2.201±0.1532.201± 0.153 0.0713±0.00390.0713± 0.0039 5.319±0.4085.319± 0.408 9.6±2.89.6± 2.8 Trans.-RNO 1.784±0.1221.784± 0.122 2.079±0.1562.079± 0.156 0.0742±0.00410.0742± 0.0041 5.015±0.3875.015± 0.387 10.6±2.910.6± 2.9 Trans.-HERO 1.197±0.0641.197± 0.064 1.389±0.0881.389± 0.088 0.0673±0.00320.0673± 0.0032 3.175±0.2493.175± 0.249 12.1±2.312.1± 2.3