A Response Theory Probe for Learned Stochastic AI Simulators, Tested on Lorenz-63

A Response Theory Probe for Learned Stochastic AI Simulators, Tested on Lorenz-63

arXiv:2610.06798v1 [math.DS] 5 Oct 2026

João Böger Department of Technology, Management and Economics Technical University of Denmark Kongens Lyngby, 2800, DK jpade@dtu.dk

Simon Driscoll Department of Applied Mathematics and Theoretical Physics University of Cambridge Cambridge, CB3 0WA, UK sd2136@cam.ac.uk

Valerio Lucarini School of Computing and Mathematical Sciences University of Leicester Leicester, LE1 7RH, UK School of Sciences Great Bay University Dongguan, P.R. China v.lucarini@leicester.ac.uk

Niccolò Zagli School of Computing and Mathematical Sciences University of Leicester Leicester, LE1 7RH, UK nz68@leicester.ac.uk

Francisco Camara Pereira Department of Technology, Management and Economics Technical University of Denmark Kongens Lyngby, 2800, DK camara@dtu.dk

Abstract Machine-learning emulators of chaotic and stochastic systems are usually validated on forecast skill and long-run statistics. Neither certifies that an emulator responds correctly to forcing, the property that projection and attribution studies rely on. Linear response theory makes this testable: the forced response follows from unperturbed correlations through a generalized fluctuation–dissipation relation, and decomposes over the stochastic Ruelle–Pollicott resonances of the Koopman generator. Building on the Koopmanism Response framework, we turn this into a calibrated, mode-resolved test for learned surrogates: each surrogate rollout passes or fails each check, and failure rates are compared with those of independent realizations of the true system. On stochastic Lorenz-63, a three-variable toy model, we evaluate SINDy, an MLP, a reservoir computer, a neural ODE and a neural SDE with learned diffusion, over up to 80 rollouts each. A sparse-regression model with the correct library passes every check at rates consistent with the true system. Invariant-statistics fidelity and response fidelity dissociate in both directions: a quarter of reservoir-computer rollouts pass every invariant-statistics check and match the static susceptibility χ(0), yet misrepresent the slow relaxation modes, while the neural ODE and SDE rarely meet the invariant-statistics floor but recover those modes in three quarters of rollouts. As expected of a time-integrated quantity dominated here by fast relaxation, χ(0) does not separate these cases. For a fixed network, the training formulation (one-step drift, flow map, or multi-step through the integrator) decides which of these properties it gets right.

1

Introduction & Motivation

Machine learning is now a standard tool for modelling dynamical systems, from sparse regression of governing equations [4] to reservoir computing [22] and neural differential equations [8]. Data-driven 40th Conference on Neural Information Processing Systems (NeurIPS 2026). Workshop: AI for Stochastic Dynamics.

weather models such as GraphCast [17], FourCastNet [23], Pangu-Weather [1] and Aurora [2] rival physics-based forecasts at a fraction of the cost. Their validation still rests mainly on forecast skill and long-run statistics, increasingly through benchmarks such as WeatherBench 2 [26]; case studies and dynamical tests of their physical behaviour are emerging [6, 12, 5], but are not yet standard. Both test fidelity to the unperturbed distribution, not how the statistics change under forcing, which is what projection, attribution and the study of tipping points require [11]: a model can reproduce the invariant measure closely and still misrepresent how it deforms under forcing [10]. Response theory makes this testable. For a large class of chaotic and stochastic systems, the linear response is a correlation function of the unperturbed dynamics, in the spirit of the generalized fluctuation–dissipation theorem [21, 27], and it decomposes over the stochastic Ruelle–Pollicott (RP) resonances of the Koopman generator and their eigenfunctions [7, 19]. The Koopmanism Response framework [20, 30] computes these objects from a single unperturbed trajectory with extended dynamic mode decomposition (EDMD), with ResDMD residuals guarding against spurious modes [29, 9], and yields a mode-resolved decomposition of the forced Green’s function [18]. “Did the model learn the physics?” then becomes a falsifiable question: does a surrogate trained only on unperturbed data carry the correct resonances, eigenfunctions and static susceptibility? Stochastic Lorenz-63 is a three-variable toy model of convection with no spatial structure, and our surrogates are not operational climate emulators: it is a controlled testbed with a known true response, used to map how network design and training formulation shape response fidelity, complementing forced-response studies of neural climate emulators [10]. We contribute (i) a calibrated evaluation protocol for learned stochastic simulators, in which invariant statistics, the static susceptibility and a mode-resolved eigenfunction test are each applied to every rollout and compared with the false-alarm rate of independent true realizations; (ii) a formulation ladder (one-step drift regression, flow map, multi-step neural ODE, neural SDE with learned diffusion) sharing architecture and data, which isolates the effects of the training target, the multi-step loss and learned noise, complemented by an exploratory causal transformer, the architecture family of current weather emulators; and (iii) evidence that invariant-statistics fidelity and response fidelity dissociate in both directions, which a scalar response summary such as χ(0) cannot detect here but the mode-resolved test does.

2

Method

System and response instrument. For a dissipative stochastic system at stationarity, the linear response to a small forcing is fixed by the unperturbed dynamics through the RP resonances of the Koopman generator and a generalized fluctuation–dissipation relation [20, 30] (Apps. A, B). Following [20], we study stochastic Lorenz-63, the Itô SDE dx = F (x) dt + ς dWt with (σ, ρ, β) = (10, 28, 8/3) and isotropic noise ς = 2, integrated by Euler–Maruyama at δt = 0.001. We force the x-equation and track x. From one unperturbed trajectory of T = 10000 time units sampled at τ = 0.1, we build a Fourier-EDMD approximation of the generator (dictionary size (2Kmax + 1)3 = 729, TSVD relative threshold 10−4 ), read off the leading resonances, reconstruct the Green’s function Gx (t) with the Method-3 estimator [30], and summarize it by the static susceptibility Z ∞ X α∗ M γk δ⟨x⟩∞ k χ(0) = = τ , (1) Gx (t) dt = ε 1 − νk 0 k

where δ⟨x⟩∞ is the long-term shift of ⟨x⟩, νk are the EDMD eigenvalues and the sum runs over all EDMD modes. Learned surrogates. All models are trained on segments of the same unperturbed master trajectory and rolled out with additive noise (App. C). (i) SINDy [4, 3] with a quadratic library containing the true equations (positive control); (ii) an MLP (2 × 64, tanh) fitted to the drift at ∆tfit = 0.01 (MLP-drift) or to the next state (MLP-flow); (iii) a reservoir computer (RC) [22] fitted as a flow map; (iv) a neural ODE (NODE) [8], the same network trained through a differentiable RK4 step on multi-step windows; (v) a neural SDE (NSDE) [14], the NODE drift plus a learned diagonal diffusion. All rollouts except the NSDE use the true ς, an oracle that favours the surrogates; the NSDE uses no true parameter. Evaluation protocol. Every surrogate passes through the identical instrument and is compared with a held-out true realization by pass/fail checks in four rungs. (1) Invariant statistics: x-marginal 2

(a)

(b)

(c) 20

Z t Gx (s) ds

0.4

Â(t) =

0.05

15

Â(0) band

!q1

I(¸)

0

½(x)

0.10

10

0.2

!p1 5 0

0.00

0.0 ¡20

0

20

0

5

Truth truth 5-95% range

10

t

x SINDy (full data) MLP-drift

MLP-flow RC

NODE NSDE

¡1.0

¡0.5

0.0

R(¸)

reference eigenvalue

Figure 1: Invariant statistics, static susceptibility and resonances, per model. (a) x-marginal density: truth median and 5–95% range over 50 realizations; per model, the rollout with median χ(0) Rt (MLP-drift not shown). (b) Accumulated response χ(t) = 0 Gx (s) ds, which tends to χ(0): truth median and 5–95% range; per model, the median over rollouts; blue band: the χ(0) band of Table 1. (c) Leading Koopman/RP resonances (Im λ ≥ 0): reference (stars) and, per model, the median of each leading group over rollouts, with its 5–95% range in Re λ; real modes spread vertically for visibility. Dashed and dotted lines: the shortest-UPO twirling frequency ωp1 and the intra-lobe turnover frequency ωq1 [20].

