Paper deep dive

Eigenanalysis framework for autoregressive neural emulators of multi-scale chaotic dynamics

Conrad Ainslie, Pedram Hassanzadeh, Michael W. Mahoney, Ashesh Chattopadhyay

Year: 2026Venue: arXiv preprintArea: cs.AIType: PreprintEmbeddings: 127

Intelligence

Status: succeeded | Model: Gemma-4-26B-A4B | Prompt: intel-v1 | Confidence: 96%

Last extracted: 8/23/2026, 1:22:56 AM

Summary

This paper introduces an eigenanalysis framework to understand the long-term instability and error growth in neural autoregressive emulators of chaotic systems. By analyzing the Jacobian of the learned one-step update map with respect to the state, the authors demonstrate that model stability is governed by the spectral radius. Direct-step models typically exhibit unstable eigenvalues (magnitude > 1), leading to rapid divergence, while integration-constrained models collapse their eigenspectrum onto the unit circle, ensuring neutral stability. The study proposes a stability-promoting loss function to regularize Jacobian-driven error amplification, validated across 29 models on the Kuramoto-Sivashinsky system.

Entities (8)

Jacobian · mathematical-object · 99%Kuramoto-Sivashinsky system · physical-system · 99%Eigenanalysis framework · method · 98%Direct-step models · model-architecture · 97%Integration-constrained models · model-architecture · 97%Stability-promoting loss · loss-function · 96%Neural autoregressive models · model-class · 95%Spectral Radius · mathematical-property · 95%

Relation Signals (6)

Eigenanalysis framework → validatedon → Kuramoto-Sivashinsky system

confidence 99% · Demonstrated across 29 models... on the Kuramoto-Sivashinsky system

Eigenanalysis framework → analyzes → Jacobian

confidence 98% · By analyzing the Jacobian of the learned one-step update map with respect to the state, we show how inference-time error growth... is governed by its spectral radius.

Direct-step models → exhibits → unstable eigenvalues

confidence 97% · Direct-step architectures... generically admit unstable eigenvalues with magnitudes exceeding one, explaining the rapid divergence

Integration-constrained models → exhibits → neutral stability

confidence 97% · integration-constrained models... collapse their eigenspectrum onto the unit circle, yielding neutral stability

Stability-promoting loss → regularizes → Jacobian-driven error amplification

confidence 96% · we introduce a stability-promoting loss that explicitly regularizes Jacobian-driven error amplification

Jacobian → governs → model stability

confidence 95% · inference-time error growth, and thus model stability, is governed by its spectral radius.

Cypher Suggestions (0)

No Cypher suggestions yet.

Abstract

Abstract:Neural autoregressive models have rapidly emerged as powerful emulators of high-dimensional chaotic systems, yet their long-term instability and error growth remain poorly understood, leading to ad-hoc solutions. Here, we develop an eigenanalysis framework that reveals the dynamical origin of this error growth. By analyzing the Jacobian of the learned one-step update map with respect to the state, we show how inference-time error growth, and thus model stability, is governed by its spectral radius. Direct-step architectures (models that predict the next state from the previous one) generically admit unstable eigenvalues with magnitudes exceeding one, explaining the rapid divergence of these widely used models. In contrast, integration-constrained models (where the time derivative is estimated and integrated with a higher-order integrator) collapse their eigenspectrum onto the unit circle, yielding neutral stability and a universal linear error-scaling law. The largest eigenvalue of this Jacobian provides an architecture-agnostic, a priori diagnostic of short-term skill, long-term stability, and spectral bias, without requiring an expensive rollout. Leveraging this theory, we introduce a stability-promoting loss that explicitly regularizes Jacobian-driven error amplification, improving both forecast accuracy and dynamical robustness. Demonstrated across $29$ models spanning two architectures, several explicit and implicit integrators, and multiple loss functions on the Kuramoto-Sivashinsky system, our results establish a theoretical foundation for the design and evaluation of neural emulators of chaotic multi-scale dynamics. More broadly, our framework is a step toward the kind of a priori stability analysis that numerical analysis provides for discretizations of differential equations and that scientific machine learning currently lacks.

PDF

Open source PDF →Open local PDF →

Trouble viewing inline? Open PDF directly →

Full Text

123,473 characters extracted from source content.

Expand or collapse full text