Wasserstein-1, x-autocorrelation distance, wing-switch ratio and lobe balance f+ , each of which must lie within the range (“floor”) reached by 50 independent true realizations. (2) Static susceptibility: χ(0) must lie within the range of 100 true realizations. (3–4) Mode-resolved test, scored separately on the two slow real modes and the four oscillatory ones (Fig. 1c): for each leading eigenfunction we measure the alignment of the surrogate’s eigenfunction φ̃ with the reference φ by the overlap O = |⟨φ, φ̃⟩µ0 |/(∥φ∥µ0 ∥φ̃∥µ0 ) ∈ [0, 1], the cosine of the angle between them weighted by the invariant measure (O = 1: same function up to scale). The mode passes if O exceeds what independent true realizations typically reach, their mean overlap minus three standard deviations (clipped to [0.5, 0.99] for the oscillatory modes). A rollout fails a rung if it fails any check in it, and a model’s failure rate on a rung is the fraction of its rollouts that fail. The failure rate of 50 fresh true realizations, scored as a perfect surrogate, is each rung’s false-alarm rate. We call a model’s failure rate consistent with the truth when a Fisher exact test against the true realizations gives p ≥ 0.05. Table 1 reports medians over 8 training seeds × 10 rollout seeds (SINDy: 10 rollouts); per-rollout failure rates are tabulated in App. C.

3

Experiments

The false-alarm rates are low: fresh true realizations fail the invariant statistics in 4% of cases, the slow modes in 6% and the oscillatory modes in 2%, and a true realization left out of the χ(0) band falls outside it in 2% of cases. SINDy trained on the full record recovers the structural coefficients to within 0.116% and fails at rates consistent with the truth on every rung, with no mode pass decided by the 0.99 cap: the instrument certifies a correct model rather than only rejecting wrong ones. The reservoir computer reproduces the invariant measure closely (W1 inside its floor in 76 of 78 rollouts; Fig. 1a), and 26 of its rollouts lie inside all four invariant-statistics floors and the χ(0) band. Of these, 20 fail the slow-mode test, against 3 of the 48 truths that pass the same checks (Fisher p = 4.9 × 10−10 ). The failure is not a disguised switching error: RC fails the slow modes in 70% of rollouts inside the switch-ratio floor and 80% outside it (p = 0.424). Spectrally, RC places its two leading oscillatory resonances within 0.133% of the reference frequencies, but its second slow mode decays faster than in any true realization in 77% of rollouts (median Re λ = −0.639 against the true range [−0.611, −0.483]; Fig. 1c), and its slow eigenfunctions pass in under a third of rollouts. 3

Table 1: Invariant statistics, χ(0) and eigenfunction overlaps, per model. Medians over scored rollouts; χ(0) as mean with the standard deviation in subscript. First row: the floor of each invariant statistic, the χ(0) band and the overlap threshold of each leading mode, labelled by its reference eigenvalue or frequency. Red: the median lies outside it. Bold: the best model in each column (lowest distance; switch ratio closest to one and f+ to one half; χ(0) mean closest to the truth’s; highest overlap). n: scored rollouts (RC: 2 diverged, left out). Slow (Re λ)

Invariant statistics Model

n

Truth

W1

ACF

switch

f+

≤0.284 ≤0.089 0.985–1.035 0.487–0.517

threshold

χ(0)

Oscillatory (ω)

−0.366 −0.508 8.42

9.85

9.98

16.68

0.162–0.309

≥0.945

≥0.922 ≥0.990 ≥0.968 ≥0.960 ≥0.990

50 0.076

0.038

1.011

0.500

0.224(0.034)

0.992

0.991

0.999

0.990

0.991

0.998

0.092 5.599 2.193 0.097 0.337 0.362

0.038 3.724 0.735 0.070 0.093 0.102

1.013 0.008 0.725 1.037 1.049 1.050

0.494 0.477 0.474 0.499 0.504 0.502

0.215(0.029) 1.045(1.278) 0.636(0.313) 0.203(0.032) 0.225(0.045) 0.221(0.039)

0.994 0.043 0.820 0.923 0.976 0.980

0.986 0.319 0.755 0.888 0.962 0.960

0.999 0.607 0.963 0.998 0.998 0.998

0.990 0.374 0.751 0.895 0.952 0.953

0.993 0.290 0.691 0.856 0.935 0.929

0.998 0.115 0.245 0.997 0.996 0.996

SINDy 10 MLP-drift 80 MLP-flow 80 RC 78 NODE 80 NSDE 80

The MLP baselines fail every rung (Table 1). Regressing ∆x/∆t at ∆tfit , the drift MLP learns the O(∆t) secant drift F + (∆t/2)JF F rather than F (its error matches this bias to within 7%), and its Euler rollout collapses onto one wing in 33 of 80 rollouts. Fitted as a flow map, the same network restores two-wing switching but not the measure; its slowest real mode decays too slowly (−0.146 against the true −0.389), so only 33% of its χ(0) has accrued by t = 1 and χ(0) overshoots. Trained through the integrator on multi-step windows, the same network (NODE) and its stochastic extension (NSDE) show the reverse of RC’s pattern: both over-switch in every rollout (median switch ratio 1.05 and 1.05) and fail invariant statistics in 93% and 100% of rollouts, yet the slow modes in only 24% and 23%. The over-switching is learned, not numerical: the true drift rolled out with the same RK4 step under-switches (median ratio 0.953). An exploratory causal transformer fine-tuned on short rollouts, as recent weather emulators are, shows the same pattern (App. C.5). For SINDy, RC, NODE and NSDE the χ(0) failure rate never exceeds 9%, while their slow-mode failure rates range from 0% to 76%. This is expected for a time-integrated quantity: for forcing and observing x, χ(0) is set mostly by fast relaxation: across the true realizations a median 89% of it has accrued by t = 1 (Fig. 1b), before the slow modes (e-folding times 1.97 and 2.73) have relaxed. Forcings or observables coupled to lobe occupation may load on exactly these modes; the resonances and eigenfunctions, unlike χ(0), do not depend on that choice.

4

Discussion

The test is calibrated, not exhaustive: one low-dimensional system, one forcing direction and one observable. Thresholds come from 16 calibration truths; rebuilding the slow-mode thresholds from the 50 fresh truths lowers RC’s slow-mode failure rate from 76% to 55%, and the 0.99 cap decides many oscillatory passes (App. A). Except for the NSDE, rollouts use the true noise intensity; as a robustness check, SINDy fitted to only 100 time units acquires symmetry-breaking terms (App. C). Neither invariant statistics nor χ(0) certifies the slow dynamics: a quarter of RC’s rollouts pass both with the wrong slow modes. A one-step drift regressed at a coarse step builds the O(∆t) secant bias into the model and breaks the attractor. One-step flow maps can match the attractor (RC) yet misplace the slow, weakly damped modes, plausibly because these change little over one step and barely enter a one-step loss. Training on multi-step windows through the integrator (NODE, NSDE) recovers the slow modes in three quarters of rollouts but over-switches; the rollout-fine-tuned transformer shows the same split. Multi-step losses, now common in weather emulators [17, 2], thus appear to carry the slow dynamics, and could be paired with invariant-measure constraints or spectral priors rather than with more flexible architectures. As slow modes govern transitions between the quasi-invariant sets of a metastable system and, for forcings that couple to them, its response, a mode-resolved check, needing one unperturbed trajectory per model, should accompany invariant statistics in validating learned simulators. A natural next step is to test whether these patterns persist in spatially extended systems and for state-of-the-art architectures of current emulators. 4

References [1] Kaifeng Bi, Lingxi Xie, Hengheng Zhang, Xin Chen, Xiaotao Gu, and Qi Tian. Accurate mediumrange global weather forecasting with 3D neural networks. Nature, 619(7970):533–538, July 2023. ISSN 1476-4687. doi: 10.1038/s41586-023-06185-3. URL https://www.nature.com/articles/ s41586-023-06185-3. [2] Cristian Bodnar, Wessel P. Bruinsma, Ana Lucic, Megan Stanley, Anna Allen, Johannes Brandstetter, Patrick Garvan, Maik Riechert, Jonathan A. Weyn, Haiyu Dong, Jayesh K. Gupta, Kit Thambiratnam, Alexander T. Archibald, Chun-Chieh Wu, Elizabeth Heider, Max Welling, Richard E. Turner, and Paris Perdikaris. A foundation model for the Earth system. Nature, 641(8065):1180–1187, May 2025. ISSN 1476-4687. doi: 10.1038/s41586-025-09005-y. URL https://www.nature.com/articles/ s41586-025-09005-y. [3] Lorenzo Boninsegna, Feliks Nüske, and Cecilia Clementi. Sparse learning of stochastic dynamical equations. The Journal of Chemical Physics, 148(24):241723, March 2018. ISSN 0021-9606. doi: 10.1063/1.5018409. URL https://doi.org/10.1063/1.5018409. [4] Steven L. Brunton, Joshua L. Proctor, and J. Nathan Kutz. Discovering governing equations from data by sparse identification of nonlinear dynamical systems. Proceedings of the National Academy of Sciences, 113(15):3932–3937, April 2016. ISSN 0027-8424, 1091-6490. doi: 10.1073/pnas.1517384113. URL https://pnas.org/doi/full/10.1073/pnas.1517384113. [5] Jochen Bröcker, Simon Driscoll, Tobias Necker, José Rodríguez, Helen Dacre, Natalie Harvey, and Zied Ben Bouallègue. Verification of AI–based environmental forecasting systems: What can we do, what do we need to do, and what are the challenges? Journal of the European Meteorological Society, 4:100032, June 2026. ISSN 29506301. doi: 10.1016/j.jemets.2026.100032. URL https://linkinghub.elsevier. com/retrieve/pii/S2950630126000013. [6] Andrew J. Charlton-Perez, Helen F. Dacre, Simon Driscoll, Suzanne L. Gray, Ben Harvey, Natalie J. Harvey, Kieran M. R. Hunt, Robert W. Lee, Ranjini Swaminathan, Remy Vandaele, and Ambrogio Volonté. Do AI models produce better weather forecasts than physics-based models? A quantitative evaluation case study of Storm Ciarán. npj Climate and Atmospheric Science, 7(1):93, April 2024. ISSN 2397-3722. doi: 10.1038/s41612-024-00638-w. URL https://www.nature.com/articles/s41612-024-00638-w. [7] Mickaël D. Chekroun, Alexis Tantet, Henk A. Dijkstra, and J. David Neelin. Ruelle–Pollicott Resonances of Stochastic Systems in Reduced State Space. Part I: Theory. Journal of Statistical Physics, 179(56):1366–1402, June 2020. ISSN 0022-4715, 1572-9613. doi: 10.1007/s10955-020-02535-x. URL https://link.springer.com/10.1007/s10955-020-02535-x. [8] Ricky T. Q. Chen, Yulia Rubanova, Jesse Bettencourt, and David K Duvenaud. Neural Ordinary Differential Equations. In Advances in Neural Information Processing Systems, volume 31. Curran Associates, Inc., 2018. URL https://proceedings.neurips.cc/paper_files/paper/2018/hash/ 69386f6bb1dfed68692a24c8686939b9-Abstract.html. [9] Matthew J. Colbrook, Lorna J. Ayton, and Máté Szőke. Residual dynamic mode decomposition: robust and verified Koopmanism. Journal of Fluid Mechanics, 955:A21, January 2023. ISSN 00221120, 1469-7645. doi: 10.1017/jfm.2022.1052. URL https://www.cambridge.org/core/product/ identifier/S0022112022010527/type/journal_article. [10] Fabrizio Falasca. Probing forced responses and causality in data-driven climate emulators: Conceptual limitations and the role of reduced-order models. Physical Review Research, 7(4):043314, December 2025. doi: 10.1103/2f4r-k8lr. URL https://link.aps.org/doi/10.1103/2f4r-k8lr. [11] Michael Ghil and Valerio Lucarini. The physics of climate variability and climate change. Reviews of Modern Physics, 92(3):035002, July 2020. doi: 10.1103/RevModPhys.92.035002. URL https: //link.aps.org/doi/10.1103/RevModPhys.92.035002. [12] Gregory J. Hakim and Sanjit Masanam. Dynamical Tests of a Deep Learning Weather Prediction Model. Artificial Intelligence for the Earth Systems, 3(3):e230090, July 2024. ISSN 2769-7525. doi: 10.1175/AIES-D-23-0090.1. URL https://journals.ametsoc.org/view/journals/aies/3/3/ AIES-D-23-0090.1.xml. [13] Herbert Jaeger. The" echo state" approach to analysing and training recurrent neural networks-with an erratum note’. Bonn, Germany: German National Research Center for Information Technology GMD Technical Report, 148, January 2001.

5

[14] Patrick Kidger, James Foster, Xuechen Li, and Terry J. Lyons. Neural SDEs as Infinite-Dimensional GANs. In Proceedings of the 38th International Conference on Machine Learning, pages 5453–5463. PMLR, July 2021. URL https://proceedings.mlr.press/v139/kidger21b.html. [15] Diederik P. Kingma and Jimmy Ba. Adam: A Method for Stochastic Optimization, January 2017. URL http://arxiv.org/abs/1412.6980. arXiv:1412.6980 [cs.LG]. [16] Milan Korda and Igor Mezić. On Convergence of Extended Dynamic Mode Decomposition to the Koopman Operator. Journal of Nonlinear Science, 28(2):687–710, April 2018. ISSN 1432-1467. doi: 10.1007/s00332-017-9423-0. URL https://doi.org/10.1007/s00332-017-9423-0. [17] Remi Lam, Alvaro Sanchez-Gonzalez, Matthew Willson, Peter Wirnsberger, Meire Fortunato, Ferran Alet, Suman Ravuri, Timo Ewalds, Zach Eaton-Rosen, Weihua Hu, Alexander Merose, Stephan Hoyer, George Holland, Oriol Vinyals, Jacklynn Stott, Alexander Pritzel, Shakir Mohamed, and Peter Battaglia. Learning skillful medium-range global weather forecasting. Science, 382(6677):1416–1421, December 2023. ISSN 0036-8075, 1095-9203. doi: 10.1126/science.adi2336. URL https://www.science.org/ doi/10.1126/science.adi2336. [18] Valerio Lucarini. Interpretable and equation-free response theory for complex systems. Philosophical Transactions of the Royal Society A: Mathematical, Physical and Engineering Sciences, 384(2322): 20250081, June 2026. ISSN 1364-503X. doi: 10.1098/rsta.2025.0081. URL https://doi.org/10. 1098/rsta.2025.0081. [19] Valerio Lucarini and Mickaël D. Chekroun. Theoretical tools for understanding the climate crisis from Hasselmann’s programme and beyond. Nature Reviews Physics, 5(12):744–765, November 2023. ISSN 2522-5820. doi: 10.1038/s42254-023-00650-8. URL https://www.nature.com/articles/ s42254-023-00650-8. [20] Valerio Lucarini, Manuel Santos Gutiérrez, John Moroney, and Niccolò Zagli. A general framework for linking free and forced fluctuations via Koopmanism. Chaos, Solitons & Fractals, 202:117540, January 2026. ISSN 09600779. doi: 10.1016/j.chaos.2025.117540. URL https://linkinghub.elsevier. com/retrieve/pii/S096007792501553X. [21] Umberto Marini Bettolo Marconi, Andrea Puglisi, Lamberto Rondoni, and Angelo Vulpiani. Fluctuation–dissipation: Response theory in statistical physics. Physics Reports, 461(4):111–195, June 2008. ISSN 0370-1573. doi: 10.1016/j.physrep.2008.02.002. URL https://www.sciencedirect.com/science/ article/pii/S0370157308000768. [22] Jaideep Pathak, Brian Hunt, Michelle Girvan, Zhixin Lu, and Edward Ott. Model-Free Prediction of Large Spatiotemporally Chaotic Systems from Data: A Reservoir Computing Approach. Physical Review Letters, 120(2):024102, January 2018. doi: 10.1103/PhysRevLett.120.024102. URL https: //link.aps.org/doi/10.1103/PhysRevLett.120.024102. [23] Jaideep Pathak, Shashank Subramanian, Peter Harrington, Sanjeev Raja, Ashesh Chattopadhyay, Morteza Mardani, Thorsten Kurth, David Hall, Zongyi Li, Kamyar Azizzadenesheli, Pedram Hassanzadeh, Karthik Kashinath, and Animashree Anandkumar. FourCastNet: A Global Data-driven High-resolution Weather Model using Adaptive Fourier Neural Operators, February 2022. URL http://arxiv.org/abs/2202. 11214. arXiv:2202.11214 [physics.ao-ph]. [24] Grigorios A. Pavliotis. Stochastic Processes and Applications: Diffusion Processes, the Fokker-Planck and Langevin Equations, volume 60 of Texts in Applied Mathematics. Springer, New York, NY, 2014. ISBN 978-1-4939-1322-0 978-1-4939-1323-7. doi: 10.1007/978-1-4939-1323-7. URL https://link. springer.com/10.1007/978-1-4939-1323-7. [25] Alec Radford, Jeffrey Wu, Rewon Child, David Luan, Dario Amodei, and Ilya Sutskever. Language Models are Unsupervised Multitask Learners. Technical report, OpenAI, 2019. URL https://cdn.openai.com/ better-language-models/language_models_are_unsupervised_multitask_learners.pdf. [26] Stephan Rasp, Stephan Hoyer, Alexander Merose, Ian Langmore, Peter Battaglia, Tyler Russell, Alvaro Sanchez-Gonzalez, Vivian Yang, Rob Carver, Shreya Agrawal, Matthew Chantry, Zied Ben Bouallegue, Peter Dueben, Carla Bromberg, Jared Sisk, Luke Barrington, Aaron Bell, and Fei Sha. WeatherBench 2: A Benchmark for the Next Generation of Data-Driven Global Weather Models. Journal of Advances in Modeling Earth Systems, 16(6):e2023MS004019, 2024. ISSN 1942-2466. doi: 10.1029/2023MS004019. URL https://onlinelibrary.wiley.com/doi/abs/10.1029/2023MS004019. [27] David Ruelle. A review of linear response theory for general differentiable dynamical systems. Nonlinearity, 22(4):855, March 2009. ISSN 0951-7715. doi: 10.1088/0951-7715/22/4/009. URL https://doi.org/ 10.1088/0951-7715/22/4/009.