Eigenanalysis framework for autoregressive neural emulators of multi-scale chaotic dynamics Conrad Ainslie1, Pedram Hassanzadeh2,3, Michael W. Mahoney4,5,6 and Ashesh Chattopadhyay 1 Thanks: aschatto@ucsc.edu Affiliation: 1Department of Applied Mathematics, University of California, Santa Cruz, Santa Cruz, 95064, CA Affiliation: 2Department of Geophysical Sciences, University of Chicago, Chicago, IL 60637 Affiliation: 3Committee on Computational and Applied Mathematics, University of Chicago, Chicago, IL 60637 Affiliation: 4International Computer Science Institute, Berkeley, CA 94720 Affiliation: 5Lawrence Berkeley National Laboratory, Berkeley, CA 94720 Affiliation: 6Department of Statistics, University of California, Berkeley, CA 94720 Abstract Neural autoregressive models have rapidly emerged as powerful emulators of high-dimensional chaotic systems, yet their long-term instability and error growth remain poorly understood, leading to ad-hoc solutions. Here, we develop an eigenanalysis framework that reveals the dynamical origin of this error growth. By analyzing the Jacobian of the learned one-step update map with respect to the state, we show how inference-time error growth, and thus model stability, is governed by its spectral radius. Direct-step architectures (models that predict the next state from the previous one) generically admit unstable eigenvalues with magnitudes exceeding one, explaining the rapid divergence of these widely used models. In contrast, integration-constrained models (where the time-derivative is estimated and integrated with a higher-order integrator) collapse their eigenspectrum onto the unit circle, yielding neutral stability and a universal linear error-scaling law. The largest eigenvalue of this Jacobian provides an architecture-agnostic, a priori diagnostic of short-term skill, long-term stability, and spectral bias, without requiring an expensive rollout. Leveraging this theory, we introduce a stability-promoting loss that explicitly regularizes Jacobian-driven error amplification, improving both forecast accuracy and dynamical robustness. Demonstrated across 29 models spanning two architectures, several explicit and implicit integrators, and multiple loss functions on the Kuramoto–Sivashinsky system, our results establish a theoretical foundation for the design and evaluation of neural emulators of chaotic multi-scale dynamics. More broadly, our framework is a step toward the kind of a priori stability analysis that numerical analysis provides for discretizations of differential equations, and that scientific machine learning currently lacks. Significance Statement Neural networks that step a physical system forward in time now rival traditional simulations of chaotic, multi-scale flows such as the atmosphere, climate, and ocean, and they can run far faster than traditional numerical methods. They can also drift or blow up over long rollouts. Remedies are currently found by trial and error, each tested by running the simulation that the network was built to replace. Here, we develop an eigenanalysis theory of these networks, which shows that one quantity, obtained from a single derivative of the trained model, predicts how fast its error grows, before any simulation is run. We find that networks that predict the next state directly amplify errors far faster than the physics allows, while networks that instead integrate a learned rate of change do not. 1 Introduction Neural autoregressive models have become effective data-driven emulators of high-dimensional chaotic dynamical systems [Fan et al. 2020, Li et al. 2022, Floryan 2024]. Their most visible success is in Earth system modeling, where they now rival or outperform operational numerical weather prediction [Pathak et al. 2022, Lam et al. 2022, Bi et al. 2023, Price et al. 2025]. Despite this progress, long-term emulation remains difficult. Many models become unstable or excessively diffusive over long rollouts [Lippe et al. 2023, Chattopadhyay et al. 2023b, Keisler 2022, Lai et al. 2025, Pedersen et al. 2025]. A growing number of models are now stable over long horizons, but this stability is typically achieved through extensive hyperparameter optimization; and even carefully tuned models can still blow up. Each such failure tends to be addressed with a new, problem-specific fix rather than by identifying the mechanism that produced it, leaving the field without a principled account of when and why these models lose stability [Watt-Meyer et al. 2025, Guan et al. 2025, Pedersen et al. 2025, Sambamurthy and Chattopadhyay 2025]. A leading explanation for this instability in multi-scale systems, such as turbulent flows with a decaying energy spectrum, is the spectral bias of the network [Chattopadhyay et al. 2023b, Yu et al. 2024], whereby a model trained to predict a single step fails to capture the high-wavenumber dynamics. This error is thought to grow during autoregressive rollout through the coupling between small and large scales until it corrupts the resolved dynamics. If spectral bias is the underlying cause, however, its consequences are realized only through the way one-step errors accumulate from step to step. The role of this error propagation, as distinct from the one-step error itself, has not been made precise, and most stabilization strategies remain disconnected from any quantitative theory of how errors grow during rollout inference. One line of work points to the integration scheme used to advance the model in time. Krishnapriyan et al. 2023 showed that constraining a model with a higher-order integrator, such as fourth-order Runge–Kutta, tends yield convergent integration analogous to classical numerical schemes. Chattopadhyay et al. 2023b combined such a hard constraint with a spectral regularizer and reported reduced error growth and longer stable rollouts, and subsequent work has scaled this idea to coupled climate and high-resolution ocean models, producing stable rollouts over centuries of integration [Chattopadhyay et al. 2024, Guan et al. 2025, Lupin-Jimenez et al. 2025]. Related stabilization strategies include refinement [Lippe et al. 2023], noise-based regularization [Stachenfeld et al. 2021, Wikner et al. 2022], and stochastic modeling [Chattopadhyay et al. 2023a, Pedersen et al. 2025, Sambamurthy and Chattopadhyay 2025]. Despite this empirical evidence, it remains unclear why an integration constraint should help. By analogy with numerical analysis, it is natural to assume that incorporating physical insight into a constraint inspired by time integration may aid stability; however, intuitions about combining physical reasoning and machine learning can be misleading [Krishnapriyan et al. 2021, Krishnapriyan et al. 2023, Sakarvadia et al. 2025], and learned spatiotemporal models can behave counterintuitively [Yu et al. 2025], including how performance scales with choice of the discrete time step [Bi et al. 2023, Chattopadhyay et al. 2022]. A general theory of inference-time stability for these models is still lacking, with rigorous results largely confined to the linear case [Floryan 2024], and the effect of the integration scheme, loss function, and architecture on accuracy and stability has not been studied systematically. This raises a set of questions. Does an integration constraint control the rate at which errors grow during rollout, and does it also improve short-term accuracy? Can a single quantity, computed before any long rollout, predict both? Useful answers should not depend on the particular architecture, integration scheme, loss function, or system, and they should not require long emulations from many initial conditions. We answer these questions through an eigenanalysis of the Jacobian of the learned one-step update map (Fig. 1). This Jacobian is taken with respect to the system state, which is the derivative that sets the stability of a discrete map in dynamical systems theory [Floryan 2024], and not with respect to the network parameters, which is the derivative returned by backpropagation and the one whose spectrum is usually studied in machine learning [Liao and Mahoney 2021]. The two derivatives have different arguments and different roles, and only the state derivative enters the error dynamics at inference. Linearizing the error dynamics shows that inference-time error growth is governed by the spectral radius of this Jacobian, so its largest eigenvalue is an a priori diagnostic that predicts error growth during inference. This can be accomplished without resorting to expensive long rollouts, and it is independent of architecture, integration scheme, and loss function. This single quantity diagnoses the amplification of error during rollout; together with the one-step error, it accounts for short-term accuracy and the growth of spectral bias (Section 3.3). The analysis also explains the empirical record. Direct-step models, which predict the next state directly, generically admit eigenvalues with magnitude well above one and have no mechanism to control them, which accounts for their rapid divergence. Integration-constrained models, in which the network estimates a time derivative that is advanced by a numerical integrator, have a Jacobian whose eigenvalues collapse onto the unit circle and yield neutral stability together with a linear error-scaling law. Building on this, we introduce a stability-promoting loss that directly regularizes the Jacobian-driven amplification of errors and improves both accuracy and stability. We develop this theory for the Kuramoto–Sivashinsky (KS) system, and we evaluate it across the full suite of models studied here, spanning two architectures, several explicit and implicit integration schemes, and loss functions with and without spectral regularization. The remainder of the paper presents the linear stability theory, the empirical results, the limits of the eigenvalue-only view and its distinctions from classical linear stability analysis, and the methods, including all the explicit integration schemes and an implicit integration scheme based on implicit-layer and deep-equilibrium models [Kawaguchi 2021]. The full model suite is listed in SI, Section 1. The results obtained with the implicit constraint, together with the derivation of its Jacobian, are reported in SI, Section 2. The Fourier-space analysis of spectral bias is given in SI, Section 3, and the normality analysis that supports the eigenvalue approximation is given in SI, Section 4. Figure 1: Schematic representation of the main contributions of this paper. An eigenvalue decomposition of the Jacobian of the learned one-step update map, taken with respect to the system state and not with respect to the network parameters, yields a linear theory for the role the Jacobian plays in short-term error growth, the rate of error accumulation over rollout, and the spectral fidelity of the emulation. Direct-step models, which predict the next state in one shot, admit eigenvalues outside the unit circle and diverge during rollout, whereas integration-constrained models, in which the network estimates a time derivative that is advanced by a numerical integrator, collapse their eigenspectrum onto the unit circle and remain neutrally stable, so that a perturbation is neither amplified nor damped from one step to the next (Section 2.3). The resulting stability-promoting loss regularizes the Jacobian-driven amplification of error, reducing the rollout error-growth rate: for some models this appears directly as eigenvalues pulled toward the unit circle (illustrated), while more generally it suppresses the realized, direction-dependent amplification of the error without necessarily changing |λmax|| _ | (Section 3.4). 2 Linear stability analysis We analyze the inference-time stability of a trained neural autoregressive model by linearizing its error dynamics during emulation, following the linear stability analysis of discrete dynamical systems. 2.1 Learned dynamics and the Jacobian Let the true system advance one step through an operator G, T(t+Δt)=⁡[T(t)],u_T(t+ t)=G [u_T(t) ], (1) with Tu_T the true state, usually available only at the initial condition t=t0t=t_0. A neural emulator approximates G by a learned operator ~ G and advances its own predicted state pu_p, p(t+Δt)=~[p(t)].u_p(t+ t)= G [u_p(t) ]. (2) Both states are functions of time. From here on we suppress the time argument whenever it is not needed and write Tu_T and pu_p, restoring it only where the rollout step matters. The learned operator is built from a neural network ⁡(⋅,θ)N(\,·\,;θ) with parameters θ (see Section 6.2). A direct-step model sets ~[T]=⁡(T,θ) G[u_T]=N(u_T;θ), predicting the next state in one shot. An integration-constrained model instead estimates a time derivative with the network and advances it with a numerical integrator H, ~[T]=T+Δt[⁡(T,θ)]. G[u_T]=u_T+ t\,H\! [N(u_T;θ) ]. (3) Equation (3) is the explicit form of the constraint. The implicit form, in which the network is evaluated at the updated state, is given in Section 6.4. The central object of our analysis is the Jacobian of the learned update map with respect to the state, =∇~[T].J=∇ G[u_T]. (4) Throughout, ∇ denotes the gradient with respect to the state Tu_T, evaluated at the current true state, and not the gradient with respect to the network parameters θ. In machine learning, the derivative of a trained model most often refers to the derivative with respect to θ, which is what backpropagation returns and what the optimizer acts on. Equation (4) is the other derivative, taken at fixed θ, and it is the one that appears in the error recursion below. The Jacobian is the sensitivity of the one-step prediction to a perturbation of its input, which is the quantity that governs how errors propagate during rollout. It is distinct from the parameter gradients used in training, and it is rarely examined when these models are evaluated. 2.2 Error dynamics Define the prediction error ⁡(t)=T(t)−p(t)e(t)=u_T(t)-u_p(t). Subtracting Eq. (2) from Eq. (1) gives ⁡(t+Δt)=⁡[T]−~[p].e(t+ t)=G [u_T ]- G [u_p ]. (5) Linearizing ~ G about the true state yields ⁡(t+Δt)=⁡[T]−~[T]⏟ϵ⁡(t)+∇~[T]⏟(t)+⁡(‖⁡(t)‖22).e(t+ t)= G [u_T ]- G [u_T ]_ ε(t)+ ∇ G [u_T ]_Je(t)+O\! (\|e(t)\|_2^2 ). (6) The first term ϵ⁡(t) ε(t) is what we call the generalization error, the one-step error of the learned operator at the true state, and the quantity that autoregressive models are trained to minimize. The second term, ⁡(t)Je(t), is what we call the propagated error, the accumulated error carried forward through the linearized map, and J in this role is the error propagator. Training and validation control only ϵ⁡(t) ε(t), whereas J, which sets the size of the propagated error, is rarely examined. At the initial condition the two error notions coincide. Since p(t0)=T(t0)u_p(t_0)=u_T(t_0), the first step incurs only the generalization error, so ⁡(t0+Δt)=ϵ⁡(t0)e(t_0+ t)= ε(t_0). Writing tn=t0+nΔt_n=t_0+n t and unrolling Eq. (6) over n steps, we obtain ⁡(tn) (t_n) =ϵ⁡(tn−1)+∑i=0n−2(∏j=i+1n−1⁡(tj))ϵ(ti)+⁡(‖⁡(tn−1)‖22), = ε(t_n-1)+ _i=0^n-2 ( _j=i+1^n-1J(t_j) ) ε(t_i)+O\! (\|e(t_n-1)\|_2^2 ), (7) with ⁡(tj)=∇~[T(tj)]J(t_j)=∇ G[u_T(t_j)] the Jacobian at step j. Total error growth therefore has two sources. The generalization error ϵ ε reflects intrinsic model quality and can be reduced through better data, architecture, or capacity. The product of Jacobians is the error propagator over the intervening steps, and it determines whether those errors are amplified or damped, thus controlling stability during model rollout. Equation (7) cleanly separates model fidelity, carried by ϵ ε, from dynamical stability, carried by J. Characterizing J across architectures, loss functions, and integration constraints, and relating it to model quality and stability during emulation, is the main contribution of this work. 2.3 Eigenvalues and the Jacobian of each model class The Jacobian admits the eigendecomposition =−1,J=V ^-1, (8) with V the eigenvector matrix and the diagonal matrix of eigenvalues λ1,…,λd _1,…, _d for a system of dimension d. The magnitude of the largest eigenvalue |λmax|| _ | sets the local stability. The dynamics are locally unstable when |λmax|>1| _ |>1, in which case a perturbation of the state is amplified at every step and grows geometrically in the step count, and neutrally stable when |λmax|=1| _ |=1, in which case a perturbation is neither amplified nor damped from one step to the next, so an error already present is carried forward at fixed size. For non-normal J, as is typical of neural-network Jacobians, the eigenvalues bound the asymptotic growth of perturbations but not their transient amplification, which also depends on the orientation of the error relative to the eigenvectors; Section 3.3 quantifies this effect. SI, Section 4 shows that the integration constraint renders J near-normal, with a departure from normality of ⁡(Δt2)O( t^2), so the eigenvalue approximation used below is valid. The two model classes have very different Jacobians. For a direct-step model, =∇(T,θ),J= (u_T;θ), (9) which places no intrinsic constraint on the spectrum. On the other hand, differentiating the integration-constrained map of Eq. (3) yields =+Δt∇[(T;θ)].J=I+ t\, \! [N(u_T;θ) ]. (10) Because Δt≪1 t 1, the eigenvalues of Eq. (10) of an integration-constrained model stay close to one, and the model is near-neutrally stable by design. Direct-step models lack this control, since their Jacobian has no explicit dependence on Δt t, and are thus prone to instability. Equation (10) is the Jacobian of the explicitly constrained map. Differentiating the implicit map instead gives =(−Δt∇[])−1J=(I- t\, [N])^-1, which reduces to Eq. (10) to first order in Δt t. SI, Section 2.2 derives this expression and shows that the inverse form places an eigenvalue outside the unit circle only when the learned tangent operator has an expanding direction. The constraint fixes the form of J but not the size of the learned tangent operator ∇[(T,θ)] [N(u_T;θ)], and the boundedness of that operator is a property the network acquires during training, rather than being an algebraic consequence of Eq. (10). For the true KS tangent operator at the resolution used here, Δtkmax4=1.07×103 t\,k_ ^4=1.07× 10^3, so an explicit Euler step applied to the exact linearization would be unstable by three orders of magnitude. The learned operators instead satisfy ∥Δt∇[]∥∼10−3\| t\, [N]\| 10^-3, read from the spectral spread of Fig. 2(d). The near-unit spectrum therefore reflects the smoothness of the learned time derivative at high wavenumbers, which is the same property that produces the spectral bias analyzed in SI, Section 3, and not the algebra of +Δt∇I+ t\, alone. 2.4 Error scaling for integration-constrained models We first derive the error growth over a single step, and we then extend it to p steps for as long as the linear approximation holds. From the linearized recursion of Eq. (6), the error advances by one step as ⁡(tk+1)=ϵ⁡(tk)+⁡(tk)(tk)+⁡(‖⁡(tk)‖22)e(t_k+1)= ε(t_k)+J(t_k)e(t_k)+O(\|e(t_k)\|_2^2), with ⁡(tk)J(t_k) the Jacobian at step k. For an integration-constrained model (tk)=+Δt∇[(T(tk))]J(t_k)=I+ t\, [N(u_T(t_k))] from Eq. (10), so (tk+1)=ϵ(tk)+(+Δt∇[(T(tk))])(tk)+(∥(tk)∥22).e(t_k+1)= ε(t_k)+ (I+ t\, \! [N(u_T(t_k)) ] )e(t_k)+O\! (\|e(t_k)\|_2^2 ). (11) At the initial condition p(t0)=T(t0)u_p(t_0)=u_T(t_0), so ⁡(t0)=e(t_0)=0 and the first step incurs only the generalization error, ⁡(t1)=ϵ⁡(t0)e(t_1)= ε(t_0). Applying Eq. (11) once more, at k=1k=1, (t2)=ϵ(t1)+(+Δt∇[(T(t1))])(t1)+(∥(t1)∥22).e(t_2)= ε(t_1)+ (I+ t\, \! [N(u_T(t_1)) ] )e(t_1)+O\! (\|e(t_1)\|_2^2 ). (12) Because Δt t is small and ~ G is continuous, the generalization error is nearly unchanged over one step, ϵ⁡(t1)≈ϵ⁡(t0)=⁡(t1) ε(t_1)≈ ε(t_0)=e(t_1). Substituting this into the first term gives (t2)=(2+Δt∇[(T(t1))])(t1)+(∥(t1)∥22).e(t_2)= (2I+ t\, \! [N(u_T(t_1)) ] )e(t_1)+O\! (\|e(t_1)\|_2^2 ). (13) The matrix acting on ⁡(t1)e(t_1) in Eq. (13) is 22I plus a correction of order Δt t, so for Δt≪1 t 1 its largest eigenvalue is close to 22 and ‖⁡(t2)‖≈2|λmax|‖⁡(t1)‖\|e(t_2)\|≈ 2\,| _ |\,\|e(t_1)\|. This is due to the tight clustering of the spectrum. Since the eigenvalues of J are nearly coincident (Fig. 2(c,d)), ‖≈|λmax()|‖\|Me\|≈| _ (M)|\,\|e\|, where M is any square matrix with tightly clustered eigenvalues at unity, holds regardless of the orientation of e, and the spectrum of (2+Δt∇[(T(t1))]) (2I+ t\, \! [N(u_T(t_1)) ] ) coincides with 2|λmax|2\,| _ | up to ⁡(Δt)O( t). The error therefore doubles at the first comparison, which is the k=1k=1 reference line in Fig. 3(a). SI, Section 4 shows that this directional insensitivity follows from the near-identity structure of Eq. (10) and not from the clustering of the eigenvalues alone. Section 3.3 examines the cases where this directional insensitivity breaks down. We now extend this to p steps. Iterating Eq. (11) from the initial condition and substituting the integration-constrained Jacobian gives ⁡(tp) (t_p) =∑i=0p−1ϵ⁡(ti)+Δt∑i=0p−2∑j=i+1p−1∇[⁡(T(tj))]ϵ(ti)+⁡(‖⁡(tp−1)‖22)+⁡(Δt2), = _i=0^p-1 ε(t_i)+ t _i=0^p-2 _j=i+1^p-1 \! [N(u_T(t_j)) ] ε(t_i)+O\! (\|e(t_p-1)\|_2^2 )+O\! ( t^2 ), (14) whose leading term is the running sum of generalization errors. The same two approximations apply for as long as the rollout error stays small. The ⁡(Δt2)O( t^2) remainder and the explicit Δt t double sum are subdominant because pΔt≪1p\, t 1 over the step counts considered here, and the continuity of ~ G keeps the generalization error nearly constant across the early rollout, ϵ⁡(ti)≈ϵ⁡(t0)=⁡(t1) ε(t_i)≈ ε(t_0)=e(t_1) for i=0,…,p−1i=0,…,p-1. The leading sum then collapses to p(t1)p\,e(t_1), and (tp)=(p+Δt∑j=1p−1j∇[(T(tj))])(t1)+(∥(tp−1)∥22)+(Δt2).e(t_p)= (p\,I+ t _j=1^p-1j\, \! [N(u_T(t_j)) ] )e(t_1)+O\! (\|e(t_p-1)\|_2^2 )+O\! ( t^2 ). (15) The bracket is now pp\,I plus a correction that is ⁡(pΔt)O(p\, t) relative to the leading term, so for pΔt≪1p\, t 1 its largest eigenvalue is close to p, giving the linear scaling law ‖⁡(tp)‖≈p|λmax|‖⁡(t1)‖,\|e(t_p)\|≈ p\,| _ |\,\|e(t_1)\|, (16) with |λmax|≈1| _ |≈ 1 the largest eigenvalue of the one-step Jacobian J. The 22I result of Eq. (13) is the p=2p=2 case. The slope grows linearly with the step count p, reproducing the slopes 22, 1111, 101101, and 10011001 measured at lead times k=1,10,100,1000k=1,10,100,1000 (that is, p=2,11,101,1001p=2,11,101,1001) in Fig. 3. Because |λmax|≈1| _ |≈ 1, models separate along these lines only through the small departures of |λmax|| _ | from one, so a larger |λmax|| _ | corresponds to faster error growth. The law holds only while both approximations remain valid, that is, while Δt≪1 t 1 and while ‖⁡(tp)‖\|e(t_p)\| stays small enough that the ⁡(‖22)O(\|e\|_2^2) terms are negligible. Explicitly, the derivation requires pΔt≪1p\, t 1, i.e., p≪1/Δt=1000p 1/ t=1000 at Δt=10−3 t=10^-3; the longer rollouts of Figs. 2(b) and 3(d) probe beyond this window, and the continued slow accumulation observed there is an empirical observation outside the guaranteed regime rather than a prediction of Eq. (16). The two models trained at Δt=5×10−2 t=5× 10^-2 and Δt=10−1 t=10^-1 exceed this window by a larger margin at the longest lead times, as discussed in SI, Section 2.4. As the rollout proceeds and the accumulated error grows, these higher-order terms become significant and the error departs from the linear prediction. Within the lead times of Fig. 3 this departure is visible first for the direct-step models, whose large per-step amplification excites the higher-order terms after only a few steps. 3 Results We evaluate the Jacobian eigenanalysis of Section 2 as an a priori diagnostic on the Kuramoto–Sivashinsky (KS) system. The suite studied here comprises 29 models, of which two are direct-step models and 27 are integration-constrained (see Methods). The main-text figures show a representative subset of these models, and the complete set of 29 is reported in Table S1. 3.1 Eigenvalues a priori diagnose the amplification of rollout error Figure 2 contrasts the two model classes using four models, the direct-step MLP and FNO and their Euler-constrained counterparts. The two Euler models stand in for all 27 integration-constrained models examined here (SI, Fig. S1 and SI, Fig. S2), and each model is trained and emulated at Δt=ΔtDNS=10−3 t= t_DNS=10^-3, unless stated otherwise. Throughout, a model whose Jacobian has the form of Eq. (10) is called integration-constrained, and every model except the two direct-step models is of this type. For the three implicit models, this holds to first order in Δt t (SI, Section 2.2). The central empirical result is in Fig. 2(a). The two direct-step models become unstable and unphysical within a few steps, with rapidly growing error, whereas the integration-constrained models remain stable and physically plausible, with error growing in proportion to the step count and orders of magnitude smaller. Figure 2(b) shows the same RMSE on log–log axes over the first 10410^4 steps of the emulation. The direct-step models saturate within a few tens to a few hundreds of steps, whereas the error of the two Euler models accumulates slowly and in proportion to the number of rollout steps, as predicted by Eq. (16), remaining orders of magnitude smaller over most of the window, with the Euler MLP approaching the decorrelation level only at its right edge. Strict proportionality to the step count does eventually give way as the accumulated error grows and the higher-order terms of Section 2.4 become significant, but the error continues to accumulate slowly rather than diverging, in contrast to the direct-step models. Equation (7) separates two budgets: the eigenspectrum sets the rate at which errors are amplified, while the one-step generalization error sets the rate at which they are injected. The integration constraint alone, however, does not guarantee long-term stability. The near-neutral spectrum controls how fast errors are amplified, but the total error also carries the accumulated one-step generalization error ϵ ε (Eq. (14)). The Euler MLP of Fig. 2(b), whose ϵ ε is the larger of the two, reaches the climatological decorrelation level by the right edge of that window, while the Euler FNO does not, and the same ordering holds for the remaining MLP and FNO pairs over the full 10510^5-step emulation (SI, Fig. S1), despite eigenvalues on the unit circle in every case. Long-term stability therefore requires both a near-neutral spectrum and a sufficiently small generalization error, with or without the spectral regularizer (Section 6.6), and the fidelity of the predicted spectrum still depends on architecture (SI, Section 3). This behavior is set by the eigenvalues of the Jacobian, evaluated at inference rather than during training, and computed once at the initial condition t0t_0. Figure 2(c) plots these eigenvalues on the Argand plane for all four models. The two direct-step models have |λmax|>1| _ |>1, while the eigenvalues of the two Euler models collapse onto the unit circle and cluster tightly around 11. The 2×10242× 1024 eigenvalues of the two integration-constrained models are nearly coincident there, so they appear as a single dense cluster and are barely visible in (c). Figure 2(d) therefore zooms into that cluster and shows a spread confined to the third decimal place, with the two models separated in the fourth, consistent with neutral stability. The measured values can be compared against the dynamics being emulated. The linear part of the KS operator has symbol k2−k4k^2-k^4, whose maximum over the resolved wavenumbers is 0.24950.2495 at k=0.69k=0.69, so the corresponding time-Δt t flow map has |λmax|=exp⁡(0.2495Δt)=1.00025| _ |= (0.2495\, t)=1.00025. The departures from unity measured for the two Euler models, 1.75×10−41.75× 10^-4 and 2.01×10−42.01× 10^-4, agree with the value implied by the linear operator, 2.5×10−42.5× 10^-4, to within a factor of 1.51.5. This estimate neglects the advective contribution to the tangent map and is a reference scale rather than a precise target. It matters because Eq. (10) would place the spectrum within ⁡(Δt)O( t) of unity for any bounded learned operator, including one that has learned the wrong dynamics, so proximity to unity is not by itself evidence that the model has learned anything. Agreement to within a factor of 1.51.5 with the value implied by the linear operator indicates instead that the trained network has recovered the size of the true tangent map and has not merely inherited the near-identity structure of the constraint. It also gives |λmax|| _ | a reference value against which a trained model can be checked before any rollout is run, since a model whose measured |λmax|| _ | departs from the value implied by the linear operator is misrepresenting the local growth rate of the system it emulates. This eigenvalue structure explains the error growth in (a) and (b). The near-unit eigenvalues of the integration-constrained models produce the slow, near-neutral accumulation derived in Section 2.4, while the large |λmax|| _ | of the direct-step models (Eq. (9)) amplifies perturbations strongly and drives the rapid divergence. The collapse onto the unit circle follows directly from Eq. (10). The identity term contributes eigenvalues exactly at 11, and the correction Δt∇[(T)] t\, [N(u_T)] is a small perturbation because Δt=10−3 t=10^-3. The qualitative spectrum, and hence the predicted neutral stability, is therefore insensitive to the architecture, the loss function, the presence or absence of spectral regularization, and the order of the explicit integrator. This insensitivity is why two Euler models can stand in for all 27 integration-constrained models in Fig. 2. SI, Fig. S1 shows the RMSE curves and the eigenspectra of the remaining explicit models, which exhibit the same behavior. 3.2 Emergence of a linear scaling law Figure 3 tests the scaling law of Section 2.4 for all 29 models, plotting ‖⁡(t+Δt)‖\|e(t+ t)\| against |λmax|‖⁡(t)‖| _ |\,\|e(t)\| at successive rollout steps, with |λmax|| _ | the largest eigenvalue of the Jacobian J evaluated at the initial condition. The two direct-step models are shown in vermilion and the 27 integration-constrained models—spanning the Euler, RK4, PEC4, and implicit Euler schemes—in blue, with architecture distinguished by marker fill (MLP filled, FNO open). Table S1 lists each model with its integration scheme, architecture, and loss. The integration-constrained models follow the predicted law ‖⁡(tp)‖≈p|λmax|‖⁡(t1)‖\|e(t_p)\|≈ p\,| _ |\,\|e(t_1)\| of Eq. (16). At rollout step p they lie on a line whose slope is the step index p itself, equal to 2 at the first step and growing to 11, 101, and 1001 at p=11,101,1001p=11,101,1001; on the log–log axes of the four panels of Fig. 3 these appear as unit-slope reference lines, one per panel, whose vertical offset grows with p. Because |λmax|≈1| _ |≈ 1 for these models, they separate along each line only through the small departures of |λmax|| _ | from one. The four panels of Fig. 3 correspond to pΔt=0.002p\, t=0.002, 0.0110.011, 0.1010.101, and 1.0011.001, so they traverse the validity window of Eq. (16) from well inside it to its boundary. The law holds tightly in (a) and (b), where pΔt≪1p\, t 1. The scatter about the reference line grows in (c), where pΔt≈0.1p\, t≈ 0.1. In (d), where pΔt≈1p\, t≈ 1, the condition is no longer met and the models depart from the line in both directions. The breakdown therefore occurs at the step count the derivation predicts, and the agreement in panels (a) to (c) is a test of the law within the regime where it is derived rather than an extrapolation beyond it. The two direct-step models track the same line for the first few steps and then peel away, growing faster than linearly. This is consistent with the law being a small-error prediction. The eigenvalues of the direct-step models are not clustered near one, so once the accumulated error grows enough to excite the dominant eigenvalue with |λmax|>1| _ |>1, the higher-order terms take over. The departure occurs after about 10 steps for the direct-step MLP and about 100 steps for the direct-step FNO, which places them above the reference line at the lead time of panel (b). Beyond that point, their error reaches the climatological decorrelation level and stops growing, so at the longer lead times of panels (c) and (d) they fall below a reference line that continues to grow in proportion to p. This saturation is the same one visible in Fig. 2(a). The slope p has a simple origin (Eq. (15)). At the first step ⁡(t1)e(t_1) is the generalization error, and the Jacobians of integration-constrained models have |λmax|≈1| _ |≈ 1 for Δt≪1 t 1, so the generalization error is nearly unchanged over the early rollout and Eq. (14) reduces to a sum of p near-identical copies of ⁡(t1)e(t_1). As p grows, the accumulated error eventually becomes large enough that the ⁡(‖22)O(\|e\|_2^2) terms are no longer negligible, and the points scatter about the line, which is what panel (d) shows. Enforcing an integration constraint is what makes this linear stability analysis usable. The constraint forces the spectrum onto the unit circle and fixes a single controlled scale, |λmax|≈1| _ |≈ 1, about which the error dynamics linearize cleanly. The direct-step models have no such control, their spectrum spans a wide range with |λmax|>1| _ |>1, and no comparable linear law holds for them over the rollout. The same holds when the constraint is imposed with an implicit rather than an explicit integrator. Three of the models are constrained with an implicit Euler scheme, two of them trained and emulated at Δt=5×10−2 t=5× 10^-2 and Δt=10−1 t=10^-1. At every Δt t tested, the implicit constraint gives a |λmax|| _ | closer to unity and a lower rate of error growth than the corresponding explicit constraint, and the ordering of |λmax|| _ | reproduces the ordering of the measured RMSE. This transfers the classical result that implicit integrators are more stable at large Δt t to the learned setting, and it follows from the inverse structure of the implicit Jacobian derived in SI, Section 2.2. SI, Section 2 reports the full comparison. The largest eigenvalue is therefore an a priori diagnostic of model behavior, obtained by differentiating the model output with respect to its input through automatic differentiation, without running a full rollout. A larger |λmax|| _ | implies faster error growth, hence poorer short-term accuracy and weaker long-term stability. This ranking holds both across the two model classes and within the integration-constrained class, where the small differences in |λmax|| _ | order the models by error growth along the lines of Fig. 3. Section 3.3 quantifies this predictive power with a linearized error law, built from |λmax|| _ | and the measured one-step error, that reproduces the measured rollout error across the explicit integration-constrained suite (Fig. 4(c,d)). Because |λmax|| _ | is a per-step multiplier, models trained at different time steps are compared through the growth rate per unit time σ=ln⁡|λmax|/Δtσ= | _ |/ t. Measured this way, every integration-constrained model satisfies σ≤0.34σ≤ 0.34, while the two direct-step models give σ=48.1σ=48.1 and σ=60.2σ=60.2. The two classes are separated by more than two orders of magnitude in the rate at which they amplify perturbations. The leading Lyapunov exponent of the KS attractor at this domain length is of order 10−110^-1, so the integration-constrained models amplify perturbations at a rate of the same order as the dynamics, while the direct-step models exceed it by a factor of roughly 500. Table S1 reports σ for every model and SI, Section 1.1 gives the definition. For high-dimensional systems the full eigendecomposition of J can be prohibitive, but the clustering of eigenvalues near one gives a cheap estimate of |λmax|| _ |. Starting from ∑i=1dλi=Tr(), _i=1^d _i=Tr(J), (17) and the triangle inequality |Tr()|=|∑i=1dλi|≤∑i=1d|λi|,|Tr(J)|= | _i=1^d _i |≤ _i=1^d | _i |, (18) the clustering λ1≈λ2≈⋯≈λd≈λmax _1≈ _2≈·s≈ _d≈ _ for integration-constrained models (Fig. 2(c,d), Eq. (10)) turns Eq. (18) into an approximate equality, so |λmax|≈1d|Tr()|,| _ |≈ 1d\,|Tr(J)|, (19) with d the dimension of J. The trace can in turn be estimated with computationally tractable random sketching methods [Meyer et al. 2021], avoiding the full spectrum. 3.3 Deviation from linear stability theory We now revisit the linear scaling law of Section 3.2 and examine where the eigenvalue-only view breaks down, that is, where a model with a larger |λmax|| _ | nonetheless exhibits lower error growth than a model with a smaller one. Figure 4(a,b) compares two FNO models constrained with the same PEC4 integrator, trained with and without the spectral regularizer μ⁡(θ)μ(θ). To the precision of the legend the two models share |λmax|=1.0002| _ |=1.0002, with Table S1 separating them only in the fifth decimal place, so the eigenvalue magnitude alone cannot account for the difference in their error growth; yet the spectrally regularized model has the visibly smaller error over the early rollout in Fig. 4(a). Equation (14) identifies two mechanisms behind this deviation. The first is the additive contribution of the generalization error, ∑iϵ⁡(T(ti)) _i ε(u_T(t_i)). The spectral regularizer, by construction, suppresses the high-wavenumber component of the one-step error and thereby reduces this sum, shown as the dashed lines in Fig. 4(a). Over the early rollout the running sum tracks the total error for each model, confirming that the accumulated generalization error is the leading term of Eq. (14); at late times the total error saturates while the running sum continues to accumulate, marking where the small-error approximation underlying the linear theory ceases to hold. The second mechanism is that the amplification depends on the direction of the error and not only on the eigenvalue magnitudes. The Jacobian acts on the vector ⁡(t)e(t), so the realized growth ‖/‖\|Je\|/\|e\| is a weighted average of the |λi|| _i| over the components of e in the eigenbasis of J, and it equals |λmax|| _ | only when e is aligned with the leading eigenvector. When the eigenvalues λi _i are tightly clustered around unity, as is characteristic of integration-constrained models, the growth of ⁡(t)e(t) is determined not only by |λmax|| _ | but also by the relative orientation (cosine similarity) between ⁡(t)e(t) and the eigenvectors iv_i of J. In Fig. 4(b) each eigenvalue is shaded by |cos⁡(,i)|| (e,v_i)|: the largest-magnitude eigenvalues of the spectrally regularized model carry very low cosine similarity, so they contribute little to the realized amplification, and the marginally larger |λmax|| _ | does not translate into larger ‖⁡(t+Δt)‖\|e(t+ t)\|. SI, Section 4 shows that these direction-dependent effects act within the ⁡(Δt)O( t) residue of the leading-order eigenvalue approximation, so they are finite-Δt t corrections to the law rather than a failure of it. Both mechanisms can be folded into a refined error approximation that remains a priori. Iterating the linearized recursion of Eq. (11) from the initial condition, with the injected one-step error held at its measured norm ‖ϵ‖\| ε\| and the per-step amplification set to the single controlled scale a=|λmax|a=| _ |, gives the geometric law e^(k)=‖ϵ‖ak−1a−1, e(k)=\| ε\|\, a^k-1a-1, (20) where k counts steps from the initial condition, tk=t0+kΔt_k=t_0+k t, so that e^(k) e(k) is the predicted error after k rollout steps. Note that this differs from the lead time of Fig. 3, which is measured from the first step t1t_1: the reference slope there at lead time k is p=k+1p=k+1, whereas here e^(k)→k‖ϵ‖ e(k)→ k\,\| ε\| as a→1a→ 1, recovering the linear scaling law of Eq. (16) with p=kp=k, and departing from it as a moves away from one. Both inputs, a and ‖ϵ‖\| ε\|, are measured at the initial condition, so Eq. (20) predicts the rollout error without running the rollout. Figure 4(c) tests this prediction after k=10k=10 rollout steps for every explicit integration-constrained model: the suite clusters on the line y=xy=x, showing that the two-term law of Eq. (14)—one-step injected error amplified by the Jacobian—predicts error growth across architectures and integration schemes with a single relation. At k=10k=10 and a=1.0002a=1.0002 the geometric factor equals 10.00910.009, within 0.1%0.1\% of k, so this panel tests the injected-error term with the amplification factor contributing negligibly. The amplification factor becomes discriminating at longer lead times and for the direct-step models, where a departs from unity by five percent. The same law also localizes the origin of spectral bias. Restricting the comparison to the high-wavenumber band (ω≥ωcω≥ _c, with ωc=100 _c=100), replacing the injection term by β, the high-ω content of the one-step error, and keeping the same amplification a=|λmax|a=| _ |, again places the suite on y=xy=x (Fig. 4(d)). The high-ω error after k rollout steps is therefore set by the high-ω content of the injected one-step error rather than by any preferential amplification of high wavenumbers, identifying the injected error spectrum ϵ^(ω) ε(ω) as the origin of spectral-bias growth. The Fourier spectra of the predicted state and of its time derivative are shown in SI, Fig. S3, and SI, Section 3 examines how |λmax|| _ | relates to the implicit diffusivity of the model. Deviations from the classical eigenvalue-only prediction therefore emerge from two complementary mechanisms, the direction dependence of error amplification, governed by the eigenvector projections of J, and the cumulative influence of the summed generalization errors, which the spectral regularizer reduces. Once the injected error is measured rather than assumed, the linearized law of Eq. (20) recovers its predictive power across the full explicit integration-constrained suite, and a complete characterization of rollout error growth must account for both the amplification (spectral) and injection (generalization-error) components of Eq. (7). 3.4 Novel stability-promoting loss function The insight about short-term performance obtained from Section 3.3 allows us to develop a novel loss function (see Eq. (38) in Section 6.7) that uses the linearized model to minimize ⁡(t+Δt)e(t+ t) using the projection of ⁡(t)e(t) on J. This constrains both the error as well as the eigenvalues to promote a more stable model. We evaluate the new loss on the FNO-based models constrained with each of the three explicit integrators (Euler, RK4, and PEC4), comparing each model trained with the baseline RMSE loss against its counterpart trained with the proposed Jacobian-based loss. In Fig. 5(a), we verify that training and inference with the new loss follow the scaling obtained in Section 3.2, and that |λmax|| _ | diagnoses model performance: for every integrator, the Jacobian-based loss yields a smaller |λmax|‖⁡(t)‖| _ |\,\|e(t)\| and a correspondingly smaller ‖⁡(t+Δt)‖\|e(t+ t)\| than the RMSE baseline. Figure 5(b) shows the RMSE over the autoregressive rollout for the same six models. To quantify the rate of error growth, the legend reports α, the slope of a linear fit of log10⁡(RMSE) _10(RMSE) against log10 _10 of the timestep over the first 100100 rollout steps, i.e., the exponent of the power law RMSE∝tαRMSE t^α. For every integrator, the model trained with the Jacobian-based loss has a smaller α than its RMSE-trained counterpart (Euler: 0.9160.916 versus 0.9560.956; PEC4: 0.9270.927 versus 0.9460.946; RK4: 0.9650.965 versus 0.9980.998), demonstrating that the new loss slows the growth of error during rollout relative to the basic RMSE loss. A reduction in the one-step generalization error alone cannot produce this effect: the leading term of Eq. (14) is the running sum of ϵ ε, so a uniformly smaller ϵ ε lowers the log-RMSE curve without changing its slope, provided the one-step error evolves similarly along the trajectory, as in the constant-injection law of Eq. (20). The smaller α therefore isolates a reduction in the per-step amplification of the accumulated error, the propagated-error term of Eq. (6). Consistent with Section 3.3, this reduction need not appear in |λmax|| _ |, which bounds but does not set the realized amplification: Table S1 shows that the Jacobian-based loss leaves the FNO |λmax|| _ | essentially unchanged, so for these models the suppression acts through the orientation of the accumulated error relative to the strongly amplified eigendirections of J, whereas for the MLP models the loss also reduces |λmax|| _ | directly. This reduction in the growth rate is an important result: it shows that the Jacobian-driven amplification of error, which the RMSE objective alone leaves uncontrolled, can be regularized directly during training. The reduced growth rate does not by itself guarantee long-term stability however. Several models trained with the Jacobian-based loss (in particular the MLP-based ones) still drift to instability over the longest rollouts because their one-step generalization error ϵ ε remains large, consistent with Eq. (14), which requires both a near-neutral spectrum and a sufficiently small ϵ ε for long-term stability. 4 Conclusion Autoregressive neural emulators are currently trained on a single time step and then used over many time steps: the training objective focuses on the one-step error, while the quantity of interest at inference is the error after hundreds or thousands of steps; and nothing in the training procedure constrains how a one-step error is carried forward. We have developed a theoretical framework to address this fundamental disconnect. At the core of our approach is a single object, the Jacobian of the learned one-step update map taken with respect to the system state and evaluated after training, Eq. (4). Our approach is distinct from the parameter gradients used during training; and, surprisingly, it is rarely examined when these learned models are evaluated. Linearizing the error dynamics about the true trajectory separates the two contributions to emulation error, Eq. (7): the one-step generalization error sets the rate at which error is injected, and the error propagator J sets the rate at which injected error is amplified into the propagated error. Training controls the generalization error and leaves the propagated error free, which is why models with comparable one-step accuracy diverge at different rates during rollout. The magnitude of the largest eigenvalue of the Jacobian, |λmax|| _ |, which governs the growth of a perturbation over one step, is an a priori diagnostic of that amplification. It requires one automatic-differentiation evaluation at a single state and no rollout. Its definition does not reference the architecture, the integration scheme, or the loss function, and across the suite tested here we observe that its predictive value does not depend on them either. Emulators fall into two classes according to how the update map is built, and the classes differ in whether the amplification is controlled at all. A direct-step model predicts the next state in one shot, so its Jacobian is that of the network itself, Eq. (9), which carries no dependence on the time step and places no constraint on the spectrum. An integration-constrained model instead uses the network N to estimate the time derivative and a numerical integrator H to advance it, so the Jacobian takes the form +Δt∇I+ t\, , with H acting on the output of N, Eq. (10). Its eigenvalues then lie within ⁡(Δt)O( t) of unity by construction. The implicit form of the same constraint, (−Δt∇)−1(I- t\, )^-1, places an eigenvalue outside the unit circle only along directions in which the learned tangent operator is expanding. Converted to a growth rate per unit time, which is the comparison that is meaningful across models trained at different time steps, the direct-step models amplify perturbations by orders of magnitude faster than the integration-constrained models. This separation, rather than the small difference in the per-step eigenvalue, is what the measured rollouts reflect, and it accounts for the divergence of direct-step architectures reported previously [Chattopadhyay et al. 2024, Chattopadhyay et al. 2023b, Bi et al. 2023, Pathak et al. 2022, Jiang et al. 2026]. The near-neutral spectrum of the integration-constrained models also yields a scaling law, Eq. (16). The amplification per step lies within ⁡(Δt)O( t) of unity by the structure of Eq. (10), and the injected error is nearly unchanged over the early rollout, which follows from the continuity of ~ G over a small step and which the measured rollouts confirm. The accumulated error after p steps is then p times the first-step error, so the error grows linearly in the step count rather than geometrically in the step count, which is what it does once |λmax|| _ | is bounded away from unity. The law holds across architectures, integration schemes, and loss functions, within the window in which it is derived. A near-neutral spectrum is necessary for slow error growth, but it is not sufficient for long-term stability. Several integration-constrained models drift over the longest rollouts because their accumulated one-step error is large, consistent with Eq. (14). Two further effects limit the eigenvalue-only view. The realized amplification is a weighted average of the eigenvalue magnitudes over the components of the accumulated error in the eigenbasis of the Jacobian, so it depends on the orientation of that error relative to the eigenvectors, which the eigenvalues bound but do not set. The high-wavenumber error after k steps is fixed by the high-wavenumber content of the injected one-step error rather than by preferential amplification of high wavenumbers, which locates the origin of spectral-bias growth in the injected error spectrum rather than in the rollout. Measuring the injected error and retaining |λmax|| _ | as the amplification scale recovers the prediction, Eq. (20). The same analysis yields a training objective, Eq. (38), that regularizes the Jacobian-driven amplification directly and reduces the measured error-growth exponent for every integrator tested. The eigenspectrum of the integration-constrained operators differs from that of a numerical solver for the same equations. The eigenvalues cluster near unity rather than decaying, which makes |λmax|| _ | cheap to estimate and makes the linearization tractable about a single controlled scale. The same clustering is why the eigenvector projections matter, since no direction is strongly preferred by magnitude alone. The direct-step Jacobians do decay, with most of their spectrum massed near the origin and a few eigenvalues beyond the unit circle, which is why neither the single-scale linearization nor the trace estimate applies to them. Our theory has been developed and tested on the Kuramoto–Sivashinsky equation, a canonical system for multi-scale chaotic dynamics, and one on which a factorial sweep over architecture, integrator, and loss is feasible. The Jacobian is evaluated at the initial condition, and although the reported values are robust over 100 random initial conditions, the diagnostic describes the local tangent map rather than its evolution along the attractor. The scaling law holds only for small time steps and over rollouts short enough that the accumulated error remains small. Extending the analysis to the evolution of the Jacobian during rollout, and from the linear theory to a nonlinear one through random matrix theory [Liao and Mahoney 2025, Hodgkinson et al. 2025, Liao and Mahoney 2021, Pennington et al. 2018] and iterated discrete recurrence relation theory [Hodgkinson and Mahoney 2020], are natural next steps, as is testing the diagnostic on higher-dimensional emulators. Numerical analysis provides established tools for deciding, before a simulation is run, whether a discretization of a differential equation will be stable, how fast its error will grow, and how that behavior depends on the scheme and the time step. Scientific machine learning has no comparable body of analysis, and stabilization strategies for learned emulators are instead assembled from heuristics and validated by long rollouts [Lippe et al. 2023, Stachenfeld et al. 2021, Wikner et al. 2022]. The framework developed here is a step toward closing that gap. It takes a trained network as a discrete dynamical system, identifies the operator that governs its behavior during inference, and derives from it the quantities that classical stability analysis provides for numerical schemes. Interpretability work for scientific neural networks has instead centered on saliency maps, layer-wise relevance propagation, and related inspection of internal structure [Molnar 2020]. The diagnostic developed here does not inspect the network, which is what makes it agnostic to architecture and loss, and diagnostics of this kind have been productive outside scientific machine learning [Martin and Mahoney 2019, Martin and Mahoney 2020, Martin et al. 2021]. 5 Acknowledgements AC, PH, and MWM designed the research. CA wrote the computational codes and executed the research. Early versions of these codes were developed by AC. AC and CA wrote the paper, and all the authors analyzed the results and edited the manuscript. AC and CA acknowledge the support from the National Science Foundation (grant no. 2425667), DARPA APaQus program, the Sloan Foundation, and Schmidt Sciences, LLC. PH is grateful to the National Science Foundation (AGS-531264) and Schmidt Sciences, LLC. MWM acknowledges DARPA, NSF, the DOE Competitive Portfolios grant, and the DOE SciGPT grant. Computational resources were provided by NSF ACCESS MTH240019, MTH250006 and NCAR CISL UCSC0008, and UCSC0009. The computational codes can be found on GitHub: https://github.com/cainsliewastaken/Spectral_Stability. 6 Methods and Systems 6.1 Kuramoto–Sivashinsky system We conduct our data-driven experiments on a high-dimensional multi-scale chaotic dynamical system, given by the Kuramoto–Sivashinsky (KS) equations. The governing equations of the KS system are given by: ∂u⁡(x,t)∂t+u⁡(x,t)∂u⁡(x,t)∂x+∂4u⁡(x,t)∂x4+∂2u⁡(x,t)∂x2=0. ∂ u(x,t)∂ t+u(x,t) ∂ u(x,t)∂ x+ ∂^4u(x,t)∂ x^4+ ∂^2u(x,t)∂ x^2=0. (21) We use a domain length (L=100L=100; x∈[−50,50]x∈ [-50,50 ]) to obtain direct numerical simulations (DNS) of the spatio-temporal evolution of this system using an initial condition given by: u⁡(x,0)=−cos⁡(2πxL)(1−sin⁡(2πxL)),u(x,0)=- \! ( 2π xL ) (1- \! ( 2π xL ) ), (22) and periodic boundary conditions. We use 10241024 spatial grid points and a ΔtDNS=10−3 t_DNS=10^-3. We train the suite of autoregressive models on 150000150000 temporal samples of u⁡(x,t)u(x,t) with Δt=ΔtDNS t= t_DNS. We then conduct emulation with the autoregressive models using a new initial condition outside the training dataset for another consecutive 100000100000 time steps. Depending on the model, such emulation will either become unstable, unphysical, or, in certain cases, remain stable and physically consistent. 6.2 Autoregressive models In order to build our suite of autoregressive models, we assume that the underlying dynamical system represented by the KS equations is unknown and the discrete dynamics are represented by: ⁡(x,t+Δt)≈⁡(x,t)+∫t+Δt⁡(⁡(x,t))dt⏟⁡[∘,θ]⏟⁡[∘].u(x,t+ t)≈ u(x,t)+ _t^t+ t F (u(x,t) )dt_N[ ,θ]_H[ ]. (23) In Eq. (23), we describe the general formulation of our autoregressive model, where N is a neural network, and where H is generally a numerical integration method implemented inside a differentiable layer. Most autoregressive models directly represent ⁡[∘]H[ ] with a deep neural network [Pathak et al. 2022, Lam et al. 2022, Bi et al. 2023] (direct step model, from now on); while, based on our previous studies [Chattopadhyay et al. 2023b, Chattopadhyay et al. 2024, Guan et al. 2025], we introduce a numerical integration-based hard constraint that parameterizes ⁡(u⁡(x,t))F (u(x,t) ) with a deep neural network and an explicitly-implemented differentiable layer to perform higher-order integration in ⁡[∘]H[ ]. For either the direct step model or the models with integration-based hard constraint, we train the parameters, θ, of the model by minimizing the difference between the predicted and true value of ⁡(x,t+Δt)u(x,t+ t) at each time step in the training dataset. The general loss function, without any spectral regularizer (see Section 6.6), is given by the standard root-mean-square error (RMSE): L⁡(θ)=1T∑t=0t=T‖⁡(x,t+Δt)−⁡(x,t)−Δt[⁡(⁡(x,t),θ)]‖2,L(θ)= 1TΣ^t=T_t=0||u(x,t+ t)-u(x,t)- tH [N (u(x,t),θ ) ]||_2, (24) where T is the total number of time steps (temporal samples) on which the model is trained. Next, we describe the few different explicit schemes and the one implicit time integration scheme that has been used, and how it has been implemented in our data-driven autoregressive models. 6.3 Explicit integration-based hard constraints Here, we describe the few explicit integration-based hard constraints that we have used in this paper. We have used a first-order Euler integrator, an RK4 integrator, and a 4th4^th-order PEC integrator. For the explicit schemes, the Δt t used to train the models and run inference is equal to ΔtDNS t_DNS, except for the large-Δt t comparison against the implicit scheme in SI, Section 2.3. For the Euler integration scheme, the predicted value of ⁡(x,t+Δt)u(x,t+ t) from Eq. (23) is written as: ⁡(x,t+Δt)=⁡(x,t)+⁡[∘,θ]Δt.u(x,t+ t)=u(x,t)+N[ ,θ] t. (25) The Euler integrator is similar to learning with the residues [Chen and Xiu 2021], which has been used in several different studies prior to this (although, in most cases, the Δt t term is neglected, which results in a non-convergent integration scheme [Krishnapriyan et al. 2023]). The integration-based constraints used in this study are a generalization of learning the residues using rigorous higher-order integrators. For the RK4 integrator, we extend Eq. (25) as: k1 k_1 =⁡[⁡(x,t),θ] =N [u(x,t),θ ] (26) k2 k_2 =⁡[⁡(x,t)+Δt2k1,θ] =N [u(x,t)+ t2k_1,θ ] (27) k3 k_3 =⁡[⁡(x,t)+Δt2k2,θ] =N [u(x,t)+ t2k_2,θ ] (28) k4 k_4 =⁡[⁡(x,t)+Δtk3,θ] =N [u(x,t)+ tk_3,θ ] (29) ⁡(x,t+Δt) (x,t+ t) =⁡(x,t)+Δt16(k1+2k2+2k3+k4)⏟⁡[∘]. =u(x,t)+ t 16 (k_1+2k_2+2k_3+k_4 )_H[ ]. (30) For the 4th4^th-order PEC method (PEC4), we use the following formulation: ~(x,t+Δt) u(x,t+ t) =⁡(x,t)+Δt[⁡(x,t),θ] =u(x,t)+ tN[u(x,t),θ] (31) ^(x,t+Δt) u(x,t+ t) =⁡(x,t)+Δt2(⁡[⁡(x,t),θ]+⁡[~(x,t+Δt),θ]) =u(x,t)+ t2 (N [u(x,t),θ ]+N [ u(x,t+ t),θ ] ) (32) ¯(x,t+Δt) u(x,t+ t) =⁡(x,t)+Δt2(⁡[⁡(x,t),θ]+⁡[^(x,t+Δt),θ]) =u(x,t)+ t2 (N [u(x,t),θ ]+N [ u(x,t+ t),θ ] ) (33) ⁡(x,t+Δt) (x,t+ t) =⁡(x,t)+Δt12(⁡[⁡(x,t),θ]+⁡[¯(x,t+Δt),θ])⏟⁡[∘]. =u(x,t)+ t 12 (N [u(x,t),θ ]+N [ u(x,t+ t),θ ] )_H [ ]. (34) For higher-order integrators such as RK4 or PEC4, we represent the predicted ⁡(x,t+Δt)=⁡(x,t)+Δt[⁡(⁡(x,t),θ)]u(x,t+ t)=u(x,t)+ tH [N (u(x,t),θ ) ], where H encapsulates the scheme as a differentiable operator. 6.4 Implicit integration-based hard constraint Implicit time-integration schemes have improved stability properties for larger Δt t, as compared with explicit integration schemes [Frank et al. 1997]. In order to implement implicit time-integration as a hard constraint, we use the theory of implicit layers [Kawaguchi 2021] in deep learning. In this approach, we can use any higher-order scheme as our ⁡[∘]H [ ] operator. The equation involving the time stepper is given by: ⁡(x,t+Δt)⏟∗=⁡(x,t)+Δt[⁡[⁡(x,t+Δt)⏟∗,θ]], u(x,t+ t)_y^*=u(x,t)+ tH [N [ u(x,t+ t)_y^*,θ ] ], (35) where ∗y^* is required to be solved via an implicit layer at each epoch (i.e. each value of θ) in the training process. In order to do that, we initialize ∗=⁡(x,t)y^*=u (x,t ), and we execute a fixed-point iteration to converge to the correct ⁡(x,t+Δt)u(x,t+ t). More details about implicit layers and deep equilibrium models can be found in Kawaguchi 2021, Bai et al. 2019. For the experiments with the implicit method, the chosen Δt t ranges from ΔtDNS t_DNS up to 100ΔtDNS100\, t_DNS (SI, Section 2.3). Here, we outline the fixed point iteration in Algorithm 1 used to train (and run inference with) the implicit numerical integrator. Differentiating Eq. (35) with respect to the state gives =(−Δt∇[(∗;θ)])−1J=(I- t\, [N(y^*;θ)])^-1, which reduces to Eq. (10) to first order in Δt t. SI, Section 2.2 gives the derivation and the resulting expansion conditions. Algorithm 1 Implicit time integration scheme 1: Input: ϵε ⊳ Convergence threshold, typically 10−810^-8 2: At each epoch, fix θ 3: Initialize 0∗=⁡(x,t)y^*_0=u(x,t), k=0k=0 4: repeat 5: k=k+1k=k+1 6: k∗=⁡(x,t)+Δt[⁡(k−1∗,θ)]y^*_k=u(x,t)+ t\,H [N (y^*_k-1,θ ) ] 7: until ‖k∗−k−1∗‖2≤ϵ||y^*_k-y^*_k-1||_2≤ε 8: Output: k∗y^*_k ⊳ k∗y^*_k is the predicted value of ⁡(x,t+Δt)u(x,t+ t) 6.5 Neural architectures We have used an MLP and an FNO to represent N. The MLP has 66 layers and 20002000 neurons in each layer with a ReLU activation. The FNO that has been used has 66 Fourier blocks and retains 512512 modes in each Fourier block along with a width of 3232. We have conducted extensive hyperparameter trials before selecting these sets of hyperparameters for our empirical analysis. Each hyperparameter trial optimizes the short-term accuracy of the model. 6.6 Spectral bias and spectral regularizer Generally, to train autoregressive models, previous studies have used a traditional L2L_2 based loss function, optimizing the difference between the true ⁡(x,t+Δt)u(x,t+ t) and its prediction, ⁡[⁡(⁡(x,t))]H [N (u(x,t) ) ], using Eq. (24). However, for multi-scale dynamical systems, e.g., turbulent flows, a spectral bias [Chattopadhyay et al. 2023b, Cao et al. 2019, Wang and Lai 2024] prevents the network from learning the high-wavenumber component of the dynamics. This often leads to instability and physical inconsistency in the models. Furthermore, spectral bias does not only present itself in the Fourier spectrum of the predicted state, but also in its derivative with respect to time. In this paper, the time-derivative is represented by ⁡(⁡(x,t))N (u(x,t) ) for the explicit methods. We propose a spectral regularizer to the L2L_2 loss function presented in Eq. (24) as: μ(θ)=1T∑t=0T∥∂^(x,t)∂t−^[[(x,t)]]∥2,μ (θ )= 1TΣ^T_t=0 ∂ u(x,t)∂ t- H [N [u(x,t) ] ] _2, (36) where [∘] [ ] represents the 1D Fourier transform in x. To approximate the true and predicted time derivatives from the training and predicted data, we use a first-order finite difference. The total loss, LtotalL_total, is given by: Ltotal=L⁡(θ)+γμ(θ),L_total=L(θ)+γμ(θ), (37) where γ is the Lagrange multiplier set to 0.10 after significant hyperparameter trials. SI, Section 3 reports the effect of this regularizer on the Fourier spectrum of the emulation. 6.7 Stability-promoting loss function Based on the insight from linear stability analysis as seen in Eq. (10), we propose a new stability-promoting loss function for the autoregressive models with integration constraints that minimizes ‖⁡(t+Δt)‖2||e(t+ t)||_2 in its linearized form, for each value of t in the training dataset: Ltotal L_total =1T∑t=0T∥((t+Δt)−(t)−Δt[[(t),θ]])∥2 = 1T _t=0^T (u(t+ t)-u(t)- t\,H\! [N\! [u(t),θ ] ] ) _2 +μ(θ)+1T∑t=0T∥(t+Δt)−(t)−Δt[[(t),θ]]∥2, + \ μ(θ) \ + \ 1T _t=0^T (t+ t)-u(t)- t\,H\! [N\! [u(t),θ ] ] _2, (38) where J is defined in Eq. (10). This novel loss function minimizes the error at the next time step, which is approximated by projecting the predicted state at the next time step on the Jacobian of the model. A more aggressive loss function that focuses more on stability could have used the projection of p(t+Δt)u_p(t+ t) on the eigenvector of J corresponding to the largest eigenvalue. Such an approximation would have been accurate for a J with a decaying eigen spectrum; but the spectrum of J for the models with integration constraints is not decaying, instead having eigenvalues that are clustered around 11 on the unit circle (see Fig. 2(c,d) in Section 3). As such, projecting J on the eigenvector corresponding to the largest eigenvalue is often a poor approximation of next time-step error. Models trained with the Jacobian term but without the spectral regularizer μ⁡(θ)μ(θ) are listed as “Jacobian” in Table S1, and those trained with both terms as “Spectral+Jac.”. Figure 2: Short-term error, rollout error accumulation, and eigenanalysis of the models. Throughout, direct-step models are shown in vermilion and integration-constrained models in blue, with MLP as solid lines and FNO as dotted lines. (a) RMSE over the first 100100 steps. The two direct-step models (Direct MLP, Direct FNO) grow rapidly and saturate near the climatological decorrelation level, where the RMSE is ⁡(1)O(1) in these normalized units. The two integration-constrained Euler models (Euler MLP, Euler FNO), whose eigenvalue structure is representative of all 2727 integration-constrained models examined here, remain stable over this window, with error orders of magnitude smaller and indistinguishable from zero on the linear axis. (b) The same RMSE on log–log axes over the first 10410^4 steps. The direct-step models saturate within a few tens of steps (MLP) to a few hundred (FNO), whereas the Euler models grow slowly and in proportion to the step count, remaining orders of magnitude more accurate over most of the window and reaching the decorrelation level only at its right edge; over the full emulation several MLP-based models eventually drift as their accumulated generalization error grows. (c) Eigenvalues of the model Jacobian J at the initial condition t0t_0 on the Argand plane for all four models, with |λmax|| _ | in the legend. The Direct MLP (∘ ) and Direct FNO (□ ) admit eigenvalues with |λmax|>1| _ |>1 (1.04931.0493 and 1.06201.0620) spread well inside and beyond the unit circle, whereas the Euler MLP (⋆ ) and Euler FNO (△ ) eigenvalues collapse onto the unit circle (dotted) and cluster tightly at λ=1λ=1 (|λmax|=1.0002| _ |=1.0002), appearing as a single dense cluster. (d) Zoomed-in view of the integration-constrained eigenvalues from (c), unit circle again dotted. Each spectrum spreads only in the third decimal place, and the two models separate in the fourth, with |λmax|=1.000175| _ |=1.000175 for the Euler MLP and 1.0002011.000201 for the Euler FNO. This neutral spectrum, consistent across architecture, explains the slow accumulation in proportion to the step count in (a) and (b) through the Jacobian structure of Eq. (10). SI, Fig. S1 shows the equivalent curves and spectra for the remaining explicit models, and Table S1 reports |λmax|| _ | for every model. Figure 3: A single linear scaling law governs error growth across the full model suite. Norm of the accumulated error ‖⁡(t+kΔt)‖2\|e(t+k t)\|_2 plotted against |λmax|‖⁡(t)‖2| _ |\,\|e(t)\|_2 for all 2929 models, where λmax _ is the largest eigenvalue of the model Jacobian J evaluated at the initial condition and ∥⋅∥2\|·\|_2 is the Euclidean (L2L_2) norm of the error vector. Here t=t1t=t_1, so ⁡(t)=⁡(t1)=ϵ⁡(t0)e(t)=e(t_1)= ε(t_0) is the first-step (generalization) error and the lead time k is measured from the first step; the reference slope at lead time k is therefore the step count p=k+1p=k+1 of Eq. (16). Panels (a)–(d) show the successive lead times k=1,10,100,1000k=1,10,100,1000, one per panel, on log–log axes. Marker color denotes model class (direct-step in vermilion, integration-constrained in blue) and fill denotes architecture (MLP filled, FNO open), while marker shape indicates the integration scheme (Euler, RK4, PEC4, implicit Euler). The gray reference line in each panel is the linear scaling law of Eq. (16) evaluated at that lead time; on log–log axes it is a unit-slope line whose vertical offset grows with the step count, from slope 22 in (a) to slopes 1111, 101101, and 10011001 in (b)–(d). The panels correspond to pΔt=0.002p\, t=0.002, 0.0110.011, 0.1010.101, and 1.0011.001, so they span the validity window of Eq. (16) from well inside it to its boundary. The integration-constrained models (|λmax|≈1| _ |≈ 1) lie on the predicted line in (a) and (b), scatter modestly about it in (c), and depart from it in both directions in (d), where pΔt≈1p\, t≈ 1 and the derivation no longer applies. The two direct-step models sit at much larger |λmax|‖⁡(t)‖2| _ |\,\|e(t)\|_2 and rise above the reference line at the lead time of (b), reflecting the faster-than-linear growth driven by their |λmax|>1| _ |>1; by (c) and (d) they have reached the climatological decorrelation level and stopped growing, so they fall below a reference line that continues to grow in proportion to p. This faster-than-linear growth arises from the ⁡(‖⁡(t)‖22)O(\|e(t)\|_2^2) higher-order terms in the Taylor expansion of ~ G, which are neglected in the linear analysis but become non-negligible once the accumulated error grows large. Each model is identified in Table S1. Two of the implicit-Euler models are trained at Δt≫ΔtDNS t t_DNS and therefore exceed the pΔt≪1p\, t 1 window by a larger margin at the longest lead time (SI, Section 2.4). Figure 4: Error growth is set by the accumulated generalization error and the error–eigenvector projection, not by |λmax|| _ | alone, and is predicted by the linearized error law. (a, b) Two FNO models constrained with the same PEC4 integrator, trained without (FNO, blue) and with (Spectral FNO, green) the spectral regularizer μ⁡(θ)μ(θ) of Eq. (36). (a) RMSE over the autoregressive rollout on log–log axes. For each model the solid line is the total error ⁡(u⁡(t))e(u(t)) and the black dashed line is the cumulative generalization error ∑i=1tϵ⁡(u⁡(i)) _i=1^t ε(u(i)), the leading term of Eq. (14), with the dash length distinguishing the two models. The spectrally regularized model has the smaller error over the early rollout. (b) Eigenvalues of the model Jacobian near λ=1λ=1 for the FNO (∘ ) and Spectral FNO (⋆ ), with the unit circle dotted; to the precision of the legend the two models share |λmax|=1.0002| _ |=1.0002. Each marker is shaded by the cosine similarity |cos⁡(,i)|| (e,v_i)| between its eigenvector iv_i and the error vector ⁡(t)e(t) (grayscale bar). (c, d) The linearized error law of Eq. (20), evaluated with a=|λmax|a=| _ | from the Jacobian at the initial condition, tested after k=10k=10 rollout steps (tk=t0+kΔt_k=t_0+k t) for every explicit integration-constrained model (marker shape: integration scheme, Euler / RK4 / PEC4; fill: MLP filled, FNO open); each point is one model and the gray line is y=xy=x. (c) Total error, with ‖ϵ‖\|ε\| the measured one-step error. (d) The same law restricted to the high-wavenumber band (ω≥ωcω≥ _c, with ωc=100 _c=100), with the injection term replaced by β, the high-ω content of the one-step error. See Section 3.3 for the analysis. The two models in (a, b) are models 11 and 12 of Table S1, and the models in (c, d) are models 1 to 24. Figure 5: The stability-promoting loss lowers error growth. FNO-based models constrained with the Euler, PEC4, and RK4 integrators, each trained with the baseline RMSE loss and compared against its counterpart trained with the proposed Jacobian-based stability-promoting loss of Eq. (38). (a) The linear scaling law of Fig. 3, ‖⁡(t+Δt)‖\|e(t+ t)\| versus |λmax|‖⁡(t)‖| _ |\,\|e(t)\|, with the RMSE baselines shown as circles and the stability-promoting loss as stars. All points lie on the predicted unit-slope line of offset 22 (gray), confirming that the new loss preserves the linear scaling law, and for every integrator the stability-promoting loss yields a smaller |λmax|‖⁡(t)‖| _ |\,\|e(t)\| and correspondingly smaller ‖⁡(t+Δt)‖\|e(t+ t)\| than the baseline. (b) RMSE over the autoregressive rollout for the same models, with the RMSE baselines shown as solid lines and the stability-promoting loss as dotted lines. The legend reports α, the slope of a linear fit of log10⁡(RMSE) _10(RMSE) against log10 _10 of the timestep over the first 100100 rollout steps (the exponent of RMSE∝tαRMSE t^α). For every integrator, the stability-promoting loss produces a smaller α, and hence lower error growth, than the corresponding RMSE baseline, because it constrains the Jacobian-driven amplification of error, which the RMSE objective alone leaves uncontrolled. The six models shown are models 7, 9, 11, 19, 21, and 23 of Table S1. Supporting Information for Eigenanalysis framework for autoregressive neural emulators of multi-scale chaotic dynamics Conrad Ainslie, Pedram Hassanzadeh, Michael W. Mahoney, and Ashesh Chattopadhyay E-mail: aschatto@ucsc.edu This PDF file includes: Supporting text, Sections 1 to 4 Figs. S1 to S3 Table S1 Algorithm S1 SI References Supporting Information Text This document contains four sets of supporting results. Section 1 reports the full 29-model suite together with its measured Jacobian spectra, and introduces the time-normalized growth rate used to compare models trained at different time steps. Section 2 describes the implicit integration-based hard constraint, derives the Jacobian of the implicit update map, and reports the comparison against explicit schemes over a range of Δt t. Section 3 presents the Fourier-space analysis of spectral bias. Section 4 establishes the near-normality of the integration-constrained Jacobian and the resulting validity of the eigenvalue approximation used in Section 2.4 of the main text. 1 The full model suite Table S1 lists all 29 models evaluated in this work. Two are direct-step models and 27 are integration-constrained. Figure S1 demonstrates the RMSE growth and the Jacobian eigenspectra of the 24 explicit integration-constrained models in a single view, providing the complete-suite reference behind the representative subsets shown in Figs. 2 and 3 of the main text. 1.1 Time-normalized growth rate The largest eigenvalue |λmax|| _ | is a per-step quantity, so models trained at different time steps cannot be ranked by |λmax|| _ | alone. Interpreting the learned update map as a time-Δt t flow map, |λmax|=exp⁡(σΔt)| _ |= (σ\, t) defines a growth rate per unit time, σ=ln⁡|λmax|Δt,σ\;=\; | _ | t, (S1) which is comparable across Δt t and carries the units of a Lyapunov exponent. Table S1 reports σ alongside |λmax|| _ |. The normalization separates the two model classes by more than two orders of magnitude. Every integration-constrained model has σ≤0.34σ≤ 0.34, and the majority lie in the range 0.010.01 to 0.250.25. The two direct-step models have σ=48.1σ=48.1 and σ=60.2σ=60.2. The direct-step models therefore amplify perturbations at a rate between 140140 and 180180 times larger than the fastest integration-constrained model, and this separation, rather than the difference in the per-step eigenvalue, is what the rollouts of Fig. 2(a,b) of the main text reflect. The integration-constrained values are of order 10−110^-1 per unit time, which is the order of the leading Lyapunov exponent of the KS attractor, so these models amplify perturbations at a rate set by the dynamics rather than by the emulator. The quantity σ inherits the precision of |λmax|| _ |. At Δt=10−3 t=10^-3 the five-decimal values of Table S1 resolve σ to one or two significant figures only, so small differences in σ within the integration-constrained class should not be over-interpreted. Equation (S1) also uses |λmax|| _ | and therefore describes the modulus of the dominant multiplier, not the full complex spectrum. Table S1: The full model suite. All 29 models evaluated in this work, comprising 27 integration-constrained models and 2 direct-step models. Every explicit model is trained and run at Δt=ΔtDNS=10−3 t= t_DNS=10^-3; the two large-Δt t implicit models (26, 27) are the exception. |λmax|(t0)| _ |(t_0) is the largest-magnitude eigenvalue of the model Jacobian at the initial condition. σ=ln⁡|λmax|/Δtσ= | _ |/ t is the corresponding growth rate per unit time of Eq. (S1), which is comparable across Δt t; its precision is limited by the five decimal places reported for |λmax|| _ |. RMSE is the short-term error at k=10k=10 steps. The spectral-bias metric is ⟨log10⁡(|u^pred|/|u^true|)⟩ω≥100 _10(| u_pred|/| u_true|) _ω≥ 100 at step 100, with positive values indicating excess high-wavenumber energy relative to truth. # Arch. Scheme Loss Δt t |λmax|(t0)| _ |(t_0) σ RMSE (k=10k=10) Spec. bias 1 MLP Euler RMSE 10−310^-3 1.000 181.000\,18 0.1800.180 1.289×10−031.289×10^-03 12.4612.46 2 MLP Euler Spectral 10−310^-3 1.000 021.000\,02 0.0200.020 3.850×10−043.850×10^-04 10.4610.46 3 MLP RK4 RMSE 10−310^-3 1.000 211.000\,21 0.2100.210 1.310×10−031.310×10^-03 12.4712.47 4 MLP RK4 Spectral 10−310^-3 1.000 071.000\,07 0.0700.070 3.970×10−043.970×10^-04 10.4410.44 5 MLP PEC4 RMSE 10−310^-3 1.000 181.000\,18 0.1800.180 1.273×10−031.273×10^-03 12.4912.49 6 MLP PEC4 Spectral 10−310^-3 1.000 301.000\,30 0.3000.300 3.577×10−043.577×10^-04 10.3510.35 7 FNO Euler RMSE 10−310^-3 1.000 201.000\,20 0.2000.200 1.142×10−041.142×10^-04 9.829.82 8 FNO Euler Spectral 10−310^-3 1.000 191.000\,19 0.1900.190 2.526×10−052.526×10^-05 8.798.79 9 FNO RK4 RMSE 10−310^-3 1.000 341.000\,34 0.3400.340 2.343×10−042.343×10^-04 9.879.87 10 FNO RK4 Spectral 10−310^-3 1.000 201.000\,20 0.2000.200 1.035×10−041.035×10^-04 8.888.88 11 FNO PEC4 RMSE 10−310^-3 1.000 191.000\,19 0.1900.190 1.153×10−041.153×10^-04 9.559.55 12 FNO PEC4 Spectral 10−310^-3 1.000 201.000\,20 0.2000.200 2.214×10−052.214×10^-05 8.818.81 13 MLP Euler Jacobian 10−310^-3 1.000 011.000\,01 0.0100.010 3.850×10−043.850×10^-04 10.8610.86 14 MLP Euler Spectral+Jac. 10−310^-3 1.000 051.000\,05 0.0500.050 8.121×10−058.121×10^-05 11.0111.01 15 MLP PEC4 Jacobian 10−310^-3 1.000 021.000\,02 0.0200.020 3.739×10−043.739×10^-04 11.2311.23 16 MLP PEC4 Spectral+Jac. 10−310^-3 1.000 081.000\,08 0.0800.080 2.335×10−042.335×10^-04 11.0311.03 17 MLP RK4 Jacobian 10−310^-3 1.000 011.000\,01 0.0100.010 3.935×10−043.935×10^-04 10.9610.96 18 MLP RK4 Spectral+Jac. 10−310^-3 1.000 011.000\,01 0.0100.010 5.071×10−055.071×10^-05 10.9510.95 19 FNO Euler Jacobian 10−310^-3 1.000 211.000\,21 0.2100.210 1.794×10−051.794×10^-05 9.839.83 20 FNO Euler Spectral+Jac. 10−310^-3 1.000 211.000\,21 0.2100.210 9.718×10−069.718×10^-06 9.809.80 21 FNO PEC4 Jacobian 10−310^-3 1.000 221.000\,22 0.2200.220 2.977×10−052.977×10^-05 9.939.93 22 FNO PEC4 Spectral+Jac. 10−310^-3 1.000 251.000\,25 0.2500.250 3.381×10−053.381×10^-05 10.2710.27 23 FNO RK4 Jacobian 10−310^-3 1.000 231.000\,23 0.2300.230 2.482×10−052.482×10^-05 10.0010.00 24 FNO RK4 Spectral+Jac. 10−310^-3 1.000 211.000\,21 0.2100.210 4.453×10−054.453×10^-05 10.2610.26 25 FNO Implicit Euler RMSE 10−310^-3 1.000 101.000\,10 0.1000.100 4.131×10−044.131×10^-04 10.2110.21 26 FNO Implicit Euler RMSE 5×10−25× 10^-2 1.000 201.000\,20 0.0040.004 1.876×10−021.876×10^-02 10.0810.08 27 FNO Implicit Euler RMSE 10−110^-1 1.018 901.018\,90 0.1870.187 3.840×10−023.840×10^-02 10.4810.48 28 MLP Direct RMSE 10−310^-3 1.049 301.049\,30 48.12348.123 5.909×10−015.909×10^-01 14.1714.17 29 FNO Direct RMSE 10−310^-3 1.062 001.062\,00 60.15460.154 3.433×10−023.433×10^-02 11.1211.12 Figure S1: Short-term error and eigenspectra of the full explicit integration-constrained suite. All 24 explicit integration-constrained models of Table S1, spanning the Euler, RK4, and PEC4 schemes, the MLP and FNO architectures, and the RMSE, spectral, Jacobian, and spectral+Jacobian losses. Every model is trained and emulated at Δt=ΔtDNS=10−3 t= t_DNS=10^-3, with MLP models drawn as solid lines and FNO models as dotted lines, and color and marker distinguishing individual models as indicated in the legend. (a) RMSE over the first 100 rollout steps on log–log axes. Every model exhibits the slow, near-linear error accumulation predicted by Eq. (16) of the main text; the curves separate in vertical offset, set by the one-step generalization error ϵ ε, rather than in slope, consistent with the near-neutral spectra shared across the suite. The window shown satisfies pΔt≤0.1p\, t≤ 0.1, so it lies inside the validity range of that scaling law. (b) Eigenvalues of the model Jacobian J, evaluated at the initial condition t0t_0, on the Argand plane, with the unit circle dotted. For every model the eigenvalues collapse onto the unit circle and cluster tightly at λ=1λ=1, with the spread of each spectrum confined to the third decimal place and the models separated in the fourth (Table S1 reports |λmax|(t0)| _ |(t_0) for each). This confirms that the qualitative spectrum is insensitive to the architecture, the explicit integration scheme, and the loss function, as expected from the Jacobian structure of Eq. (10) of the main text. 2 Implicit integration-based hard constraints 2.1 Implementation Implicit time-integration schemes have improved stability properties for larger Δt t, as compared with explicit integration schemes [Frank et al. 1997]. In order to implement implicit time-integration as a hard constraint, we use the theory of implicit layers [Kawaguchi 2021] in deep learning. In this approach, we can use any higher-order scheme as our ⁡[∘]H [ ] operator. The equation involving the time stepper is given by ⁡(x,t+Δt)⏟∗=⁡(x,t)+Δt[⁡[⁡(x,t+Δt)⏟∗,θ]], u(x,t+ t)_y^*=u(x,t)+ t\,H [N [ u(x,t+ t)_y^*,θ ] ], (S2) where ∗y^* is required to be solved via an implicit layer at each epoch, that is, at each value of θ, in the training process. In order to do that, we initialize ∗=⁡(x,t)y^*=u (x,t ), and we execute a fixed-point iteration to converge to the correct ⁡(x,t+Δt)u(x,t+ t). More details about implicit layers and deep equilibrium models can be found in Kawaguchi 2021, Bai et al. 2019. For the experiments with the implicit method, the chosen Δt t ranges from ΔtDNS t_DNS up to 100ΔtDNS100\, t_DNS. Algorithm S1 outlines the fixed-point iteration used to train and to run inference with the implicit numerical integrator. Algorithm S1 Implicit time integration scheme 1: Input: ϵε ⊳ Convergence threshold, typically 10−810^-8 2: At each epoch, fix θ 3: Initialize 0∗=⁡(x,t)y^*_0=u(x,t), k=0k=0 4: repeat 5: k=k+1k=k+1 6: k∗=⁡(x,t)+Δt[⁡(k−1∗,θ)]y^*_k=u(x,t)+ t\,H [N (y^*_k-1,θ ) ] 7: until ‖k∗−k−1∗‖2≤ϵ||y^*_k-y^*_k-1||_2≤ε 8: Output: k∗y^*_k ⊳ k∗y^*_k is the predicted value of ⁡(x,t+Δt)u(x,t+ t) 2.2 The Jacobian of the implicit update map The implicit models do not share the Jacobian of the explicitly constrained models, and the difference is what produces their distinct eigenspectrum. Differentiating Eq. (S2) with respect to the state Tu_T, and using the fact that ∗y^* itself depends on Tu_T, gives ∂∗∂T=+Δt∇[(∗;θ)]∂∗∂T, ^* _T=I+ t\, \! [N(y^*;θ) ] ^* _T, (S3) so that the Jacobian of the implicit update map is =(−Δt∇[(∗;θ)])−1.J= (I- t\, \! [N(y^*;θ) ] )^-1. (S4) Expanding Eq. (S4) in a Neumann series gives =+Δt∇[(∗;θ)]+(Δt2)J=I+ t\, [N(y^*;θ)]+O( t^2), so the implicit Jacobian agrees with Eq. (10) of the main text to first order in Δt t. The conclusion about explicit schemes carries over to the implicit models, and the implicit models are integration-constrained in that sense. The distinction matters only at ⁡(Δt2)O( t^2) and above, which is the regime probed by the large-Δt t experiments of Section 2.4. The inverse structure of Eq. (S4) controls the spectrum. Let νi _i denote the eigenvalues of ∇[(∗,θ)] [N(y^*;θ)], which carry the units of a rate. The eigenvalues of the implicit Jacobian are λi=11−Δtνi, _i= 11- t\, _i, (S5) whereas the corresponding explicit scheme yields λi=1+Δtνi _i=1+ t\, _i. Writing νi=ai+ibi _i=a_i+i\,b_i and evaluating the moduli exactly gives the two conditions for local expansion, |λi|>1⟺ai>Δt2|νi|2(implicit),|λi|>1⟺ai>−Δt2|νi|2(explicit).| _i|>1\; \;a_i> t2\,| _i|^2 (implicit), | _i|>1\; \;a_i>- t2\,| _i|^2 (explicit). (S6) The implicit scheme places an eigenvalue outside the unit circle only when the learned tangent operator has an expanding direction, ai>0a_i>0. The explicit scheme places an eigenvalue outside the unit circle even for mildly contracting directions with −Δt2|νi|2<ai<0- t2| _i|^2<a_i<0. At fixed Δt t, the implicit spectrum therefore lies closer to and inside the unit disk than the explicit one. This is the origin of the smaller |λmax|| _ | and the longer left tail of the implicit eigenspectra in Fig. S2(d), and it is the spectral counterpart of the classical statement that implicit schemes have larger stability regions. 2.3 Short-term and long-term performance From classical numerical analysis, we know that implicit numerical time integrators, while more computationally expensive, have larger stability regions than explicit time integration schemes, and in practice are more accurate at large Δt t for stiff problems. We use our proposed a priori diagnostic metric, |λmax|| _ |, to investigate whether similar properties emerge for neural autoregressive models. We consider models that are hard-constrained with implicit and explicit integration schemes and are trained and tested with Δt t ranging from ΔtDNS t_DNS up to 100ΔtDNS100 t_DNS. To do so, we have performed experiments with Δt=ΔtDNS=10−3 t= t_DNS=10^-3, Δt=50ΔtDNS=5×10−2 t=50 t_DNS=5× 10^-2, and Δt=100ΔtDNS=10−1 t=100 t_DNS=10^-1. Unless otherwise stated, all models in this work are trained and emulated at Δt=ΔtDNS t= t_DNS. The large-Δt t models of this section are the exception, and the explicit models trained at Δt>ΔtDNS t> t_DNS are used only for this comparison and are not part of the 29-model suite of Table S1. As shown in Fig. S2(a), similar to the explicit integration-based constraint with the Euler scheme, the implicit Euler scheme-based model is also long-term stable and physically consistent. However, the RMSE error curves in Fig. S2(c) clearly demonstrate that for all values of Δt t, the implicit time integration constraint-based models have lower error growth, as compared to explicit ones. This is expected from traditional numerical analysis. This is further validated by investigating |λmax|| _ | of each of the implicit schemes and confirming that they are smaller, that is, closer to unity, than the explicit schemes. Similar to explicit constraints, for implicit constraints as well, |λmax|| _ | acts as an a priori diagnostic of the error-growth rate of the models. This is of particular importance when developing novel autoregressive architectures for complex systems. The integration constraints allow us to translate our knowledge about the performance and stability of different integration schemes to neural network-based modeling in scientific computing. The comparison across Δt t requires the normalization of Eq. (S1), since |λmax|| _ | is a per-step multiplier and the three experiments use different step sizes. Expressed as growth rates per unit time, the implicit constraint gives σ=0.100σ=0.100, 0.0040.004, and 0.1870.187 at Δt=10−3 t=10^-3, 5×10−25× 10^-2, and 10−110^-1, against σ=0.200σ=0.200, 0.1550.155, and 0.2430.243 for the corresponding explicit models. The implicit scheme has the lower rate at every Δt t tested. The raw |λmax|| _ | values reported in the legend of Fig. S2(c) order the implicit and explicit models correctly within each Δt t pair, but they cannot be compared across pairs, and the normalized rate is what makes the three comparisons commensurate. 2.4 Large time steps, expanding eigenvalues, and stable emulation Model 27 of Table S1 is an implicit Euler FNO trained and emulated at Δt=10−1 t=10^-1. Its Jacobian has |λmax|=1.0189| _ |=1.0189, which places an eigenvalue outside the unit circle, and the linear analysis of Section 2.3 of the main text therefore predicts local expansion. The model nonetheless remains stable and physically consistent over the full emulation. The quantity |λmax|| _ | is a multiplier per time step, and models trained at different Δt t take different numbers of steps to cover the same physical time, so their per-step values cannot be compared with one another. Equation (S1) converts |λmax|| _ | to a growth rate per unit time, which is comparable. Model 27 gives σ=0.187σ=0.187. The two direct-step models give σ=48.1σ=48.1 and σ=60.2σ=60.2, larger by factors of roughly 260260 and 320320, even though their per-step eigenvalues, 1.0491.049 and 1.0621.062, are numerically close to 1.0191.019. The closeness of the three raw eigenvalues follows from the difference in Δt t and carries no information about the models. A rate of σ≈0.19σ≈ 0.19 is what an emulator of the KS attractor should give. The attractor has positive Lyapunov exponents, so a tangent operator that reproduces it must have expanding directions, and an emulator whose spectrum lay strictly inside the unit circle would be over-diffusive rather than accurate. Equation (S5) converts |λmax|=1.0189| _ |=1.0189 at Δt=10−1 t=10^-1 into an expanding eigenvalue of the learned tangent operator of 0.1850.185 per unit time, which agrees with σ=0.187σ=0.187 to within one percent, the residual difference being of order Δt t and expected from the two definitions. The expansion described by |λmax|>1| _ |>1 is local and is bounded by nonlinear saturation on the attractor, so it does not imply unbounded divergence, and the emulation confirms this. Sections 2.3 and 2.4 of the main text assume Δt≪1 t 1. At Δt=10−1 t=10^-1 that condition is marginal, the ⁡(Δt2)O( t^2) terms neglected in Section 4 are no longer small, and the quantitative predictions of the linear theory carry a correspondingly wider tolerance for this model. Equation (16) of the main text holds for pΔt≪1p\, t 1, and this restriction governs where the implicit models sit in Fig. 3 of the main text. At the longest lead time, p=1001p=1001, models 26 and 27 reach pΔt=50p\, t=50 and pΔt=100p\, t=100, so the condition fails by two orders of magnitude. Their departure from the reference line in that panel is what the stated validity window predicts and is not evidence against the scaling law. Figure S2: Short-term and long-term performance of neural autoregressive models constrained with implicit time integration. Color encodes the time step, Δt=0.1 t=0.1 in blue, 0.050.05 in vermilion, 0.0010.001 in black. Implicit schemes are solid lines and filled pentagons, explicit schemes dotted lines and open stars. (a) Long-term emulation of the KS state u⁡(x,t)u(x,t) by the FNO model constrained with first-order implicit Euler at Δt=0.05=50ΔtDNS t=0.05=50\, t_DNS, stable and physically plausible over 10510^5 steps of length ΔtDNS t_DNS. (b) Pointwise difference between the emulation in (a) and the reference numerical solution. The error stays small and develops slowly, with structure inherited from the flow rather than diverging, which quantifies the physical consistency. (c) RMSE growth on linear axes for the FNO models with implicit and explicit Euler at the three time steps. At each Δt t the implicit scheme gives lower error growth, and the linear scale makes the separation visible late in the rollout, where accumulated error dominates. The legend reports |λmax|| _ |, closer to unity for the implicit scheme in each pair. Within each Δt t pair the ordering of |λmax|| _ | matches the ordering of the RMSE curves, so |λmax|| _ | acts as an a priori diagnostic of the error-growth rate. Comparison across pairs requires the normalized rate σ of Eq. (S1), reported in Section 2.3. (d) Eigenvalues on the Argand plane for both schemes at Δt=0.1 t=0.1 and 0.050.05, with the unit circle dotted. At each Δt t the implicit spectrum has a smaller maximum and a more diffusive left tail than the explicit one, consistent with its lower error growth. Equation (S6) accounts for both. The implicit map places an eigenvalue outside the unit circle only for expanding directions of the learned tangent operator, the explicit map even for mildly contracting ones. 3 Spectral bias and the Fourier spectrum of the emulation This section presents the Fourier-space analysis of spectral bias, connecting the eigenvalues of the model Jacobian to the fidelity of the predicted spectrum and its time derivative. We analyze whether |λmax|| _ | can explain different values of spectral bias in neural autoregressive models. As shown previously [Chattopadhyay et al. 2023b, Guan et al. 2025], a spectral regularizer in the form of Eq. (36) of the main text can reduce spectral bias in neural autoregressive models and improve stability. However, for a decaying spectrum like the one in the KS system, the regularizer does not completely eliminate spectral bias either in the state or in its derivative (Fig. S3(b,c)). This is because the amplitude of the modes of the state in a system with a decaying spectrum decreases as a power law. Figure S3(b) shows that after wavenumber 100100, the slope of the spectrum increases with diminished amplitude of the high-wavenumber modes. Despite the spectral regularizer, spectra with high slopes are more challenging to capture, an effect mirrored in the spectra of state derivatives (Fig. S3(c)) as well. A detailed explanation of the impact of the slope of the system spectrum on the spectral bias of the autoregressive model can be found in Chattopadhyay et al. 2023b. Nonetheless, the extent of spectral bias varies across architectures. For instance, neural operators have been shown to better capture the spatial Fourier spectra of chaotic systems [Azizzadenesheli et al. 2024, Oommen et al. 2025]. The magnitude of |λmax|| _ |, representing the distance of the dominant Jacobian eigenvalue from the origin, reflects the implicit diffusivity of a model. Smaller values indicate higher diffusion and often correlate with reduced spectral bias. Figure S3(c) shows that the FNO model equipped with the PEC-based integrator and a spectral regularizer has lower spectral bias, and hence more diffusion, than the MLP-based model. This is captured by the smaller |λmax|| _ | of the FNO model as compared to the MLP model, as shown in Fig. S3(a). The quantity |λmax|| _ | therefore complements the one-step error as an a priori indicator. It sets the amplification scale, while ϵ ε sets the injection. For the matched pair studied here it also tracks the implicit diffusivity and the fidelity of the high-wavenumber spectrum. This correspondence is established for a single matched pair of models and should be read as an empirical association rather than as a general law. Section 3.3 of the main text establishes the more general statement, namely that the high-wavenumber error after k rollout steps is set by the high-ω content of the injected one-step error rather than by preferential amplification of high wavenumbers. Figure S3: Spectral bias and its connection to the eigenvalues of the model Jacobian. (a) Eigenvalues on the Argand plane for the MLP-based model (vermilion plus signs, |λmax|=1.0003| _ |=1.0003) and the FNO-based model (blue open stars, |λmax|=1.0002| _ |=1.0002), both constrained with the PEC4 integrator and trained with the spectral regularizer of Eq. (36) of the main text, with the unit circle dotted. The FNO-based model has the smaller |λmax|| _ | and a longer left tail, that is, eigenvalues further inside the unit circle, indicating greater implicit diffusivity. This is reflected in (b) the spatial Fourier spectrum |u^(ω)|| u(ω)| of the predicted state after 100100 time steps and (c) the spectrum of its time derivative d|u^|/dtd| u|/dt. Up to ω≈100ω≈ 100 both models track the reference numerical data (dotted black). In the shaded high-wavenumber band (ω≳100ω 100), where the true spectrum continues to decay, both learned models instead saturate at a flat noise floor set by the smallest scale each can represent rather than by the true dynamics. The floor is markedly lower for the FNO (blue) than for the MLP (vermilion) in both the state and its derivative, so the FNO more faithfully captures the high-wavenumber tail. For this matched pair, |λmax|| _ | thus tracks the implicit diffusivity and the high-ω noise floor, with the smaller |λmax|| _ | corresponding to greater implicit diffusion and more faithful capture of the high-wavenumber tail of the spectrum. 4 Near-normality of the integration-constrained Jacobian and validity of the eigenvalue approximation The scaling arguments of Section 2.4 of the main text replace the action of the Jacobian on the error vector by its largest eigenvalue, ‖‖≈|λmax|‖\|J\,e\|≈| _ |\,\|e\|. For a general non-normal matrix this replacement can fail. Eigenvalues bound only the asymptotic growth of perturbations, while transient, orientation-dependent amplification is controlled by the singular values, which for a strongly non-normal matrix can exceed |λmax|| _ | substantially. Here we show that the integration-constrained structure, together with the smallness of Δt t, renders the constrained Jacobian near-normal, so that the eigenvalue approximation is accurate to the same order as the theory itself. This argument applies only to the integration-constrained models. No analogous control is available for the direct-step Jacobian. Write the constrained Jacobian of Eq. (10) of the main text as =+ΔtJ=I+ t\,A, with =∇[⁡(T,θ)]A= [N(u_T;θ)]. Its departure from normality is measured by the self-commutator, ∗−∗=Δt2(∗−∗),JJ^*-J^*J\;=\; t^2 (AA^*-A^*A ), (S7) since the identity commutes with every matrix and the terms linear in Δt t cancel. The departure from normality is therefore ⁡(Δt2)O( t^2), one order smaller than the ⁡(Δt)O( t) corrections already retained in the scaling analysis. To first order in Δt t the constrained Jacobian is normal, regardless of how non-normal the network Jacobian A itself may be. The same structure controls the singular values directly. The squared singular values of J are the eigenvalues of ∗=+Δt(+∗)+Δt2∗,J^*J\;=\;I+ t (A+A^* )+ t^2\,A^*A, (S8) so that, for bounded ‖\|A\|, every singular value satisfies σi()=1+⁡(Δt) _i(J)=1+O( t). Consequently, for an error vector of any orientation, ‖‖∈[σmin(),σmax()]= 1+⁡(Δt),while|λmax|=1+⁡(Δt), \|J\,e\|\|e\|\;∈\; [ _ (J),\, _ (J) ]\;=\;1+O( t), | _ |=1+O( t), (S9) so that ‖‖=|λmax|‖+⁡(Δt)‖\|J\,e\|=| _ |\,\|e\|+O( t)\,\|e\|. Replacing the vectorial action of J by the scalar |λmax|| _ | is therefore exact at leading order, and the error incurred is of the same ⁡(Δt)O( t) size as the corrections already neglected in Section 2.4 of the main text. At Δt=ΔtDNS=10−3 t= t_DNS=10^-3 this residue is small, consistent with the fourth-decimal spread of the measured spectra in Fig. 2(d) of the main text and in Table S1. The symbol σi _i in Eqs. (S8) and (S9) denotes a singular value of J and is distinct from the growth rate σ of Eq. (S1). It is this near-identity structure, rather than the clustering of the eigenvalues per se, that justifies the directional insensitivity invoked in Section 2.4 of the main text. Clustered eigenvalues alone would not constrain the transient amplification of a strongly non-normal matrix. The ⁡(Δt)O( t) residue left open by the leading-order approximation is the range within which the direction-dependent effects of Section 3.3 of the main text act. The error-eigenvector projections and the accumulated one-step error modulate the realized amplification within the ⁡(Δt)O( t) band around unity. The deviations documented in Section 3.3 of the main text are thus the expected finite-Δt t corrections to the leading-order law, rather than a breakdown of the framework. For the direct-step models, by contrast, the Jacobian ∇(T,θ) (u_T;θ) carries no Δt t scaling, its departure from normality is uncontrolled, and eigenvalue-based estimates should be read only as bounds on asymptotic growth. The same Δt t scaling delimits where the argument applies. At Δt=10−1 t=10^-1, used for model 27 of Table S1, the ⁡(Δt2)O( t^2) departure from normality is four orders of magnitude larger than at Δt=ΔtDNS t= t_DNS, and the near-normality guarantee is correspondingly weaker. The implicit Jacobian of Eq. (S4) admits the same expansion, =+Δt+⁡(Δt2)J=I+ t\,A+O( t^2) with =∇[⁡(∗,θ)]A= [N(y^*;θ)], so Eqs. (S7) to (S9) hold for the implicit models at leading order as well. References Azizzadenesheli et al. [2024] K. Azizzadenesheli, N. Kovachki, Z. Li, M. Liu-Schiaffini, J. Kossaifi, and A. Anandkumar. Neural operators for accelerating scientific simulations and design. Nature Reviews Physics, 6(5):320–328, 2024. Bai et al. [2019] S. Bai, J. Z. Kolter, and V. Koltun. Deep equilibrium models. Advances in neural information processing systems, 32, 2019. Bi et al. [2023] K. Bi, L. Xie, H. Zhang, X. Chen, X. Gu, and Q. Tian. Accurate medium-range global weather forecasting with 3D neural networks. Nature, pages 1–6, 2023. Cao et al. [2019] Y. Cao, Z. Fang, Y. Wu, D.-X. Zhou, and Q. Gu. Towards understanding the spectral bias of deep learning. arXiv preprint arXiv:1912.01198, 2019. Chattopadhyay et al. [2022] A. Chattopadhyay, M. Mustafa, P. Hassanzadeh, E. Bach, and K. Kashinath. Towards physics-inspired data-driven weather forecasting: integrating data assimilation with a deep spatial-transformer-based U-NET in a case study with ERA5. Geoscientific Model Development, 15(5):2221–2237, 2022. Chattopadhyay et al. [2023a] A. Chattopadhyay, J. Pathak, E. Nabizadeh, W. Bhimji, and P. Hassanzadeh. Long-term stability and generalization of observationally-constrained stochastic data-driven models for geophysical turbulence. Environmental Data Science, 2:e1, 2023a. Chattopadhyay et al. [2023b] A. Chattopadhyay, Y. Q. Sun, and P. Hassanzadeh. Challenges of learning multi-scale dynamics with ai weather models: Implications for stability and one solution. arXiv e-prints, pages arXiv–2304, 2023b. Chattopadhyay et al. [2024] A. Chattopadhyay, M. Gray, T. Wu, A. B. Lowe, and R. He. Oceannet: a principled neural operator-based digital twin for regional oceans. Scientific Reports, 14(1):21181, 2024. Chen and Xiu [2021] Z. Chen and D. Xiu. On generalized residual network for deep learning of unknown dynamical systems. Journal of Computational Physics, 438:110362, 2021. Fan et al. [2020] H. Fan, J. Jiang, C. Zhang, X. Wang, and Y.-C. Lai. Long-term prediction of chaotic systems with machine learning. Physical Review Research, 2(1):012080, 2020. Floryan [2024] D. Floryan. On instabilities in neural network-based physics simulators. arXiv preprint arXiv:2406.13101, 2024. Frank et al. [1997] J. Frank, W. Hundsdorfer, and J. G. Verwer. On the stability of implicit-explicit linear multistep methods. Applied Numerical Mathematics, 25(2-3):193–205, 1997. Guan et al. [2025] H. Guan, T. Arcomano, A. Chattopadhyay, and R. Maulik. Lucie: A lightweight uncoupled climate emulator with long-term stability and physical consistency. Journal of Advances in Modeling Earth Systems, 17(11):e2025MS005152, 2025. Hodgkinson and Mahoney [2020] L. Hodgkinson and M. W. Mahoney. Multiplicative noise and heavy tails in stochastic optimization. Technical Report Preprint: arXiv:2006.06293, 2020. Hodgkinson et al. [2025] L. Hodgkinson, Z. Wang, and M. W. Mahoney. Models of heavy-tailed mechanistic universality. Technical Report Preprint: arXiv:2506.03470, 2025. Jiang et al. [2026] R. Jiang, X. Zhang, K. Jakhar, P. Y. Lu, P. Hassanzadeh, M. Maire, and R. Willett. Hierarchical implicit neural emulators. Advances in Neural Information Processing Systems, 38:73718–73751, 2026. Kawaguchi [2021] K. Kawaguchi. On the theory of implicit deep learning: Global convergence with implicit layers. arXiv preprint arXiv:2102.07346, 2021. Keisler [2022] R. Keisler. Forecasting global weather with graph neural networks. arXiv preprint arXiv:2202.07575, 2022. Krishnapriyan et al. [2021] A. S. Krishnapriyan, A. Gholami, S. Zhe, R. M. Kirby, and M. W. Mahoney. Characterizing possible failure modes in physics-informed neural networks. Technical Report Preprint: arXiv:2109.01050, 2021. Krishnapriyan et al. [2023] A. S. Krishnapriyan, A. F. Queiruga, N. B. Erichson, and M. W. Mahoney. Learning continuous models for continuous physics. Communications Physics, 6(1):319, 2023. Lai et al. [2025] C.-Y. Lai, P. Hassanzadeh, A. Sheshadri, M. Sonnewald, R. Ferrari, and V. Balaji. Machine learning for climate physics and simulations. Annual Review of Condensed Matter Physics, 16(1):343–365, 2025. Lam et al. [2022] R. Lam, A. Sanchez-Gonzalez, M. Willson, P. Wirnsberger, M. Fortunato, A. Pritzel, S. Ravuri, T. Ewalds, F. Alet, Z. Eaton-Rosen, et al. Graphcast: Learning skillful medium-range global weather forecasting. arXiv preprint arXiv:2212.12794, 2022. Li et al. [2022] Z. Li, M. Liu-Schiaffini, N. Kovachki, K. Azizzadenesheli, B. Liu, K. Bhattacharya, A. Stuart, and A. Anandkumar. Learning chaotic dynamics in dissipative systems. Advances in Neural Information Processing Systems, 35:16768–16781, 2022. Liao and Mahoney [2021] Z. Liao and M. W. Mahoney. Hessian eigenspectra of more realistic nonlinear models. Advances in Neural Information Processing Systems, 34:20104–20117, 2021. Liao and Mahoney [2025] Z. Liao and M. W. Mahoney. Random matrix theory for deep learning: Beyond eigenvalues of linear models. Technical Report Preprint: arXiv:2506.13139, 2025. Lippe et al. [2023] P. Lippe, B. S. Veeling, P. Perdikaris, R. E. Turner, and J. Brandstetter. Pde-refiner: Achieving accurate long rollouts with neural pde solvers. arXiv preprint arXiv:2308.05732, 2023. Lupin-Jimenez et al. [2025] L. Lupin-Jimenez, M. Darman, S. Hazarika, T. Wu, M. Gray, R. He, A. Wong, and A. Chattopadhyay. Simultaneous emulation and downscaling with physically-consistent deep learning-based regional ocean emulators. arXiv preprint arXiv:2501.05058, 2025. Martin and Mahoney [2019] C. H. Martin and M. W. Mahoney. Traditional and heavy-tailed self regularization in neural network models. In Proceedings of the 36th International Conference on Machine Learning, pages 4284–4293, 2019. Martin and Mahoney [2020] C. H. Martin and M. W. Mahoney. Heavy-tailed Universality predicts trends in test accuracies for very large pre-trained deep neural networks. In Proceedings of the 20th SIAM International Conference on Data Mining, 2020. Martin et al. [2021] C. H. Martin, T. Peng, and M. W. Mahoney. Predicting trends in the quality of state-of-the-art neural networks without access to training or testing data. Nature Communications, 12(1):4122, 2021. Meyer et al. [2021] R. A. Meyer, C. Musco, C. Musco, and D. P. Woodruff. Hutch++: Optimal stochastic trace estimation. In Symposium on Simplicity in Algorithms (SOSA), pages 142–155. SIAM, 2021. Molnar [2020] C. Molnar. Interpretable machine learning. Lulu. com, 2020. Oommen et al. [2025] V. Oommen, A. Bora, Z. Zhang, and G. E. Karniadakis. Integrating neural operators with diffusion models improves spectral representation in turbulence modelling. Proceedings of the Royal Society A, 481(2309):20240819, 2025. Pathak et al. [2022] J. Pathak, S. Subramanian, P. Harrington, S. Raja, A. Chattopadhyay, M. Mardani, T. Kurth, D. Hall, Z. Li, K. Azizzadenesheli, et al. FourCastNet: A global data-driven high-resolution weather model using adaptive Fourier neural operators. arXiv preprint arXiv:2202.11214, 2022. Pedersen et al. [2025] C. Pedersen, L. Zanna, and J. Bruna. Thermalizer: Stable autoregressive neural emulation of spatiotemporal chaos. arXiv preprint arXiv:2503.18731, 2025. Pennington et al. [2018] J. Pennington, S. Schoenholz, and S. Ganguli. The emergence of spectral universality in deep networks. In International Conference on Artificial Intelligence and Statistics, pages 1924–1932. PMLR, 2018. Price et al. [2025] I. Price, A. Sanchez-Gonzalez, F. Alet, T. R. Andersson, A. El-Kadi, D. Masters, T. Ewalds, J. Stott, S. Mohamed, P. Battaglia, et al. Probabilistic weather forecasting with machine learning. Nature, 637(8044):84–90, 2025. Sakarvadia et al. [2025] M. Sakarvadia, K. Hegazy, A. Totounferoush, K. Chard, Y. Yang, I. Foster, and M. W. Mahoney. The false promise of zero-shot super-resolution in machine-learned operators. Technical Report Preprint: arXiv:2510.06646, 2025. Sambamurthy and Chattopadhyay [2025] A. Sambamurthy and A. Chattopadhyay. Lazy diffusion: Mitigating spectral collapse in generative diffusion-based stable autoregressive emulation of turbulent flows. arXiv preprint arXiv:2512.09572, 2025. Stachenfeld et al. [2021] K. Stachenfeld, D. B. Fielding, D. Kochkov, M. Cranmer, T. Pfaff, J. Godwin, C. Cui, S. Ho, P. Battaglia, and A. Sanchez-Gonzalez. Learned coarse models for efficient turbulence simulation. arXiv preprint arXiv:2112.15275, 2021. Wang and Lai [2024] Y. Wang and C.-Y. Lai. Multi-stage neural networks: Function approximator of machine precision. Journal of Computational Physics, 504:112865, 2024. Watt-Meyer et al. [2025] O. Watt-Meyer, B. Henn, J. McGibbon, S. K. Clark, A. Kwa, W. A. Perkins, E. Wu, L. Harris, and C. S. Bretherton. Ace2: accurately learning subseasonal to decadal atmospheric variability and forced responses. npj Climate and Atmospheric Science, 8(1):205, 2025. Wikner et al. [2022] A. Wikner, B. R. Hunt, J. Harvey, M. Girvan, and E. Ott. Stabilizing machine learning prediction of dynamics: Noise and noise-inspired regularization. arXiv preprint arXiv:2211.05262, 2022. Yu et al. [2024] A. Yu, D. Lyu, S. H. Lim, M. W. Mahoney, and N. B. Erichson. Tuning frequency bias of state space models. arXiv preprint arXiv:2410.02035, 2024. Yu et al. [2025] A. Yu, D. C. Maddix, B. Han, X. Zhang, A. F. Ansari, O. Shchur, C. Faloutsos, A. G. Wilson, M. W. Mahoney, and Y. Wang. Understanding the implicit biases of design choices for time series foundation models. Technical Report Preprint: arXiv:2510.19236, 2025.

Instant research discovery