6

[28] Ashish Vaswani, Noam Shazeer, Niki Parmar, Jakob Uszkoreit, Llion Jones, Aidan N Gomez, Ł ukasz Kaiser, and Illia Polosukhin. Attention is All you Need. In Advances in Neural Information Processing Systems, volume 30. Curran Associates, Inc., 2017. URL https://proceedings.neurips.cc/paper_ files/paper/2017/hash/3f5ee243547dee91fbd053c1c4a845aa-Abstract.html. [29] Matthew O. Williams, Ioannis G. Kevrekidis, and Clarence W. Rowley. A Data–Driven Approximation of the Koopman Operator: Extending Dynamic Mode Decomposition. Journal of Nonlinear Science, 25 (6):1307–1346, December 2015. ISSN 0938-8974, 1432-1467. doi: 10.1007/s00332-015-9258-5. URL http://link.springer.com/10.1007/s00332-015-9258-5. [30] Niccolò Zagli, Matthew J. Colbrook, Valerio Lucarini, Igor Mezić, and John Moroney. Bridging the Gap Between Koopmanism and Response Theory: Using Natural Variability to Predict Forced Response. SIAM Journal on Applied Dynamical Systems, 25(1):196–229, March 2026. ISSN 1536-0040. doi: 10.1137/24M1699206. URL https://epubs.siam.org/doi/10.1137/24M1699206. [31] Haoyi Zhou, Shanghang Zhang, Jieqi Peng, Shuai Zhang, Jianxin Li, Hui Xiong, and Wancai Zhang. Informer: Beyond Efficient Transformer for Long Sequence Time-Series Forecasting. Proceedings of the AAAI Conference on Artificial Intelligence, 35(12):11106–11115, May 2021. ISSN 2374-3468. doi: 10.1609/aaai.v35i12.17325. URL https://ojs.aaai.org/index.php/AAAI/article/view/17325.

A

Operator theory

This appendix collects the operator-theoretic objects that the response diagnostic in the main text is built on, following the Koopmanism Response framework [30]. Standard references for the underlying stochastic-analysis facts are [24] for the Kolmogorov operator and [7] for stochastic Ruelle-Pollicott resonances.

A.1

Stochastic dynamics and the invariant measure

We consider an autonomous Itô SDE at statistical stationarity, x ∈ M ⊆ Rd ,

dx = F(x) dt + ς dWt ,

(2)

with (generally nonlinear) drift F, isotropic additive noise of intensity ς > 0, and a d-dimensional Wiener process Wt . The one-time statistics are carried by the probability density ρ(x, t), which solves the Fokker-Planck (forward Kolmogorov) equation 2

∂t ρ = L0 ρ := −∇· (Fρ) + ς2 ∇2 ρ.

(3)

Under hypoellipticity of L0 and mild growth conditions on F, the process (2) admits a unique, ergodic invariant measure dµ0 (x) = ρ0 (x) dx that is absolutely R continuous with respect to Lebesgue and has a smooth density [24]. Stationary expectations are ⟨f ⟩µ0 = M f dµ0 , and the two pillars of standard surrogate validation (attractor geometry, long-run marginals) are exactly statements about µ0 .

A.2

The Markov semigroup and its generator

The statistics of observables evolve under the Markov semigroup P = (Pt )t≥0 , Z (Pt f )(x) = f (y) p(y, t | x) dy = E[ f (x(t)) | x(0) = x ] ,

(4)

M

where p(y, t | x) is the transition density. On L2µ0 = {f : semigroup, so it can be written in exponential form

R

|f |2 dµ0 < ∞} this is a strongly continuous

Pt = etK0 ,

(5)

with (generally unbounded) generator K0 , the Kolmogorov (backward) operator. It is the second-order differential operator 2

K0 = F(x) · ∇ + ς2 ∇2 ,

(6)

and it advances expectations of observables, ∂t f = K0 f with f (x, t) = (Pt f )(x). In the deterministic limit ς → 0, K0 = F · ∇ reduces to the usual Koopman generator. We retain the operator name “Koopman generator” for K0 in the body, following the convention of [30].

7

A.3

Why the noise: deterministic vs. stochastic spectra

The regularization ς > 0 in (2) is operator-theoretic, not cosmetic, and it is what makes the response machinery well posed. The deterministic Lorenz-63 flow carries a singular SRB measure, and its Koopman operator is unitary on L2µ0 : the dynamics has already settled onto the attractor, so no decay toward it can be observed and the spectrum sits on the imaginary axis (equivalently, on the unit circle for the time-τ map). Adding non-degenerate noise makes L0 hypoelliptic, yielding a unique smooth ρ0 and a quasi-compact semigroup, so that a discrete set of relaxation modes emerges and the data-driven response reconstruction becomes well posed [7, 30].

A.4

Ruelle-Pollicott resonances and the spectral decomposition

Assume P is quasi-compact, i.e. the eigenvalues of K0 with real part above the essential growth bound ress < 0 are discrete with finite multiplicity. By the spectral mapping theorem, if λj is an eigenvalue of K0 with eigenfunction φj , then eλj t is an eigenvalue of Pt with the same eigenfunction. These discrete eigenvalues {λj }j≥0 , ordered as λ0 = 0 > Re λ1 ≥ Re λ2 ≥ · · · , are the (stochastic) Ruelle-Pollicott resonances, and they give a physically meaningful decomposition of the semigroup, X λ t etK0 = e j Πj + R(t), (7) j≥0

where Πj projects onto the eigenspace of λj and R(t) collects the essential-spectrum contribution, whose norm decays and becomes negligible as t → ∞ under quasi-compactness. The resonances are intrinsic relaxation rates: a nonzero imaginary part produces oscillatory correlations, the signature of nonequilibrium dynamics, while the spectral gap γ = − Re λ1 controls the slowest decay. An accumulation of resonances toward λ0 = 0 signals diverging sensitivity (critical slowing down). This is the structure a learned surrogate must reproduce, and it is precisely what neither forecast skill nor invariant-measure agreement constrains. For zero-mean observables f, g ∈ L2µ0 , the time-lagged correlation function inherits (7), Z X
Cf,g (t) = etK0 f (x) g ∗ (x) dµ0 (x) = αj (f, g) eλj t , M

(8)

j≥1

a sum of exponentially decaying (and possibly oscillating) modes; the rates λj are intrinsic to the system, whereas the coefficients αj (f, g) encode the projection of the chosen observables onto the Kolmogorov modes.

A.5

Extended Dynamic Mode Decomposition (EDMD)

We estimate the spectral objects in (7)-(8) from a single long unperturbed trajectory via EDMD [29]. Given snapshot pairs D = {(xi , yi = xi+1 )}M i=1 sampled at step τ and a dictionary Ψ(x) = (ψ1 (x), . . . , ψN (x)) spanning FN ⊆ L2µ0 , form the Gram and one-step (action) matrices X ∗ X ∗ 1 1 G= M Ψ (xi )Ψ(xi ), A= M Ψ (xi )Ψ(yi ), K = G+ A, (9) i

i

so that Gab = ⟨ψa , ψb ⟩µ0 and Aab = ⟨ψa , K0 ψb ⟩µ0 (here G+ is a pseudoinverse). The matrix K is the Galerkin approximation of the semigroup on FN . Writing its eigenpairs as Kξk = νk ξk , the Kolmogorov eigenvalues and eigenfunctions follow from the spectral mapping theorem, λk =

ln νk , τ

φk (x) = Ψ(x) ξk ,

(10)

and the Koopman-mode coefficients of any observable f = Ψh ∈ FN are f = W′ h, with W the matrix of left P eigenvectors of K. Correlations then reconstruct as Cf,g (t) = k f˜k eλk t , where the projection coefficients are read off from the EDMD Gram matrix, ⟨φl , φk ⟩µ0 = (Ξ′ GΞ)lk with Ξ = [ξ1 · · · ξN ]. As N → ∞ with {ψi } a basis of L2µ0 , EDMD converges to K0 in the strong operator topology [16], which justifies the reconstruction of correlation and response operators from data. In the experiments we take Ψ to be a three-dimensional Fourier dictionary of size (2Kmax + 1)3 on a bounded rescaling of the attractor domain, and regularize G+ by a truncated SVD (TSVD).

A.6

Certifying the spectrum: ResDMD residuals

Galerkin truncation can produce spurious (spectrally polluted) modes, so a finite EDMD spectrum must be verified rather than trusted. ResDMD supplies a data-driven residual for each candidate eigenpair using one additional matrix computed from the same snapshots [9], X ∗ 1 (11) L= M Ψ (yi )Ψ(yi ), Lab = ⟨K0 ψa , K0 ψb ⟩µ0 . i

8

Table 2: Per-mode eigenfunction overlaps. Each cell: median overlap O (Eq. (13)) over scored rollouts (top) and the percentage of rollouts that pass the mode (bottom); shading: redder = more rollouts fail. Oscillatory passes use the applied threshold clip(mean − 3 std, 0.5, 0.99) from the 16 calibration truths (first row; uncapped value in parentheses); slow modes use the uncapped threshold. †: more than half of the passes are decided by the 0.99 cap, i.e. would fail under the uncapped threshold. Merged: the joint subspace of the two near-degenerate pairs (9.85 and 9.98), with null 0.9982 ± 0.0003 over the 16 calibration truths; reported here only. Bottom block: SINDy trained on 100 time units (App. C.1) and the true drift rolled out on the NODE/NSDE integrator paths (App. B.8). n: scored rollouts. Slow (Re λ) Model

n

-0.366

-0.508

Oscillatory (ω) 8.42

9.85

9.98

16.68

merged

0.9450 0.9216 0.9900 0.9681 0.9601 0.9900 0.9900 (0.9450) (0.9216) (0.9984) (0.9681) (0.9601) (0.9960) (0.9974)

threshold

Truth

50

0.992 96%

0.991 94%

0.999 100%

0.990 98%

0.991 100%

0.998 100%

0.998 100%

SINDy

10

MLP-drift

80

MLP-flow

80

RC

78

NODE

80

NSDE

80

0.994 100% 0.043 0% 0.820 8% 0.923 27% 0.976 86% 0.980 89%

0.986 100% 0.319 0% 0.754 5% 0.887 29% 0.962 78% 0.960 78%

0.999 100% 0.607 0% 0.963 0% 0.998 99%† 0.998 100%† 0.998 100%†

0.990 90% 0.373 0% 0.751 0% 0.895 0% 0.952 30% 0.953 29%

0.993 90% 0.290 0% 0.691 0% 0.856 0% 0.935 21% 0.929 24%

0.998 100% 0.115 0% 0.245 0% 0.997 97% 0.996 100%† 0.996 100%†

0.998 100% 0.231 0% 0.486 0% 0.996 97%† 0.996 82%† 0.996 90%†

Transformer (expl.) 80

0.981 84%

0.974 86%

0.997 100%†

0.949 29%

0.954 42%

0.996 100%†

0.996 85%†

SINDy-100

10

G0 , NODE path

10

G0 , NSDE path

10

0.977 100% 0.988 100% 0.990 100%

0.969 100% 0.975 90% 0.978 100%

0.995 100%† 0.994 100%† 0.993 100%†

0.904 10% 0.978 70% 0.982 60%

0.748 0% 0.967 70% 0.971 70%

0.995 100%† 0.993 100%† 0.993 100%†

0.993 100%† 0.988 10%† 0.744 30%†

For a candidate pair (λ, ξ) with g = Ψξ, the relative residual is res(λ, ξ) =

(K0 − λ) g µ

0

∥g∥µ0

=

!1/2 ξ∗ L − λA∗ − λ̄A + |λ|2 G ξ , ξ∗ Gξ

(12)

which is exact in the L2µ0 norm and requires no knowledge of the true spectrum. Residuals of the leading candidate pairs are stored for every scored rollout and used as a diagnostic, not as an acceptance gate. In the main text this is what lets us treat the reference resonances as ground truth and, symmetrically, what makes the mode-resolved surrogate test meaningful. We compare a surrogate’s leading eigenfunctions against the reference not by eigenvalue proximity alone but by the µ0 -weighted subspace overlap O(φ, φ̃) =

|⟨φ, φ̃⟩µ0 | , ∥φ∥µ0 ∥φ̃∥µ0

(13)

generalized to the principal angles between eigenspaces for near-degenerate clusters (e.g. the near-degenerate 9.85/9.98 pair). This overlap is the discriminating quantity in Table 1: it is sensitive to eigenfunction structure, which the eigenvalue-only self-consistency check Cφk (t) ∝ eλk t and the invariant statistics both miss. RC, NODE and NSDE miss the individual eigenvectors of the near-degenerate oscillatory pair in most rollouts while recovering their joint subspace; the MLPs recover neither, and only SINDy resolves the split (Table 2; Fig. 2a). Table 2 gives every model’s per-mode overlaps; the thresholds of the ω ≈ 8.42 and 16.68 pairs sit at the 0.99 cap (uncapped 0.9984 and 0.9960), so passes on these modes can be decided by the cap, as for 75 of RC’s 77 passes at 8.42 and 63 and 71 of the 80 and 80 NODE and NSDE passes there, against 0 for SINDy.

9

(b)

oscillatory (direct)

slow Truth

0.96

0.94

1.00

0.98

1.00

1.00

SINDy (full data)

1.00

1.00

1.00

0.90

0.90

1.00

MLP-flow

0.07

0.05

0.00

0.00

0.00

0.00

RC

0.27

0.29

0.99 †

0.00

0.00

0.97

NODE

0.86

0.78

1.00 †

0.30

0.21

1.00 †

NSDE

0.89

0.78

1.00 †

0.29

0.24

1.00 †

real 1

real 2

8.42

9.85

9.98

16.68

Truth

1.05

slow-mode margin m

(a)

passes both

fails invariant only

1.00 0.95 0.90 0.85 0.80 fails modes only

0.20

0.50

fails both

1.00

2.00

5.00

invariant-statistics exceedance e SINDy (full data)

RC

NODE

NSDE

Â(0) outside band

Figure 2: Per-mode passes and the dissociation map. (a) Fraction of rollouts passing each leading mode (slow modes uncapped, oscillatory modes matched directly); †: more than half of the passes are decided by the 0.99 cap. MLP-drift passes no mode in any rollout (not shown). (b) One point per scored rollout: e is the largest exceedance over the four invariant-statistics floors (W1 and ACF as value over floor; switch ratio and f+ as distance from the floor interval’s centre over its half-width), m the smaller of the two slow-mode overlaps over their thresholds; e ≤ 1 is inside every floor and m ≥ 1 passes both slow modes. Hollow markers: χ(0) outside the band. The MLP models are not shown. Table 3: Slow-mode threshold sensitivity. The two slow-mode thresholds rebuilt from the 50 scored fresh truths by a leave-one-out minimum rule: each truth is judged against the minimum overlap of the other truths (its row is the leave-one-out false-alarm rate), every model against the minimum over all truths (λ ≈ −0.366: 0.9257 instead of 0.9450; λ ≈ −0.508: 0.8693 instead of 0.9216). The stored per-mode overlaps are re-thresholded; nothing is re-run, and Table 8 keeps the calibration-set thresholds. Rollouts inside all four invariant-statistics floors and the χ(0) band that fail the slow modes under the rebuilt thresholds: RC 14 of 26, truth 1 of 48 (Fisher exact p = 2.6 × 10−7 ). slow modes failed (%) n

pass, rebuilt thresholds (%)

calibration set rebuilt λ ≈ −0.366 λ ≈ −0.508

Truth (perfect surrogate) SINDy (full data) MLP-drift MLP-flow RC NODE NSDE SINDy-100 G0 (NODE path) G0 (NSDE path)

50 10 80 80 78 80 80 10 10 10

6 0 100 97 76 24 23 0 10 0

2 0 100 94 55 7 9 0 0 0

98 100 0 10 45 95 95 100 100 100

98 100 0 19 62 94 91 100 100 100

Transformer (expl.)

80

19

9

91

97

Rebuilding the two slow-mode thresholds from the 50 fresh truths (leave-one-out minimum: 0.9257 and 0.8693 instead of 0.9450 and 0.9216) lowers RC’s slow-mode failure rate from 76% to 55% and NODE’s from 24% to 7%, and the RC exhibit becomes 14 of 26 rollouts against 1 of 48 truths (Fisher p = 2.6 × 10−7 ; Table 3).

B

Response theory

We now perturb the stationary dynamics and recall how linear response theory fixes the forced behaviour entirely from unperturbed statistics, through a generalized fluctuation-dissipation relation and the same resonances {λj } of Appendix A. This is the content that turns “did the surrogate learn the physics?” into a falsifiable spectral question. The formulation follows [30].

10

B.1

Perturbation and the perturbed Fokker-Planck equation

Perturb the drift by a small, time-modulated phase-space pattern, F → F + ε X(x) T (t), with ε ≪ 1, time profile T (t), and spatial profile X. The density then obeys

∂t ρ = L0 + ε T (t) Lp ρ, ρ(x, 0) = ρ0 (x), Lp · = −∇· X(x) · , (14) where the initial condition encodes that the forcing switches on after the system has reached stationarity. Because the drive is nonautonomous, no invariant measure of (14) exists in general; we work to first order in ε.

B.2

Linear response and the generalized FDT

To first order the expectation of any observable f obeys
⟨f ⟩µ (t) = ⟨f ⟩µ0 + ε Gf ⋆ T (t) + o(ε), with the causal Green’s function given by the generalized fluctuation-dissipation theorem Z Lp ρ0 tK0
L p ρ0 Gf (t) = Θ(t) Γ(x) = e f dµ0 = Θ(t) Cf,Γ (t), , ρ ρ0 0 M

(15)

(16)

where Θ is the Heaviside step. Equation (16) is the crux: the forced response is a specific unperturbed correlation function between the observable f and the conjugate observable Γ. Away from equilibrium there is no firstprinciples form for ρ0 , so Γ is intrinsic and not known a priori; nonetheless (15)-(16) let us estimate the response purely from unperturbed trajectory data.

B.3

Mode-resolved response and the static susceptibility

Substituting the spectral decomposition (7) into (16) gives the response as a sum over the same resonances, X Gf (t) = Θ(t) βj (f, Γ) eλj t , (17) j≥1

with coefficients βj (f, Γ) that project both f and Γ onto the Kolmogorov eigenfunctions. This is a genuinely interpretable diagnostic: rather than a scalar error, the response decomposes into contributions from distinct modes of natural variability, so a surrogate can be interrogated mode by mode. We summarize the response of the tracked observable (f = x for Lorenz-63) by its zero-frequency limit, the static susceptibility, Z ∞ δ⟨x⟩∞ χ(0) = Gx (t) dt = , (18) ε 0 i.e. per unit R ∞ λ shift P the steady-state P forcing under a switched-on profile T (t) = Θ(t). From (17), χ(0) = k t dt = − β (x, Γ) e k k k βk (x, Γ)/λk , convergent because Re λk < 0. Evaluated from the sampled 0 EDMD spectrum (νk = eλk τ ), the tail integral becomes a geometric series over the sampling grid t = nτ , X X X n X α ∗ M γk k χ(0) ≈ τ Gx (nτ ) = τ βk νk = τ , (19) 1 − νk n≥0

k

n≥0

k

which is the estimator used in the body. Here αk∗ and γk are the EDMD projection coefficients of the observable x and of the conjugate observable Γ onto the k-th Kolmogorov eigenfunction, and M denotes the eigenfunction overlap (Gram) matrix ⟨φl , φk ⟩ contracting the two, not the snapshot count M of Appendix A. Because χ(0) is a mode sum weighted by 1/(1 − νk ), it inherits nontrivial variance at finite trajectory length, which is why we calibrate a ground-truth band by re-evaluating the instrument on independent noise realizations of the true system.

B.4

Instrument diagnostics and finite-length variance

Several checks quantify how far the finite-dictionary, finite-trajectory instrument departs from its exact counterpart, and they are worth reporting because they bound what the main-text comparisons can resolve. With Γ = −∂x ln ρ0 and forcing along x, the value of the Green’s function at zero lag is fixed analytically: from (16), Z Z Z
Gx (0) = x − ∂x ln ρ0 dµ0 = − x ∂x ρ0 dx = ρ0 ∂x x dx = 1, (20) M

M

M

the boundary term vanishing for the smooth, decaying ρ0 . Nothing in the pipeline normalizes Gx , so the measured deviation from unity is a free, per-series estimate of dictionary-truncation error: |Gx (0) − 1| is at most 0.83% over the truths and at most 1.92% over all scored rollouts.

11

Table 4: Where χ(0) sits. Signed share of χ(0) carried by the leading resonances (median [min, max]), their share of the summed absolute mode contributions, and the fraction of χ(0) accumulated by t = 1. n reference truth SINDy-100 RC NODE NSDE

leading |share| accumulated by t = 1

leading share

1 0.34 10 0.33 [-0.15, 0.60] 10 -0.06 [-0.52, 0.40] 10 0.29 [-0.06, 0.60] 10 0.24 [-0.36, 0.50] 10 0.08 [-0.41, 0.57]

– 0.06 0.04 0.07 0.04 0.03

0.65 0.89 [0.60, 1.32] 0.97 [0.76, 1.30] 0.87 [0.65, 1.65] 0.87 [0.62, 1.54] 0.87 [0.62, 1.51]

Table 5: Null calibration. Leave-one-out out-rates on the independent true ensemble (floors, E1b band) and the pooled band; out-of-range counts of the 50 runner truths and the one-sided binomial test against the leave-one-out rate. The E1b band test fails (p < 0.05), so the band is the pooled range [0.162, 0.309] over 100 truths; no floor test fails, so all four floors stay as built. Runner and E1b χ(0) variances: F = 1.46, p = 0.191. W1 ACF switch leave-one-out out-rate 0.02 0.02 runner truths outside 0 0 one-sided binomial p 1 1

B.4.1

0.04 2 0.6

f+

all four band (E1b) band (pooled)

0.04 0 1

0.10 – –

0.04 6 0.0144

0.02 – –

χ0 diagnostics

P χ(0) is evaluated in closed form through (19), not by truncating τ n Gx (nτ ) at a finite lag, and MLP-flow accrues slowly because its slowest real mode decays too slowly (−0.146 against the true −0.389): only 33% of its χ(0) has accrued by t = 1, against a median 89% for the truths. The invariant-statistics floors are the maxima of each distance over the 50 independent true realizations scored against the held-out reference realization (seed 101) rather than the master trajectory, which is every model’s training data; no floor fails its leave-one-out null test (Table 5). The χ(0) band is the range of the 100 pooled true realizations, because the range of the independent-IC ensemble alone left 6 of the 50 truths scored through the runner outside (one-sided binomial P = 0.0144 against its leave-one-out rate); Table 6 gives the in-band counts under each alternative band. The slowest real eigenvalue varies little across the 50 truths (coefficient of variation 0.0508), while the share of χ(0) carried P by the leading resonances varies strongly (coefficient of variation 1.09 over 10 truths). Since χ(0) = − k βk /λk and the λk barely move, the finite-length variance enters mainly through the projection amplitudes βk , hence through the Method-3 estimate of Γ (App. B.5), which is a time average along a single trajectory.

B.5

Data-driven estimation of Γ: Method 3 (integration by parts)

The one nontrivial ingredient is the conjugate observable Γ = Lp ρ0 /ρ0 . [30] give three routes to Gf . Methods 1 and 2 build Γ by first estimating ρ0 (e.g. by binning) and then dividing by it; this is fragile whenever µ0 concentrates near a (near-)fractal attractor, exactly the Lorenz-63 situation. Method 3 avoids density estimation altogether. It seeks the µ0 -weighted least-squares projection of Γ onto the dictionary, Z Z


1 T c = G+ ∆, ∆i = ⟨ψi , Γ⟩µ0 = X·∇ψi∗ ρ0 dx ≈ X x(t) ·∇ψi∗ x(t) dt, (21) T 0 M where the middle equality is an integration by parts that moves the derivative off ρ0 (boundary terms vanish for the smooth, decaying ρ0 ), and the last step replaces the µ0 -average by a time average along the same unperturbed trajectory used for EDMD. The moments ∆i therefore require no estimate of ρ0 and no division by it; if the dictionary is closed under differentiation (as for Fourier), they use no information beyond EDMD itself, and the method is markedly data-efficient. The coefficients of Γ in the eigenfunction basis follow from c via the left eigenvectors, feeding directly into (17).

B.6

Lorenz-63 specialization

We force the x-equation only, X(x) = (1, 0, 0)⊤ , with T (t) = Θ(t), and track f = x. The conjugate observable then simplifies to Γ(x) = − ∂x ln ρ0 (x), (22)

12

Table 6: χ(0) in-band counts under each band (sensitivity rows), and diverged rollouts (non-finite or a coordinate beyond 1.5× its true maximum: |x|, |y|, |z| > 31.7, 43.6, 76.1), excluded from all rates. Model

n

Truth (perfect surrogate) SINDy (full data) MLP-drift MLP-flow RC NODE NSDE SINDy-100 G0 (NODE path) G0 (NSDE path)

50 10 80 80 78 80 80 10 10 10

pooled [0.162, 0.309] E1b [0.163, 0.266] N = 18 [0.18, 0.324] runner [0.162, 0.309] diverged 50 10 7 1 72 73 73 7 10 7

44 10 5 0 67 61 62 7 5 3

44 9 8 2 59 66 67 5 10 8

50 10 7 1 72 73 73 7 10 7

0 0 0 0 2 0 0 0 0 0

Table 7: Discretization sensitivity (G0): the true drift through the NODE/NSDE RK4 step h, failure rates per rung and the shift of the median χ(0) from the truth median. n Truth 50 G0, NODE path 10 G0, NSDE path 10

invariant (%) χ(0) (%) slow (%) oscill. (%) χ(0) median 4 100 100

0 0 30

6 10 0

2 40 40

0.221 0.265 0.275

shift – 0.0443 0.0549

RT and the Method 3 moments (21) reduce to ∆i ≈ T1 0 ∂x ψi∗ (x(t)) dt. We force x because it is the most strongly coupled variable, driving both y and z through the nonlinear terms, so a single-component forcing maximally probes the cross-variable response and is the cleanest first stress test of the framework. The reference response is the instrument reconstruction on the held-out true realization (seed 101).

B.7

Why the response is a stringent diagnostic

Equations (16)-(17) show that the entire forced response is determined by unperturbed objects: the resonances {λj }, their eigenfunctions φj , and the projection of the intrinsic conjugate observable Γ. None of these is constrained by a short-horizon forecast loss or by agreement of the invariant measure. A surrogate trained only on unperturbed data can therefore match the attractor, the wing-switching statistics, and even the scalar χ(0), while still carrying the wrong eigenfunction structure, a discrepancy that only the mode-resolved test of Appendix A.6, applied through this response machinery, resolves. Response theory thus supplies a rigorous, physically grounded criterion for whether a learned simulator respects dynamical structure rather than merely fitting the attractor.

B.8

Discretization control

To separate learning error from the error of the NODE/NSDE integrator, we roll out the true drift through the same RK4 path at h = ∆tfit = 0.01 with oracle noise, on both code paths (NODE: constant ς; NSDE: the diffusion network replaced by the constant ς), 10 draws each (Table 7). The switch ratio falls below one in every draw (median 0.953 and 0.948, maximum 0.975 and 0.968), whereas NODE and NSDE lie above one in 80 of 80 and 80 of 80 rollouts, so their over-switching is learned rather than numerical. The step shifts the median χ(0) by 0.0443 and 0.0549 from the truth median and leaves the slow modes largely intact (slow-mode failure rate at most 10%).

C

Surrogate models

All five surrogates, the exploratory transformer of App. C.5 and the controls of App. B are evaluated through the identical instrument of App. A: same dictionary (Kmax = 4, size 729), same TSVD threshold, same sampling step τ = 0.1, same forcing direction; the surrogates are trained on unperturbed trajectories only. Table 8 gives the per-rollout failure rates of every model on each rung. The transformer is a first step from the formulation ladder toward the architectures of current weather emulators; it was run with a fixed configuration and no tuning, so it is reported separately from the ladder. Table 9 lists the training and rollout settings of every model.

C.1

SINDy

We fit the drift by sequentially thresholded least squares (STLSQ) on a polynomial library of degree 2, with threshold 0.1 on standardized columns. Under Euler-Maruyama the one-step increment of (2) at the raw step

13

Table 8: Per-rollout failure rates, in percent of scored rollouts (redder = fails more often); a rollout fails a rung (invariant statistics, χ(0), slow modes, oscillatory modes) if it fails any check in it. Truth: fresh true realizations scored as a perfect surrogate, whose failure rate is the false-alarm rate; SINDy: trained on the full record. Invariant statistics: any of the four checks, then each check alone, outside its floor (the maximum over the 50 independent true realizations scored against the reference seed 101): W1 ≤ 0.284, ACF distance ≤ 0.0885, switch ratio in [0.985, 1.03], f+ in [0.487, 0.517]. χ(0): mean with the standard deviation in subscript; fails outside [0.162, 0.309], the range of the 100 pooled true realizations (50 independent-IC plus 50 scored through the runner); the truth row is inside by construction, and its leave-one-out out-rate is given in parentheses. Slow modes: the two real reference modes (λ ≈ −0.366 and −0.508), overlap below the uncapped threshold mean − 3 std over 16 calibration truths (0.9450, 0.9216). Oscillatory modes: the pairs at ω ≈ 8.42, 9.85, 9.98, 16.68, matched directly with thresholds clip(mean − 3 std, 0.5, 0.99); the 8.42 and 16.68 thresholds sit at the cap. MLP-drift: 33 of 80 rollouts collapse onto one wing (scored). n counts scored rollouts; diverged rollouts (non-finite, or a coordinate beyond 1.5× its maximum over the true realizations) are left out of the rates: RC 2. Transformer: exploratory, fixed configuration, not tuned (App. C.5). χ(0)

Invariant statistics (%)

Modes (%)

Model

n (div.)

any

W1

ACF

switch

f+

mean(sd)

out (%)

slow

oscill.

Truth

50 (0)

4

0

0

4

0

0.224(0.034)

0 (2)

6

2

SINDy MLP-drift MLP-flow RC NODE NSDE

10 (0) 80 (0) 80 (0) 78 (2) 80 (0) 80 (0)

20 100 100 63 93 100

10 100 100 3 58 60

0 100 93 27 56 69

0 100 100 58 75 79

20 96 100 6 63 61

0.215(0.029) 1.045(1.278) 0.636(0.313) 0.203(0.032) 0.225(0.045) 0.221(0.039)

0 91 99 8 9 9

0 100 98 76 24 23

10 100 100 100 86 83

Transformer (expl.)

80 (0)

79

51

29

53

64

0.200(0.033)

13

19

73

Table 9: Parameters. System, instrument, seeds (top) and per-model training and rollout (bottom). Every non-NSDE surrogate is rolled out with the true noise amplitude s (oracle noise). System Truth scheme Instrument Mode thresholds Seeds

σ, ρ, β = 10, 28, 2.667; noise s = 2; forcing on x Euler–Maruyama, dtraw = 0.001; burn-in 100 t.u.; T = 104 t.u. τ = 0.1, Kmax = 4 (729 functions), TSVD 10−4 , nlead = 6 clip(mean − 3 std, 0.5, 0.99) over 16 calibration truths; grid 20000 points ref 101; grid 909; calibration 1000–1015; floor truths 5000–5049; runner truths 40000, 42001–42049 E2 train 0–7 × rollout 0–9; NODE init 31000–31007; NSDE init 31100–31107; G0 40001, 41001–41009 and 40002, 41011–41019; Transformer init 32000–32007

Model

∆tfit

rollout step

train span (t.u.)

SINDy (full data) SINDy-100 MLP-drift MLP-flow RC NODE NSDE

0.001 0.001 0.01 0.01 0.02 0.01 0.01

0.001 0.001 0.001 0.01 0.02 0.01 0.01

10000 100 10000 10000 2000 10000 10000

Transformer (expl.)

0.01

0.01

10000

rollout noise x += dt F(x) + s sqrt(dt) ξ x += dt F(x) + s sqrt(dt) ξ x += dt F_theta(x) + s sqrt(dt) ξ x = MLP(x) + s sqrt(dt) ξ x = ESN(x, r) +√s sqrt(dt) ξ (oracle s) RK4 step of Fθ + s h ξ√ RK4 step of Fθ + σϕ (x) h ξ (learned) √ xn + ∆θ (xn−K+1 , . . . , xn ) + s h ξ (oracle s)

satisfies E[∆x/∆t | x] = F (x), so regressing increments on state recovers the drift with the noise averaged out [3]; at a coarser step its conditional mean is the secant drift (App. C.2). The library contains the true Lorenz-63 equations, so SINDy is a positive control on the whole pipeline rather than a competitive baseline: if the instrument cannot certify a model whose functional form is exactly right, nothing else it reports is trustworthy. The positive control is SINDy trained on the same 10000 time units of the master as the neural models, at the raw step. As a robustness check on record length we also report a fit on 100 time units (SINDy-100); for transparency, this shorter fit was the control before the master trajectory was refitted, and the full-record fit replaced it after the results were seen. Trained on 100 time units instead, SINDy acquires symmetry-breaking constant terms that skew the lobes (f+ = 0.416); it fails invariant statistics in every rollout yet passes both slow modes in every rollout, the converse of the dissociation in Sec. 3. Table 10 gives both fits’ coefficients and Table 11 the drift RMSE on the grid of 20000 points drawn from the held-out true realization (seed 909); the MLP-drift row is one representative fit (initialization seed 5).

14

Table 10: SINDy coefficients: true, SINDy-100 and SINDy on the full training span (the same recipe; large terms are those above the STLSQ threshold). term

true

SINDy-100

SINDy (full)

σ (y in ẋ) 10.000 10.520 10.001 −(x in ẋ) 10.000 10.755 10.003 ρ (x in ẏ) 28.000 27.734 27.967 β (−z in ż) 2.667 2.402 2.667 const in ẋ 0.000 -3.423 0.000 const in ẏ 0.000 3.139 0.307 const in ż 0.000 -2.169 0.000 active terms 7 30 18 large terms 7 14 8 training span (t.u.) / pairs – 100 / 100000 10000 / 10000000 f+ of the training window – 0.464 0.500

Table 11: Drift RMSE against the true drift on the 20000-point grid drawn from the held-out true realization (seed 909).

SINDy (full data) SINDy-100 NODE (median over 8 fits) MLP-drift secant bias of the MLP-drift target (exact) secant bias, leading order (∆t/2)JF F

C.2

x

y

z

0.014 0.717 0.577 2.048 2.192 2.439

0.070 0.537 1.125 4.399 4.264 4.740

0.043 0.647 0.961 4.213 4.541 5.056

MLP

The MLP is a map R3 → R3 with 2 hidden tanh layers of width 64 and standardized inputs and outputs, trained with Adam (learning rate 10−3 , batch 2048) [15] on the mean-squared error, with early stopping on a random 10% validation split (patience 5, at most 25 epochs). It is trained on 10000 time units of the master at ∆tfit = 0.01, either on the increment ∆x/∆t (MLP-drift) or on the next state (MLP-flow); the two share data, architecture and seed, so the target is the only difference. A pre-refit sweep over width and depth did not reduce the drift error. At ∆tfit the conditional mean of ∆x/∆t is not F but the secant drift F + (∆t/2) JF F + O(∆t2 ) (the Itô 2 term ς2 ∇2 F vanishes because each component of the Lorenz-63 drift is at most bilinear), and the drift MLP learns it: its error against F matches the secant bias of its target to within 7.21% per component (Table 11). Integrated with Euler–Maruyama at δt = 0.001, MLP-drift collapses onto one wing in 33 of 80 rollouts. The flow map restores two-wing switching (0 of 80 rollouts one-wing; median switch ratio 0.725) but not the measure (invariant-statistics failure rate 100%, median f+ = 0.474); its slowest real mode decays too slowly (−0.146 against the true −0.389), so only 33% of its χ(0) has accrued by t = 1 (App. B.4).

C.3

Reservoir Computing

The reservoir is a leaky echo-state network [13] with N = 1000 units, spectral radius 1, leak rate 0.3, input scaling 0.5, washout 1000 steps and ridge readout regularization 10−4 , trained on 2000 time units of the master at ∆tfit = 0.02 and advanced as a native flow map with injected noise. Like every other surrogate, it is trained on the Euler–Maruyama master, i.e. on noisy data from the stochastic system; its hyperparameters were not retuned for it. The configuration was selected on valid-prediction horizon against a reference trajectory of the deterministic Lorenz-63 system, subject to an attractor gate (bounded, bimodal, plausible wing-switch ratio); response metrics were never used for selection. Two implementation faults had to be fixed first, both of which silently degrade the attractor rather than failing loudly: an off-by-one in the closed-loop rollout seed, and a washout discontinuity in which warm-up ended at the last training state while the rollout restarted from the first. Current RC results are in Tables 1, 8 and 2.

C.4

Neural ODE and neural SDE

The NODE uses the MLP network of App. C.2 as a drift Fθ , trained through a differentiable RK4 step at h = ∆tfit = 0.01 on the same data with the same validation fraction: from each window start the model is unrolled for L steps, the mean-squared error against the master is taken in standardized coordinates and minimized with Adam, with early stopping on the held-out multi-step validation MSE. We swept L ∈ {5, 20} and learning rates {10−3 , 3 × 10−4 } and selected by validation MSE at a common horizon L = 20 (L = 5, learning rate 10−3 ; Table 12); the NODE is rolled out with the same RK4 step and oracle noise. The NSDE

15

Table 12: NODE sweep (one seed): validation MSE (×10−3 , standardized) at horizons L = 1, 5, 20; selection by the lowest MSE at the common horizon L = 20 (first row). Drift RMSE on the grid is reported only. L

MSE1 MSE5 MSE20 epochs train (s) drift RMSE x/y/z

lr

5 0.001 5 0.0003 20 0.001 20 0.0003

0.547 0.548 0.549 0.548

1.543 1.547 1.555 1.551

6.225 6.242 6.268 6.239

25 25 14 25

177 177 342 601

0.490/0.921/0.877 0.648/1.078/0.986 0.841/1.483/1.479 0.800/1.229/1.260

Table 13: NSDE learned noise relative to the true s, at h = 0.01; a = (σ, 1, β).

NSDE mean σϕ /s (median of 80) NSDE constant σ̂/s true p drift, step h, constant σ̂/s (1 − e−2ah )/(2ah)

x

y

z

0.9587 0.9567 0.9567 0.9520

0.9990 0.9990 0.9982 0.9950

0.9883 0.9891 0.9885 0.9868

freezes the NODE drift and adds a diagonal diffusion σϕ (x), a network with one hidden layer of 32 tanh units and a softplus output, fitted by the Gaussian negative log-likelihood of the frozen drift’s one-step residuals; it is rolled out with its learned noise and uses no true parameter. The NSDE’s learned noise matches the one-step noise measured with the true drift to within 0.204% per component (largest in x); in x both lie 4.33% below ς, close to the 4.8% predicted by in-step relaxation under the linear damping alone (Table 13). Both models use 8 training seeds × 10 rollout seeds; the pilot runs used to set up the pipeline are excluded from every table. Refitting the seed-0 NODE and NSDE from scratch changes every reported metric by at most 0 of its range across seeds: on the same GPU the refits and rollouts were bit-identical; reproducibility across hardware was not tested.

C.5

Exploratory: a causal transformer flow map

Transformer weather emulators such as Pangu-Weather [1], FourCastNet [23] and Aurora [2] are, in the taxonomy of Sec. 2, learned one-step flow maps rolled out autoregressively, with attention mixing spatial tokens. A three-variable system has no spatial tokens, so the analogue we can test here is attention over time. As an exploratory model, with a configuration fixed before the run and no tuning, we add a decoder-only causal transformer [28, 25, 31] that maps the last 16 states to the next-state increment (embedding dimension 64, 2 layers, 4 heads, 68,547 parameters). It is trained on the MLP-flow data with a teacher-forced one-step loss, then fine-tuned on 8-step autoregressive rollouts, as recent emulators are [17, 2], and rolled out with the same oracle noise as MLP-flow. Because the true state is fully observed and Markov, the context cannot improve the one-step prediction: the model is a sequence model with memory, like RC, not a stand-in for spatial-attention emulators. It passed the smoke-test gate fixed in advance, and all 80 rollouts over 8 training × 10 rollout seeds were scored, none diverging. Its pattern matches NODE and NSDE (Table 8). It fails invariant statistics in 79% of rollouts but the slow modes in only 19%, and the oscillatory modes in 73%, mostly through the individual near-degenerate eigenvectors, while recovering their joint subspace (merged-subspace failure rate 15%; Table 2). Of its rollouts, 17 lie inside all four invariant-statistics floors and the χ(0) band, and 2 of these fail the slow-mode test. It differs from MLP-flow in both architecture and the rollout fine-tuning, so its improvement over MLP-flow cannot be attributed to either alone. That it, like the NODE, was trained on multi-step rollouts is consistent with the training formulation, rather than the architecture, setting which properties a model gets right.

16

Related documents

• A Response Theory Probe for Learned Stochastic AI Simulators, Tested on Lorenz-63#1217566
• An Introduction to Chaotic Dynamics: Classical and quantum#726066
• Autonomous Chaotic Time Series Prediction using Physical Neuromorphic Networks#973975
• Noise-Induced Predictability Redistribution Across Forecast Horizons of Extreme Events in Chaotic Dynamics#1064421
• Learning Hierarchical Causal Representations of the Effects of Forcings on Temperature in Climate Models#1105554