Paper deep dive
A neural operator framework for data-driven discovery of stability and receptivity in physical systems
Chengyun Wang, Liwei Chen, Nils Thuerey
Intelligence
Status: succeeded | Model: Gemma-4-26B-A4B | Prompt: intel-v1 | Confidence: 95%
Last extracted: 6/21/2026, 5:41:11 AM
Summary
The paper introduces a data-driven neural operator framework designed to discover the stability and receptivity properties of complex physical systems without requiring explicit governing equations. By training a neural network as a dynamics emulator using a rollout training strategy, the authors leverage automatic differentiation to extract the system's local Jacobian. This Jacobian is then used to perform both linear stability analysis (via eigenvalue decomposition) and resolvent analysis (via weighted SVD in a reduced-order subspace). The method is demonstrated on canonical chaotic models and high-dimensional fluid flows, proving effective for identifying dominant instability modes and optimal forcing responses in strongly nonlinear regimes.
Entities (11)
Relation Signals (5)
Neural Network Emulator → enables → Linear Stability Analysis
confidence 100% · The Jacobian of the trained NN can be computed at any point in the state space via automatic differentiation... opening the door to an equation-free, data-driven modal analysis for nonlinear systems.
Neural Network Emulator → enables → Resolvent Analysis
confidence 100% · By leveraging a neural network-based emulator, we readily obtain a nonlinear representation of system dynamics while additionally retrieving intricate dynamical patterns...
Automatic Differentiation → extracts → Jacobian
confidence 100% · using automatic differentiation to extract its Jacobian, we can compute eigenmodes and resolvent modes directly from data.
Neural Operator Framework → uses → Neural Network Emulator
confidence 100% · By training a neural network as a dynamics emulator and using automatic differentiation to extract its Jacobian...
Dynamic Mode Decomposition → iscomparedto → Neural Operator Framework
confidence 90% · Furthermore, by comparing with traditional data-driven methods like DMD, we provide insights on how to select the appropriate tool for data-driven analysis.
Cypher Suggestions (0)
No Cypher suggestions yet.
Abstract
Abstract:Understanding how complex systems respond to perturbations, such as whether they will remain stable or what their most sensitive patterns are, is a fundamental challenge across science and engineering. Traditional stability and receptivity (resolvent) analyses are powerful but rely on known equations and linearization, limiting their use in nonlinear or poorly modeled systems. Here, we introduce a data-driven framework that automatically identifies stability properties and optimal forcing responses from observation data alone, without requiring governing equations. By training a neural network as a dynamics emulator and using automatic differentiation to extract its Jacobian, we can compute eigenmodes and resolvent modes directly from data. We demonstrate the method on both canonical chaotic models and high-dimensional fluid flows, successfully identifying dominant instability modes and input-output structures even in strongly nonlinear regimes. By leveraging a neural network-based emulator, we readily obtain a nonlinear representation of system dynamics while additionally retrieving intricate dynamical patterns that were previously difficult to resolve. This equation-free methodology establishes a broadly applicable tool for analyzing complex, high-dimensional datasets, with immediate relevance to grand challenges in fields such as climate science, neuroscience, and fluid engineering.
Tags
Links
- Source: https://arxiv.org/abs/2604.19465v2
- Canonical: https://arxiv.org/abs/2604.19465v2
Trouble viewing inline? Open PDF directly →
Full Text
115,523 characters extracted from source content.
Expand or collapse full text
A neural operator framework for data-driven discovery of stability and receptivity in physical systems Chengyun Wang1, Liwei Chen2 & Nils Thuerey1,∗ 1School of Computation, Information and Technology, Technical University of Munich, Boltzmannstrasse 3, D-85748 Garching bei München, Germany 2Beijing Institute of Astronautical Systems Engineering, Beijing, China ∗Corresponding author: nils.thuerey@tum.de Significance Investigating how complex systems respond to perturbations remains a central challenge in both science and engineering. Traditional analysis methods are limited to simple, linearized systems or require explicit governing equations. We introduce a general, data-driven framework that uses neural networks to automatically discover the stability and receptivity properties of any system. Our approach identifies which perturbations will grow unstably near an equilibrium state and, crucially, pinpoints the optimal way to force the system to elicit the most amplified responses, even for strongly nonlinear systems. This framework is relevant to applications in fields ranging from climate science and neuroscience to fluid engineering, particularly where first-principles models are often intractable but sufficiently rich trajectory data is available. Abstract Understanding how complex systems respond to perturbations, such as whether they will remain stable or what their most sensitive patterns are, is a fundamental challenge across science and engineering. Traditional stability and receptivity (resolvent) analyses are powerful but rely on known equations and linearization, limiting their use in nonlinear or poorly modeled systems. Here, we introduce a data-driven framework that automatically identifies stability properties and optimal forcing responses from observation data alone, without requiring governing equations. By training a neural network as a dynamics emulator and using automatic differentiation to extract its Jacobian, we can compute eigenmodes and resolvent modes directly from data. We demonstrate the method on both canonical models and high-dimensional fluid flows, employing model reduction where needed and successfully identifying dominant instability modes and input–output structures even in strongly nonlinear regimes. By leveraging a neural network–based emulator, we readily obtain a nonlinear representation of system dynamics while additionally retrieving intricate dynamical patterns that were previously difficult to resolve. This equation-free methodology establishes a broadly applicable tool for analyzing complex, high-dimensional datasets, with immediate relevance to grand challenges in fields such as climate science, neuroscience, and fluid engineering. 1 Introduction Real-world systems are inherently nonlinear and exhibit rich, multi-scale behavior in both space and time. Dynamical systems theory offers a unified mathematical framework to describe, analyze, and predict the evolution of such systems, capturing the complex interactions among quantities that evolve over time. From its foundations in the seminal work of Poincaré on celestial mechanics, this field has expanded to encompass a vast range of phenomena across the engineering, physical, and life sciences. Within this universal framework, a fundamental challenge persists: to understand and predict how a system responds to perturbations. Answering this involves determining if a system’s state is stable, identifying the conditions that trigger transitions to instability, and discovering which external inputs will most effectively steer the system’s evolution. Given that nonlinearity remains a primary challenge in analyzing and controlling dynamical systems, a leading perspective is to investigate the system’s local behavior through linearization around a fixed point to extract the dominant modes of perturbation evolution (Guckenheimer and Holmes, 2013). This strategy has given rise to a powerful suite of modal analysis techniques, many of which were developed within the hydrodynamics and cybernetics community (Simon and Mitter, 1968; Holmes, 2012; Proctor et al., 2016; Rowley and Dawson, 2017). These techniques can be broadly divided into two categories (Taira et al., 2017). The first consists of operator-based methods, such as linear stability analysis (Theofilis, 2003, 2011) and resolvent analysis (Trefethen et al., 1993; McKeon and Sharma, 2010), which are built upon the linearized operator derived from the system’s governing equations. Their profound insights, however, are contingent on this explicit, model-based knowledge. To address systems not amenable to first principles, a second category of data-driven methods has emerged. These approaches, including dynamic mode decomposition (DMD) (Schmid, 2010; Kutz et al., 2016), proper orthogonal decomposition (POD) (Lumley, 1967; Holmes, 2012), and sparse identification of nonlinear dynamics (SINDy)(Brunton et al., 2016), rely purely on observational data to uncover underlying dynamical structures, offering a pathway to analyze systems where governing equations are unknown. Linear stability analysis and resolvent analysis, as two representative operator-based methods, have been demonstrated to be highly effective in providing physical insights for system instability at an equilibrium state. The former linearizes the governing equations around the base state and assumes exponential perturbation growth. This assumption elegantly converts the linear initial-value problem into an eigenvalue problem for the Jacobian, whose eigenvalues can indicate whether the base state is stable. However, this framework is limited to describing the asymptotic perturbation evolution and does not account for short-term transient behavior (Schmid, 2007). To remedy this shortcoming, resolvent analysis was developed to investigate the linearized system’s receptivity to external forcings. Instead of focusing on internal instabilities, it can identify the optimal forcing patterns that elicit the most amplified responses, which usually represent the natural starting point for control design (Jovanović, 2021). As an input-output response framework, it provides not only a transfer-function viewpoint of analyzing the system’s receptivity characteristics, but also reveals the transient growth of perturbations, thereby explaining the transient energy amplification that arises in linearly stable systems. Therefore, linear stability analysis and resolvent analysis complement each other to provide a comprehensive picture of the system’s linearized dynamics. Despite the effectiveness of operator-based methods, their practical application to many dynamical systems of interest remains very challenging. In numerous scientific frontiers, such as neuroscience, epidemiology, and ecology, there is a fundamental lack of known physical laws from which to derive the governing equations of motion. Even in systems where such equations are well-established, such as turbulence and combustion, the construction of the linearized operator often requires specialized numerical solvers or adjoint simulations (Rolandi et al., 2024). Meanwhile, the rapid progress in experimental measurements and numerical simulations has led to data generation at an unprecedented scale. This data deluge is driving a paradigm shift to data-driven methods that can extract meaningful information directly from datasets without prior knowledge. Among which, DMD has been widely applied beyond its original domain of fluid dynamics (Schmid, 2010), in large part due to its simple framing in terms of linear regression and strong connections to nonlinear dynamical systems via Koopman spectral theory (Rowley et al., 2009). In addition, these two classical categories are being complemented by a promising hybrid perspective that seeks to perform operator-based analysis in a purely data-driven manner. A leading example is a DMD-based approach to resolvent analysis for linearly stable systems, in which DMD is applied to time-resolved flow snapshots to approximate the system’s eigenmodes. The resolvent operator is then constructed within the reduced-order subspace spanned by these DMD modes (Herrmann et al., 2021). This idea is further refined by incorporating physics-informed constraints to restrict the admissible DMD model space to certain matrix manifolds that preserve the desired physical properties (Baddoo et al., 2023). However, the linear assumption of DMD becomes increasingly restrictive as nonlinear effects grow. A notable recent improvement is the linear and nonlinear disambiguation optimization (LANDO) framework (Baddoo et al., 2022), which uses kernel learning to identify nonlinear dynamics while isolating an interpretable linear component, and has already shown promising results in hydrodynamic stability applications (Hernández et al., 2023). Neural networks (NNs) offer a pathway beyond these limitations. As universal function approximators, N emulators can learn complex nonlinear dynamics from spatiotemporal measurements of the system, serving as surrogate models that capture both linear and nonlinear behaviors (Floryan and Graham, 2022; Yu and Wang, 2024; Page et al., 2024). They can be trained either in a purely data-driven manner (Morton et al., 2018; Thuerey et al., 2020; Chen and Thuerey, 2023) or augmented with physical information (Bar-Sinai et al., 2019; Raissi et al., 2019; Kochkov et al., 2021; List et al., 2022; Chen and Thuerey, 2024). Crucially, these emulators, which can be viewed as learned approximations of classical numerical methods such as finite difference (Koehler et al., 2024) or pseudo-spectral schemes (Li et al., 2020; McCabe et al., 2023), are fully differentiable. This property is the key to our approach: the local Jacobian of the trained N can be computed at any point in the state space via automatic differentiation (Chen and Thuerey, 2024), opening the door to an equation-free, data-driven modal analysis for nonlinear systems. To this end, we propose a data-driven approach to learn a N emulator from data generated by truly nonlinear systems, with a schematic illustration presented in Fig. 1. The Jacobian of the learned N emulator serves as an approximation of the system’s local linear operator, which is subsequently used to perform operator-based modal analysis. To demonstrate the power and generality of this framework, we diagnose four representative dynamical systems of increasing complexity, successfully extracting their essential stability characteristics and input-output resolvent modes. For high-dimensional systems, we additionally introduce a reduced-order formulation that makes the learning process more efficient while retaining the leading spectral information. Furthermore, by comparing with traditional data-driven methods like DMD, we provide insights on how to select the appropriate tool for data-driven analysis. 2 Methods Considering a forced nonlinear dynamical system (temporally continuous but spatially discretized into N degrees of freedom), the governing equation can be expressed in terms of the state variable ∈ℂN q ^N in the compact operator form as ddt=()+, d qdt=N( q)+ f, (1) where N is the nonlinear operator collecting all the terms except for the time derivative and ∈ℂN f ^N is an external forcing term. To linearize this ordinary differential equation, the state variable can be decomposed as a sum of base state b∈ℂN q_b ^N and infinitesimal perturbation ′∈ℂN q ^N, which yields dbdt+d′dt=(b+′)+=(b)+∂|b′+O(|′|2)+. d q_bdt+ d q dt=N( q_b+ q )+ f=N( q_b)+ ∂ q |_ q_b q +O(| q |^2)+ f. (2) While the base state could be chosen as a time-averaged solution (Edstrand et al., 2016; Barkley and Henderson, 1996), we focus on the stability around steady-state solutions, and thus b q_b is an equilibrium point of Eq. 1 satisfying db/dt=(b)=0d q_b/dt=N( q_b)=0. After neglecting the higher-order infinitesimals in Eq. 2, we obtain the linearized equation of the perturbation ′ q as d′dt=∂|b′+=′+, d q dt= ∂ q |_ q_b q + f=A q + f, (3) where ∈ℂN×NA ^N× N represents the discretized linear operator of this continuous dynamical system at base state b q_b, i.e., the local Jacobian. Based on this Jacobian A, operator-based modal analyses such as linear stability and resolvent analysis can be performed; for details see SI Appendix, section 1. Neural network operator Figure 1: Schematic of our data-driven algorithm. (1) Time-resolved state snapshots collected from a dynamical system are used. (2) An N emulator is trained using a rollout strategy to accurately learn the system’s evolution map. (3) The Jacobian of the trained emulator is then extracted via automatic differentiation and used to approximate the local Jacobian. This data-driven operator is the key ingredient for (4) N-based modal analysis to perform both linear stability analysis via an eigenvalue decomposition to find unstable modes, and for resolvent analysis via a weighted SVD on the projected resolvent to identify optimal forcing modes, response modes, and the gain distribution. We adopt a time-stepping framework using a N emulator to evolve the original unforced dynamical system of Eq. 1. The state n q_n at time t is mapped to the subsequent state n+1 q_n+1 over a time increment Δt t as: n+1=fθ(n), q_n+1=f_θ( q_n), (4) where fθf_θ is an arbitrary N architecture with weights θ. Let P represent a classical numerical solver that also advances the system state with a time step Δt t. Then the learning task is to find the optimal weights θopt _opt of the N to best approximate the numerical solver, i.e., fθ≈Pf_θ≈ P. When trained on a sufficiently diverse dataset, the N emulator can capture the underlying latent dynamics and generalize effectively to previously unseen system states. Since the network is trained by minimizing the difference between predicted and true trajectories, their interplay during the learning process is crucial. A straightforward way of training is one-step training, which only minimizes the one-step prediction error (Tran et al., 2021). By contrast, rollout training by autoregressively applying fθf_θ for multiple time steps can learn long-term dependencies and correct for compounding errors (Um et al., 2020). This can significantly improve the long-term prediction performance, especially when dealing with high-dimensional nonlinear dynamical systems. Suppose the ground truth trajectories produced by P are 0(j),1(j),⋯,m(j) q^(j)_0, q^(j)_1,·s, q^(j)_m, where the superscript j∈1,⋯,nj∈ \1,·s,n \ indexes the trajectories corresponding to n distinct initial conditions, and m is the maximum sampling time step per trajectory (Part 1, Fig. 1). Given an arbitrary state i(j) q^(j)_i extracted from the ground truth datasets, the trajectories predicted by iteratively applying the learned model fθf_θ for t steps are denoted as i(j),fθ1(i(j)),⋯,fθt(i(j)) q^(j)_i,f^1_θ( q^(j)_i),·s,f^t_θ( q^(j)_i). The total loss function for training with t-step rollouts averaged on n trajectories is: ℒ(θ)=1n∑j=1n∑i=0m−t∑k=1t[ℒ2(fθk(i(j)),i+k(j))],L(θ)= 1n _j=1^n _i=0^m-t _k=1^t [L_2 (f_θ^k( q^(j)_i), q^(j)_i+k ) ], (5) where ℒ2L_2 represents L2-norm loss function. Full gradients can be computed owing to the differentiability of the N emulator, which enables backpropagation through t rollout steps and subsequent updates of the network weights θ (Part 2, Fig. 1). Relationship between N-based and operator-based Jacobians Since the N emulator operates in a temporally discrete manner, whereas operator-based analysis relies on temporally continuous equations, an essential step is to bridge the gap between these two representations of the dynamical system. Here, we can also linearize the temporally discrete dynamical system (Eq. 4) based on the base-perturbation decomposition: b+n+1′=fθ(b+n′), q_b+ q _n+1=f_θ( q_b+ q _n), (6a) b+n+1′=fθ(b)+∂fθ∂|bn′+O(|n′|2). q_b+ q _n+1=f_θ( q_b)+ ∂ f_θ∂ q |_ q_b q _n+O(| q _n|^2). (6b) Since the base state solution b q_b satisfies b=fθ(b) q_b=f_θ( q_b), we can eliminate both terms in Eq. 6b, and neglect the high-order term to obtain: n+1′=∂fθ∂|bn′=n′, q _n+1= ∂ f_θ∂ q |_ q_b q _n=N q _n, (7) where the Jacobian N of our learned neural operator fθf_θ can be calculated via automatic differentiation evaluated at b q_b (Part 3, Fig. 1). The differentiability of the N emulator actually allows for the use of automatic differentiation to evaluate N at any input state. In the absence of external forcing, the general solution of Eq. 3 takes an exponential form, and thus two successive perturbation states ′ q should satisfy: n+1′=exp(Δt)n′. q _n+1= (A t) q _n. (8) Hence, we can derive the relationship between the N-based Jacobian N and the operator-based Jacobian A as =exp(Δt).N= (A t). (9) Benefiting from this relationship, if a learned surrogate model can accurately capture the underlying system dynamics, it can then be utilized to evaluate the stability of the original system without knowing its governing equations. N-based linear stability and resolvent analysis According to Eq. 9, the eigenvalues of the N-based and operator-based Jacobians, denoted by the diagonal matrices _N and _A, respectively, are expected to satisfy the exponential relationship: =exp(Δt) _N= ( _A t). The associated eigenvectors corresponding to each eigenvalue should be identical, i.e., =V_N=V_A. Provided the learned model accurately approximates the original system, the extracted eigenvalues and eigenvectors from N serve as good approximations of those from A, whose spectral characteristics can thereby be inferred indirectly (Part 4a, Fig. 1). This constitutes the core procedure of N-based linear stability analysis. In essence, the entire workflow can be interpreted as a data-driven eigendecomposition, enabling stability assessment and extraction of dominant unstable modes without requiring explicit governing equations. As for resolvent analysis, a natural idea is to replace the Jacobian A by log()/Δt (N)/ t to compute the resolvent operator (ω)=(−iω−)−1H(ω)=(-i -A)^-1. However, it is not advisable to directly take the logarithm of the Jacobian N, as the logarithmic function is a multivalued function in the complex domain, and as such, is strongly underdetermined. N typically corresponds to a large number of A that satisfy Eq. 9. Instead, we can leverage the results of the linear stability analysis to implement the dimensionality reduction (Reddy and Henningson, 1993). In a nutshell, we can project the forced system Eq. 3 into the subspace expanded by a truncated r eigenvectors (Schmid and Henningson, 2012; Herrmann et al., 2021). Supposing the original state variable ′ q and f in Eq. 3 can be expressed in the subspace spanned by the first r linearly independent eigenvectors ,rV_A,r, we have d,rdt=,r+,r, dV_A,r xdt=AV_A,r x+V_A,r y, (10) where =[x1,x2,⋯,xr]T∈ℂr x=[x_1,x_2,·s,x_r]^T ^r and =[y1,y2,⋯,yr]T∈ℂr y=[y_1,y_2,·s,y_r]^T ^r are the vectors of expansion coefficients in the eigenvector coordinate. Based on the definition of eigenvalue ,r=,r,rAV_A,r=V_A,r _A,r, we can obtain the new governing equation in the eigenvector coordinate: ddt=,r+. d xdt= _A,r x+ y. (11) Given then weighting matrix Q used to define the state norm in the original coordinate as ‖^‖2=^∗^|| q||^2_Q= q^*Q q, the corresponding weighting matrix in the new coordinate becomes ~=,r∗,r Q=V_A,r^*QV_A,r, which can be factorized as ~=~∗~ Q= F^* F via Cholesky factorization. Under this new coordinate system, we can carry out a weighted SVD of the projected resolvent: ~(−iω−,r)−1~−1=~(ω)(ω)~∗(ω), F(-i - _A,r)^-1 F^-1= _ F(ω) (ω) ^*_ F(ω), (12) where the diagonal entries of the singular value matrix Σ(ω) (ω) are the resolvent gains and the column vectors in ~(ω) _ F(ω) and ~(ω) _ F(ω) are the corresponding response and forcing modes at frequency ω. The conversion of resolvent modes back to the original space is then performed via =,r~−1~(ω) =V_A,r F^-1 _ F(ω) and =,r~−1~(ω) =V_A,r F^-1 _ F(ω). In practice, the derivation of ,r _A,r and ,rV_A,r requires the knowledge of the corresponding governing equation, so ,r _N,r and ,rV_N,r acquired from data are employed instead (Part 4b, Fig. 1). However, it is important to avoid indiscriminately increasing the dimensionality of the subspaces (Herrmann et al., 2021). A reliable projection of a high-dimensional system onto a reduced subspace requires that the chosen eigenvectors be linearly independent. But especially in non-normal systems, where eigenvectors are generally non-orthogonal, care must be taken not to select an overly large number of modes r. An excessively large r may lead to a redundant coordinate system, resulting in a non-unique representation of the system state in terms of eigenvector coordinates. This will deteriorate the validity of resolvent operator estimation instead of enriching the basis. Stability in reduced-order models For high-dimensional systems, it is natural to first represent the state snapshots in a POD basis and then learn the dynamics in the correspondingly smaller space of coefficients. Collecting n snapshots from m trajectories into a matrix ∈ℂN×mnQ ^N× mn, and noting that typically mn≪Nmn N, we employ the method of snapshots (Sirovich, 1987) to obtain the compact decomposition =∗Q=U W^*, where ∈ℂN×lU ^N× l has orthonormal columns and spans the POD subspace, with l=rank()≤mnl=rank(Q)≤ mn. Hence, every snapshot n q_n in the dataset can be represented exactly as n=∗n,n=n. a_n=U^* q_n, q_n=U a_n. (13) According to Eq. 4, projecting the dynamics onto the POD coordinates gives the evolution operator gθg_θ for the induced reduced dynamics n+1=gθ(n)=∗fθ(n). a_n+1=g_θ( a_n)=U^*f_θ(U a_n). (14) Let b q_b be a base state and b=∗b a_b=U^* q_b its POD representation. The Jacobian of the reduced dynamics PODN_POD evaluated at b a_b is POD=∂gθ∂|b=∗∂fθ∂|b=∗,N_POD= ∂ g_θ∂ a |_ a_b=U^* ∂ f_θ∂ q |_ q_bU=U^*NU, (15) where N is the aforementioned original N-based Jacobian. Since U defines an exact coordinate transform on the POD subspace, PODN_POD is the representation of N in the POD coordinates restricted to this subspace. Now let ,s=[1,⋯,s]V_N,s=[ v_1,·s, v_s] denote the dominant eigenvectors associated with the leading s eigenvalues collected in ,s _N,s, such that ,s=,s,s.NV_N,s=V_N,s _N,s. (16) If the POD subspace contains this leading invariant subspace, span(,s)⊆range()span(V_N,s) (U), then there exists a full-column-rank matrix ∈ℂl×sT ^l× s such that ,s=.V_N,s=UT. (17) Substituting this into the eigenvalue relation and left-multiplying by ∗U^* yields POD=,s.N_PODT=T _N,s. (18) Therefore, the reduced dynamics expressed in the full POD coordinates preserve the eigenvalues of the original dynamics on the POD subspace, and the corresponding eigenvectors are related through the POD projection and lifting. In practice, retaining a suitably chosen number of POD modes is usually adequate to recover the leading eigenvalues accurately, with the required number depending on the dimensionality of the specific problem. The reason is that POD ranks the modes according to their contribution to the snapshot energy, while the dominant linear structures are typically the most weakly damped or most amplified components and therefore dominate the variance of the trajectory data. As a result, these dynamically most important structures are preferentially captured by the leading POD modes, so that mode truncation usually has little effect on the dominant linear characteristics of interest. From a mathematical standpoint, the POD‑reduced system provides a Ritz approximation of the eigenspectrum of the original high‑dimensional system. Hence, learning the truncated POD coefficient dynamics is sufficient to preserve the dominant linear stability characteristics while substantially reducing the cost of training and analysis. 3 Results We present four examples of different complexity, ranging from the Lorenz system and the complex Ginzburg-Landau equation to two- and three-dimensional channel flows, which were chosen to highlight different features of the N-based modal analysis. To ensure a robust Jacobian estimation, we adopt a model ensemble strategy (Chua et al., 2018; Allen et al., 2022) by averaging the Jacobians predicted from five trained models initialized with different random seeds. This can not only mitigate the effects of weight initialization diversity but also lead to a more representative average gradient. The network architecture settings for each case, plotted in SI Appendix Fig. S9, were selected following the systematic experiments conducted in our prior benchmark study (Koehler et al., 2024). In addition, several ablation studies on data quality, as well as additional application examples, are also provided in SI Appendix due to space constraints. Lorenz-63 Base State b q_b N-based Jacobian N Analytical Results exp(Δt) (A t) equilibrium points [000] bmatrix0\\ 0\\ 0 bmatrix [0.91850.0947−0.00010.26571.00360.00050.0002−0.00010.9735] bmatrix0.9185&0.0947&-0.0001\\ 0.2657&1.0036&0.0005\\ 0.0002&-0.0001&0.9735 bmatrix [0.91790.095100.26631.00350000.9737] bmatrix0.9179&0.0951&0\\ 0.2663&1.0035&0\\ 0&0&0.9737 bmatrix [626227] bmatrix6 2\\ 6 2\\ 27 bmatrix [0.90540.0945−0.00410.00740.9861−0.08360.07890.08800.9705] bmatrix0.9054&0.0945&-0.0041\\ 0.0074&0.9861&-0.0836\\ 0.0789&0.0880&0.9705 bmatrix [0.90520.0946−0.00410.00600.9869−0.08320.08000.08730.9700] bmatrix0.9052&0.0946&-0.0041\\ 0.0060&0.9869&-0.0832\\ 0.0800&0.0873&0.9700 bmatrix [−62−6227] bmatrix-6 2\\ -6 2\\ 27 bmatrix [0.90500.09480.00390.00710.98650.0835−0.0812−0.08690.9697] bmatrix0.9050&0.0948&0.0039\\ 0.0071&0.9865&0.0835\\ -0.0812&-0.0869&0.9697 bmatrix [0.90520.09460.00410.00600.98690.0832−0.0800−0.08730.9700] bmatrix0.9052&0.0946&0.0041\\ 0.0060&0.9869&0.0832\\ -0.0800&-0.0873&0.9700 bmatrix non-equilibrium points [1.46898.7045−7.6905] bmatrix1.4689\\ 8.7045\\ -7.6905 bmatrix [0.92210.09410.00010.33881.0075−0.01800.08550.02380.9734] bmatrix0.9221&0.0941&0.0001\\ 0.3388&1.0075&-0.0180\\ 0.0855&0.0238&0.9734 bmatrix [0.92150.0952−0.00070.33931.0071−0.01450.08470.01870.9736] bmatrix0.9215&0.0952&-0.0007\\ 0.3393&1.0071&-0.0145\\ 0.0847&0.0187&0.9736 bmatrix Table 1: Comparison of the N-based Jacobian matrix N and the operator-based Jacobian matrix A around three equilibrium points and one non-equilibrium point (the initial point in Fig. 2). We first consider the canonical Lorenz system as a simple model based on ordinary differential equations (ODEs). Although it is a chaotic system, NNs have been proven capable of characterizing its temporal dynamics, including accurately advancing the solution in time for unseen initial conditions (Brunton and Kutz, 2022). The state of the Lorenz system given by =[x,y,z]T q=[x,y,z]^T is governed by dxdt dxdt =σ(y−x), =σ (y-x ), (19) dydt dydt =x(ρ−z)−y, =x (ρ-z )-y, dzdt dzdt =xy−βz. =xy-β z. This system has three equilibrium points, one at the origin ([0,0,0]T[0,0,0]^T) and two symmetric off-origin points for ρ>1ρ>1 ([±β(ρ−1),±β(ρ−1),ρ−1][± β(ρ-1),± β(ρ-1),ρ-1]). By applying a local linearization to the Lorenz system around a base state b=[xb,yb,zb]T q_b=[x_b,y_b,z_b]^T, the analytical local Jacobian is obtained: =∂|b=[−σ0ρ−zb−1−xbybxb−β].A= ∂ q |_ q_b= bmatrix-σ&σ&0\\ ρ-z_b&-1&-x_b\\ y_b&x_b&-β bmatrix. (20) Here, the standard parameter values of σ=10σ=10, ρ=28ρ=28, and β=8/3β=8/3 are chosen, and thus the system exhibits chaotic motion confined to a strange attractor. The training data is constructed from high-fidelity simulations of the Lorenz system via the fourth-order Runge–Kutta method. 100100 different initial conditions are advanced in time for a total of m=800m=800 snapshots with a fixed 0.010.01 time units. The N architecture is a three-layer network with 10 neurons per hidden layer adopted from (Brunton and Kutz, 2022). Figure 2: Time evolution of the true Lorenz system (red solid) compared with the N emulator prediction (green dotted) from an unseen initial point (blue circle). (a) Trajectories in phase space; (b) Component-wise time evolution series. Even trained on the trajectory dataset with a one-step rollout strategy, it can be used to autoregressively generate new trajectories as far into the future as desired, even for initial conditions not present in the training dataset. Fig. 2(a) presents the trajectory comparison between the N emulator’s prediction and ground truth starting from a randomly selected initial condition over an interval of 8 time units. Fig. 2(b) offers a clearer view of the approximation’s performance by showing the time evolution of three components separately. The close agreement between both simulations in the early stage suggests that the trained N emulator captures the system dynamics with high fidelity, effectively mimicking the update rule of a classical Runge–Kutta scheme over a time increment of 0.010.01. Despite a noticeable deviation observed for the N emulator at the end of the simulation, the performance remains commendable given the lack of rollout steps at training time. Meanwhile, we leverage the automatic differentiation to compute the N-based Jacobians around the three equilibrium points and compare them with the analytical Jacobians derived from Eq. 20 and Eq. 9. As shown in Tab. 1, the N-based and analytical Jacobians exhibit an almost perfect element-wise match. This suggests that the following modal analysis performed using the N-based Jacobian can reproduce theoretical results, despite the fact that the N itself is trained purely on data and is agnostic to the underlying governing equations. This paves the way for data-driven modal analysis approaches for both linear and nonlinear scenarios in higher-dimensional settings. Furthermore, we found that for this simple ODE system, the Jacobian agreement is not limited to equilibrium points but instead extends to arbitrary points along the trajectory (see the last row of Tab. 1). This indicates that the trained neural operator can also faithfully represent local dynamical behavior at any points and may facilitate tasks such as computing finite-time Lyapunov exponents (FTLE) (Shadden et al., 2005). This is not the primary focus of our investigation. Complex Ginzburg-Landau system The second example considered is the one-dimensional complex Ginzburg–Landau equation (Bagheri et al., 2009). It is widely recognized as a reduced-order model for studying diverse instability behaviors in spatially evolving flows by varying parameter configurations. As its linearized operator can be obtained analytically through the manipulation of the differential matrices, it has become a benchmark for the evaluation and development of stability analysis techniques (Chen and Rowley, 2011; Herrmann et al., 2021; Martini et al., 2021). The nonlinear complex Ginzburg–Landau equation defined on an infinite interval x∈(−∞,∞)x∈(-∞,∞) is written as ∂t=(−ν∂x+γ∂2∂x2+μ(x))−a||2=−a||2, ∂ q∂ t= (-ν ∂ x+γ ∂^2∂ x^2+μ(x) ) q-a| q|^2 q\\ =A q-a| q|^2 q, (21) where −a||2-a| q|^2 q is the cubic nonlinear term. Since =0 q=0 is an equilibrium point, the linearized system around this point simply involves omitting the cubic term and thus A is the linearized operator. The advection and dispersion properties are controlled by the complex terms ν=U+2icuν=U+2ic_u and γ=1+icdγ=1+ic_d, respectively. The real-valued term μ(x)=μ0−cu2+μ2x2/2μ(x)= _0-c_u^2+ _2x^2/2 is defined as a quadratic function to model exponential instabilities (Bagheri et al., 2009). A detailed description of the parameters can be found in Tab. 2. The homogeneous boundary conditions are considered at x→±∞x→±∞, so that perturbations are allowed to grow and decay throughout the entire domain. Variable Description Value U Mean advection velocity 2.0 cuc_u Most unstable wavenumber 0.2 cdc_d Dispersion parameter -1.0 μ0 _0 Bifurcation parameter 0.23(0.41) μ2 _2 Degree of non-parallelism -0.01 a Nonlinearity 0(1) Table 2: Parameter descriptions and settings for the complex Ginzburg–Landau system. Values corresponding to the nonlinear scenario are indicated in parentheses. Figure 3: Comparison of the real part of the Jacobians between the N-based operator and the operator-based ground truth. The first and second rows correspond to the linear and nonlinear scenarios, respectively. We aim to demonstrate the applicability of the N emulator approach to both linear and nonlinear regimes even in complex number space. For this purpose, we consider two datasets derived respectively from the linear (Herrmann et al., 2021) and nonlinear (Bagheri et al., 2009) forms of this system. The sample trajectories and the specific parameter settings used in each scenario are presented in SI Appendix, Fig. S7 and Tab. 2. For both scenarios, we conduct 30 independent simulations, each initiated with distinct initial conditions derived from the superposition of multiple random Gaussian profiles. For each simulation, m=100m=100 snapshots are collected at regular time intervals of 0.50.5. We employ a complex-valued, six-layer perceptron with 512512 neurons per hidden layer. A truncation level of r=24r=24 is used for the eigenbasis in this example. All results are obtained using a one-step rollout training strategy. More details about this case are provided in SI Appendix, section 2. As shown in Fig. 3, we compare the Jacobian heatmaps of the learned emulator with those of the analytical operator for both linear and nonlinear scenarios. We can find that the network perfectly recovers the operator in the linear scenario, whereas for the nonlinear scenario, it successfully captures the main features but exhibits minor discrepancies in finer details. These deviations are most likely caused by nonlinear mappings being inherently more difficult to approximate than linear ones. Figure 4: Linear stability analysis results of the linear and nonlinear Ginzburg-Landau systems. (a)(b) Eigenspectra obtained from the N-based Jacobian (∘ ), DMD-based Jacobian (++) and operator-based ground truth (∙ [rgb]0,0,1 [named]pgfstrokecolorrgb0,0,1 ). The dashed line denotes the stability boundary in the complex plane. (c)(d) The leading three eigenmodes of the N-based Jacobians in both scenarios, where solid and dashed lines show the real part and magnitude of the modes. The thick gray lines in the background show operator-based ground-truth for comparison. In the context of linear stability analysis, Fig. 4(a) and (b) compare the eigenspectra of the N-based and operator-based Jacobians for linear and nonlinear scenarios, respectively. For reference, the eigenspectra derived from the DMD model are also included. Fig. 4(c) and (d) show the leading three eigenmodes identified by the N emulator for both scenarios. For the linear scenario, it is anticipated that the DMD-based results show excellent agreement with the analytical solution, as DMD is fundamentally designed to approximate linear operators. Even so, the N-based results also show excellent alignment of the leading eigenvalues and eigenvectors with the analytical reference. For the nonlinear scenario, DMD-based approach totally breaks down, as demonstrated by the identification of more than one unstable mode located in the right half of the complex plane (see Fig. 4(b)). In contrast, the N-based Jacobian still yields an accurate leading unstable eigenvalue, even though it differs in finer structures according to Fig. 3. This suggests that the N emulator can effectively identify the most dominant mode regardless of the dynamical properties of datasets. While the remaining higher-order modes deviate more noticeably, they decay rapidly and contribute minimally to the long-term behavior of perturbations. Figure 5: Resolvent analysis results of the complex Ginzburg–Landau systems. The first and the second row correspond to the linear and nonlinear scenarios, respectively. (a)(c) Resolvent gain distribution for the first three modes with respect to the frequency. (b)(d) The first three forcing and response modes at the peak gain frequency. The thick gray lines in the background show operator-based ground-truth for comparison. In the context of resolvent analysis, Fig. 5(a) and (c) show the gain curves of the first three resolvent modes as a function of the forcing frequency for the two scenarios. The corresponding dominant response and forcing modes at the peak gain frequency are presented in Fig. 5(b) and (d). For the linear scenario, our N-based approach can yield highly accurate gain curves and mode shapes. For the nonlinear scenario, it also succeeds in predicting the gain curve near the dominant frequency, and accurately reconstructs the response and forcing modes corresponding to this frequency. Though applicable to both linear and nonlinear regimes, the higher-order mode shapes predicted by the N emulator become increasingly influenced by noise and deviate from ground truth more noticeably, as shown in Fig. 5(c) and (d). There are two important aspects to note for the gain distributions of Fig. 5(c). First, the N emulators have worse performance for gain prediction at frequencies far from the peak gain. This phenomenon can be interpreted via dyad expansion of the resolvent operator (Symon et al., 2018). The resolvent operator H(ω)H(ω) can be expressed as a sum of dyadic products of left eigenvectors ~j g_j and right eigenvectors ~j h_j, each weighted by the distance between the eigenvalue and the imaginary axis 1/(iω−λj)1/(iω- _j) as H(ω)=(iω−)−1=∑j=1n1iω−λj~j~j∗≈1iω−λ1~1~1∗,H(ω)=(iω I- Q)^-1= _j=1^n 1iω- _j g_j h_j^*≈ 1iω- _1 g_1 h_1^*, (22) where we can see that near the frequency of maximal gain, i.e., the forcing frequency ω is close to the imaginary part of the dominant eigenvalue, the resolvent operator admits a good rank-1 approximation. Since the leading eigenmode is accurately recovered (see Fig. 4), the corresponding resolvent operator is also highly accurate for those frequency intervals satisfying the rank-1 approximation. Second, the N-based method tends to exhibit more noise in the predicted forcing modes compared to the response modes. However, this behavior is in line with other data-driven approaches and can be attributed to the same cause: in non-normal systems, the spatial structure of the direct eigenvectors resembles that of the response modes, but can differ significantly for the forcing modes (Herrmann et al., 2021). 2D transitional channel flow We then consider the two-dimensional plane Poiseuille flow with Re = 2000 (based on the centerline velocity and the channel half-height) governed by the incompressible Navier-Stokes equations. The analytical derivation of the linearized operator for this example remains difficult, but by considering perturbations in the wavenumber space, the problem can be effectively reduced to the Orr–Sommerfeld (OS) equation (Schmid and Henningson, 2012). The dynamics of the velocity perturbations is governed by the OS operator, which we use for the ground truth reference (Details in SI Appendix, section 3). In this case, the N emulator aims to model the temporal evolution of infinitesimal perturbations around the base parabolic velocity profile U(y)=1−y2U(y)=1-y^2, and hence, snapshots of the velocity perturbation field in two spatial directions are collected during the transient evolution stage for training. We employ a spectral solver (Mortensen, 2017, 2018) on a computational domain of 2π×22π× 2 in the streamwise and wall-normal directions, where periodicity is imposed along the streamwise direction. The spatial domain is discretized by Ny=64N_y=64 Chebyshev collocation points in the wall-normal direction and Nx=32N_x=32 Fourier collocation points in the streamwise direction. To simplify the analysis, the flow perturbations with a single wavenumber kx=1k_x=1 are advanced in time with a time step of 0.0020.002 time units. A total of 6060 trajectories are simulated with m=400m=400 snapshots saved every 0.50.5 time units. The initial velocity perturbation fields are of the form u(x,y) u(x,y) =ϵ⋅ℜ[ϕ′(y)⋅eikxx], =ε· [φ (y)· e^ik_xx ], (23) v(x,y) v(x,y) =−ϵ⋅ℜ[i⋅ϕ(y)⋅eikxx] =-ε· [i·φ(y)· e^ik_xx ] where ϕ(y)φ(y) is the initial profile constructed as a superposition of multiple Chebyshev‑polynomial profiles with random amplitudes, with special treatment at the wall boundaries to satisfy both homogeneous Dirichlet and Neumann boundary conditions. And ϵε serves as a scaling factor for the initial perturbation magnitude, determining the strength of the nonlinear effects. Although the dataset is normalized prior to training, the specific value of ϵε has a strong influence on the fluid behavior during the transient phase (see SI Appendix, Fig. S8). We consider a weakly nonlinear dataset with ϵ=10−6ε=10^-6 and a strongly nonlinear dataset with ϵ=10−1ε=10^-1. For network training, we employ a classic U-Net architecture (Ronneberger et al., 2015) with four encoding and decoding layers. We set the number of channels in the top layer to 3232 channels for optimal trade-offs between runtime efficiency and predictive performance. A truncation level of r=30r=30 is used for the eigenbasis of the projection subspace, and ten steps of temporal unrolling are used at training time. Figure 6: Comparison of the eigenspectra obtained from the N-based Jacobian (∘ ), DMD-based Jacobian (++) and operator-based ground truth (∙ [rgb]0,0,1 [named]pgfstrokecolorrgb0,0,1 ) for the 2D channel flow system. The dashed line denotes the stability boundary in the complex plane. (a) Weakly nonlinear dataset. (b) Nonlinear dataset. Figure 7: Comparison of the first-order A and P eigenmodes of the streamwise velocity field for the 2D channel flow system. The first column is the reference results obtained from the OS operator. The second and third columns are the N-based Jacobian obtained from weakly nonlinear and nonlinear datasets. For the chosen Reynolds number of 2000, the channel system is stable (Schmid and Henningson, 2012), which can be verified by the eigenspectrum of the OS operator from Fig. 6. All operator-based eigenvalues fall in the stable complex half-plane and they are located on three branches forming a Y-shaped pattern. These branches have been labeled A (cr→0c_r→ 0), P (cr→1c_r→ 1), and S (cr≈2/3c_r≈ 2/3) by (Mack, 1976). The branches describe distinct features of the perturbation dynamics: P modes (center modes) on the P-branch with much higher phase speeds capture the fast dynamics in the center of the channel, and A modes (wall modes) with rather small phase velocities capture the slow dynamics near the channel wall. The S modes approaching a phase speed 2/32/3 are highly damped. In Fig. 6, we note that the first leading eigenvalues of the A-branch and P-branch agree well between the results of the OS operator and the N emulator trained on both weakly nonlinear and nonlinear datasets. This indicates that the N emulator can successfully identify the two most dominant perturbation structures near the wall and centerline, i.e., the least-damped A and P modes (see Fig. 7), despite the strong nonlinearities of the dataset. While the remaining higher-order A and P modes appear to be spurious, and the S modes with high decay rates were not recovered, their contribution to the system’s dynamics is negligible. We attribute this to limitations in the numerical precision of the dataset, as well as the inherent noise in the trained emulator, which potentially obscures the evolution of these highly damped modes. Additionally, the randomly dispersed DMD-based eigenvalues in Fig. 6(b) again highlight its inherent flaws when handling nonlinear datasets. It is noteworthy that another spurious eigenvalue appears on the imaginary axis close to the stability boundary for the case of the nonlinear dataset (Fig. 6(b)), which is absent for the operator-based ground truth and the case of the weakly nonlinear dataset. This discrepancy may be attributed to the preprocessing applied to the training data. Given the negligible magnitude of streamwise perturbations, the raw data were preprocessed by subtracting the base flow and subsequently normalized. While such treatment is reasonable in the linear or weakly nonlinear regime, it can cause distortions in nonlinear scenarios where the linear superposition principle does not hold. Consequently, zero-frequency modes associated with non-oscillatory components were introduced in the dataset. These modes are subsequently identified by the N emulator during feature extraction, as evidenced by their appearance along the imaginary axis in Fig. 6(b). Figure 8: Resolvent analysis results of the 2D channel flow system. (a) Resolvent gain curves of the analytical operator and the N-based resolvent operators derived from weakly nonlinear and nonlinear datasets. (b) The first-order forcing and response modes of the streamwise velocity field at the ω=0.31ω=0.31. As for the resolvent analysis, Fig. 8(a) displays the gain curves of three distinct operators in response to harmonic forcings with frequencies in the range of [0,0.5][0,0.5]. The peak frequency at ω=0.31ω=0.31 is accurately predicted by the N emulator despite the nonlinear nature of this scenario. We note that the corresponding resolvent gain curve exhibits an extra peak at zero frequency for the nonlinear scenario. This phenomenon aligns with the presence of non-oscillatory spurious modes discussed above, and the additional gain peak comes from resonating with them. In our tests, it can be substantially mitigated either by filtering the corresponding spurious modes in the reconstructed spectrum or by training on the wall-normal velocity dataset alone, for which no such subtraction is required. Nevertheless, the N emulator accurately predicts the physically relevant peak frequency of ω=0.31ω=0.31 originating from a near-resonance with the first A-branch eigenmode for both scenarios. The resolvent modes of streamwise and wall-normal velocity perturbations at ω=0.31ω=0.31 are illustrated in Fig. 8(b), revealing that the N emulator detects the area that is most sensitive for perturbation amplification for channel flow at the current Reynolds number. Since the leading eigenvectors of the N-based Jacobian align closely with those of the true system (Fig. 7), it naturally leads to accurate resolvent mode predictions. Similar to the observations in the Ginzburg–Landau case, the reconstructed response mode shows near-exact agreement with the ground truth, whereas the reconstructed forcing mode yields a less accurate approximation, which is also attributed to the non-normal nature of the channel flow. This behavior is more prominent in the context of nonlinear datasets, indicating that nonlinearity may exacerbate the challenges associated with training under the current neural network architecture. Reduced 3D transitional channel flow To further demonstrate the applicability of our method to high‑dimensional systems, we consider the three‑dimensional plane Poiseuille flow, which is closer to realistic flow configurations due to the complex dynamics arising in three spatial dimensions. Although a U‑Net can in principle learn the snapshot‑to‑snapshot mapping directly in physical space, for high‑dimensional problems the corresponding Jacobian matrix has a dimension equal to the total number of degrees of freedom of all flow variables over the full discretized field. The resulting matrix makes direct eigendecomposition and the subsequent construction of the resolvent operator computationally expensive. Hence, reduced-order models are attractive to make such problems tractable. Here, a typical linear compression method, proper orthogonal decomposition (POD), is applied to the stored data snapshots to obtain the orthogonal spatial basis and corresponding temporal coefficient trajectories. Then, we train the neural emulator on these coefficient trajectories and apply the following modal analysis in an encoder-decoder approach. Figure 9: Comparison of the (a) eigenspectra and (b) resolvent gain curves from the operator-based Jacobian and N-based Jacobian for the 3D channel flow system. Figure 10: Comparison of the first-order A and P eigenmodes of the streamwise velocity field for the 3D channel flow system. The problem setting is in line with that used in the two-dimensional case. The computational domain is now extended to 2π×2×22π× 2× 2, where an additional periodic spanwise direction is discretized using Nz=32N_z=32 Fourier collocation points. The streamwise and spanwise wavenumbers are fixed at kx=kz=1k_x=k_z=1 for simplicity. The ground-truth reference is obtained from the three-dimensional OS operator. To satisfy the incompressibility constraint, the initial velocity perturbation fields in the three-dimensional domain are prescribed as u(x,y,z) u(x,y,z) =ϵ⋅ℜ[ϕ(y)ei(kxx+kzz)], =ε· [φ(y)e^i(k_xx+k_zz) ], (24) v(x,y,z) v(x,y,z) =ϵ⋅ℜ[ikxϕ′(y)−ikzϕ(y)kx2+kz2ei(kxx+kzz)], =ε· [ ik_x\,φ (y)-ik_z\,φ(y)k_x^2+k_z^2e^i(k_xx+k_zz) ], w(x,y,z) w(x,y,z) =ϵ⋅ℜ[ikzϕ′(y)+ikxϕ(y)kx2+kz2ei(kxx+kzz)], =ε· [ ik_z\,φ (y)+ik_x\,φ(y)k_x^2+k_z^2e^i(k_xx+k_zz) ], where ϕ(y)φ(y) is the profile function introduced previously, and ϵε controls the magnitude of the initial perturbation. In this three-dimensional case, we only generate a strongly nonlinear dataset with ϵ=10−1ε=10^-1. Figure 11: Comparison of the first-order forcing and response modes of the streamwise velocity field at the ω=0.38ω=0.38 for the reduced 3D channel flow system. The original three-dimensional data correspond to a state vector of dimension close to 2×1052× 10^5. Here, the first 100 POD modes and their associated temporal coefficient trajectories are retained. Since the temporal coefficients form a one-dimensional vector, we employ a four-layer multilayer perceptron with 256 neurons in each hidden layer. The truncation level used for the projection eigensubspace in the resolvent analysis is chosen to be identical to the number of POD modes, and thirty-step temporal unrolling is used during training. Fig. 9 compares the eigenspectra and resolvent gain curves obtained from the OS operator with those predicted by the N-based method. Similar to the two-dimensional case, the dominant modes can still be categorized into near-wall A modes and center-region P modes, while the dominant resonant peak now shifts to ω=0.38ω=0.38. Fig. 10 and Fig. 11 further compare the dominant eigenmodes and resolvent modes, respectively. Owing to the additional spanwise dimension, however, the velocity structures are no longer simple two-dimensional wave packets propagating along the streamwise direction. Instead, they appear as a set of inclined wave-like structures traversing the flow field, namely, the oblique waves (Jovanović and Bamieh, 2005). Such oblique waves are known to play an important role because they can experience transient amplification through the coupling between wall-normal velocity and vorticity perturbations via a vortex-stretching mechanism. We notice that the resulting neural emulator trained on reduced POD coefficients still captures the leading linear stability characteristics of the original full-state system. In addition to the accurately identified first A and P eigenvalues, several higher‑order eigenvalues also fall along the branch of the analytical results. A similar phenomenon is observed in the reduced complex Ginzburg–Landau study in the SI Appendix, section 5, where a larger portion of the higher‑order modes is able to approach the analytical results. Although higher‑order modes correspond to rapidly decaying dynamics, POD separates modes into different energy levels beforehand. Rather than forcing the neural emulator to discover these less pronounced structures directly from the full state, learning the dynamics within the POD modal subspace is apparently more efficient. The accurate recovery of the leading eigenmodes naturally leads to excellent agreement in the corresponding input-output behavior in Fig. 11. Moreover, the extra peak at zero frequency becomes more pronounced, which is expected since a larger number of non‑oscillatory spurious modes appear compared with the two‑dimensional case. In addition to reducing the dimensionality of the learning problem, POD reduction offers several further advantages. Traditional U-net architectures can only enlarge the receptive field gradually through stacked encoder blocks, while POD modes already provide global spatial patterns of the dominant flow structures. The neural network thus only needs to learn the dynamics within the subspace spanned by these physically relevant modes, which is particularly important for fluid systems with global spatial interactions. Moreover, the extracted POD modes for channel flow usually inherit the symmetry or central symmetry with respect to the channel centerline. Therefore, the reconstructed stability and resolvent modes, expressed as combinations of POD modes, provide a powerful inductive bias to preserve the expected symmetry. 4 Discussion and Outlook In this work, we have developed a data-driven framework for extracting dominant local linear characteristics from a differentiable neural emulator. Rather than requiring explicit access to the governing equations or a tractable linearized operator, the method learns the nonlinear flow map directly from trajectory data and then recovers local Jacobian-based information through automatic differentiation. In this way, linear stability analysis and resolvent analysis can be extended to equation-free settings, where conventional operator-based approaches are difficult or unavailable. Its effectiveness is demonstrated on four canonical systems of increasing complexity and nonlinearity, underscoring its potential for analyzing complex dynamics. The fact that such results are achieved purely through the representational capacity of primitive network architectures is highly encouraging. We also compared our framework with the typical data-driven approach, DMD. The results show a clear regime dependence. When the training data are generated by linear or weakly nonlinear dynamics, DMD yields an accurate eigenspectrum reconstruction, which is expected because it is fundamentally designed to approximate linear evolution operators directly from data. As the data move further into a fully nonlinear regime, DMD loses even the dominant spectral structure, whereas our framework can still retain sufficient local dynamical information to recover the leading Jacobian-based characteristics. Although an exact one-to-one reconstruction of the full eigenspectrum is challenging, the dominant local linear characteristics are usually most important and physically relevant. As a novel data-driven approach, previous extensions and modifications for traditional methods such as DMD may offer valuable insights for improving our approach. One promising direction lies in incorporating sparsity‑promoting constraints to automatically select the few most dynamically relevant modes, in a manner analogous to sparsity‑promoting DMD(Jovanović et al., 2014). Moreover, since DMD can handle non-uniformly sampled datasets to circumvent the Nyquist criterion (Guéniat et al., 2015), our N emulator could also be extended to such data by adopting higher-order integration schemes rather than the direct mapping used in this work. As with other data-driven approaches, N-based modal analysis has the potential to significantly lower the barrier to entry for conducting operator-based modal analysis across diverse fields. Promising application domains include systems in which repeated trajectories are available, such as neural recordings around baseline-to-event transitions (Zhao et al., 2025) or climate reanalysis data around recurrent circulation regimes (Van der Avoird et al., 2002), although the assessment of our framework under noise, partial observability, and sampling limitations requires a dedicated future study. On the other hand, N-based modal analysis also provides a perspective for assessing the fidelity of a trained N emulator in capturing the system’s true dynamics, or for diagnosing prediction divergence through Jacobian‑based stability analysis, especially when reference modal results are available. This synergy between machine learning and classical modal analysis opens a promising avenue toward interpretable, physics-informed modeling and advances the integration of data-driven and theoretical approaches. Supporting Information Appendix A Modal analysis A.1 Linear stability analysis When there is no exogenous forcing applied to the original dynamical system, the governing equation of the perturbation ′ q degrades to a homogeneous system: d′dt=′, d q dt=A q , (25) where we can assume the perturbation in the form of normal mode ′=^eσt q = qe^σ t with σ=λ+iωσ=λ+iω. The real part λ and the imaginary part ω are the growth rate and characteristic frequency, respectively, and q describes the mode shape. By substituting the normal mode into Eq. 25, we finally obtain an eigenvalue problem as follows: ^=σ^.A q=σ q. (26) Since the Jacobian matrix A depends on the base state b q_b, it indicates that for a dynamic system with multiple fixed points, they may have different stability properties. It depends on whether each corresponding A has at least one eigenvalue with a positive real part λ. However, linear stability analysis provides only a partial description of the system dynamics because it solely focuses on unstable modes. In practice, a system can remain asymptotically stable with all eigenvalues residing in the stable region, yet still exhibit significant transient growth of perturbations. This motivates the use of such tools as resolvent analysis, in order to obtain a input-output viewpoint of the perturbation evolution. A.2 Resolvent analysis Since we aims to examine the response outputs of perturbation ′ q for different forcing inputs and identify the optimal input mode which yields the most amplified energy by internal dynamics, the external forcing should be taken into account. Thus, the evolution of the perturbation ′ q is described by a forced linear dynamical system as: d′dt=′+. d q dt=A q + f. (27) At a steady base state, we assume the response of perturbations ′ q and the external forcing f to be harmonic, which can be expressed as ′(t)=^e−iωt+c.c. q (t)= qe^-iω t+c.c. and (t)=^e−iωt+c.c. f(t)= fe^-iω t+c.c., where c.c.c.c. represents the complex conjugate and ω is the angular driving frequency. Based on this assumption, Eq. 27 can be rewritten as ^=(−iω−)−1^=(ω)^, q=(-i -A)^-1 f=H(ω) f, (28) where (ω)∈ℂN×NH(ω) ^N× N is defined as the resolvent operator, i.e., transfer function matrix between the forcing inputs and response outputs at frequency ω. Since the goal is to obtain the optimal energy amplification over all possible forcing modes q, a physically meaningful norm for energy measurement needs to be selected. A state variable norm can be defined as ‖^‖2=^∗^|| q||^2_Q= q^*Q q, where ()∗()^* denotes the Hermitian transpose and =∗Q=F^*F is a positive-definite weighting matrix that may account for grid discretization size (Yuan et al., 2023), compressible flow energy (Chu, 1965) and other factors. The Q norm is related to the standard Euclidean norm as ‖^‖2=‖^‖22|| q||^2_Q=||F q||^2_2. The largest input-output gain can be quantified via the ratio of their respective norms σ2(ω)=max^‖′^‖2‖^‖2=max^‖(ω)^‖2‖^‖2=‖(ω)−1‖22.σ^2(ω)= _ f || q ||^2_Q|| f||^2_Q= _ f ||H(ω) f||^2_Q|| f||^2_Q=||FH(ω)F^-1||^2_2. (29) Eq. 29 is an optimization problem, whose solution is given by the weighted SVD of the resolvent operator (ω)−1=(ω)(ω)∗(ω),FH(ω)F^-1= _F(ω) (ω) ^*_F(ω), (30) where ∈ℝn ^n is the singular value matrix whose diagonal entries are the optimal gains associated with the leading n resolvent modes, and the column vectors in −1(ω)F^-1 _F(ω) and −1(ω)F^-1 _F(ω) are the corresponding response and forcing modes. Appendix B Complex Ginzburg-Landau equation B.1 Supercritical and subcritical behaviors The Ginzburg–Landau equation can describe two fundamentally different behaviors of real flows, depending on the comparison between the bifurcation parameter μ0 _0 and the critical value for global instability μc _c (μc=0.4 _c=0.4 for our case used in the paper). The subcritical regime (0<μ0<μc0< _0< _c) is often referred to as a noise amplifier, in which the advection is so strong that all perturbations will be convected downstream before amplification. This corresponds to the convectively unstable flows, such as boundary layers (Ehrenstein and Gallaire, 2005) and mixing layers (Ho and Huerre, 1984). The supercritical regime (μ0>μc _0> _c) corresponds to self-sustained oscillatory behavior, where perturbations will grow exponentially until the nonlinear cubic term causes them to saturate and oscillate, as shown in Fig. S7(b). This can be used to model global unstable flows, such as hot jets (Lesshafft et al., 2006) and cylinder wakes (Lauga and Bewley, 2004; Cohen et al., 2005). B.2 Discrete representation To discretize the linearized operator A of the Ginzburg-Landau equation, we first need to create a domain discretization using Hermite functions (Bagheri et al., 2009). A total of N=220N=220 collocation points x1,x2,⋯,xNx_1,x_2,·s,x_N are given by the roots of Hn(χx)H_n(χ x), where HnH_n is the n-th Hermite polynomial and χ=(−μ2/(2γ))1/4χ=(- _2/(2γ))^1/4 is a scaling factor. With this choice of N, the discretized domain stretches from x1=−84.99x_1=-84.99 to xN=84.99x_N=84.99. Then the corresponding Hermite differentiation matrices provided by (Weideman and Reddy, 2000) are adopted to formulate the first derivative ∂/∂x∂/∂ x and second derivative ∂2/∂x2∂^2/∂ x^2 to construct the discrete matrix A of size N×N× N. Appendix C Hydrodynamic stability C.1 Linearized Navier-Stokes equations The motion of a two-dimensional incompressible flow is governed by the Navier-Stokes equations in nondimensional form (characteristic velocity U and length L are used to nondimensionalize) as: ∂t+⋅∇=−∇p+Re−1∇2,∇⋅=0. ∂ u∂ t+ u·∇ u=-∇ p+Re^-1∇^2 u,~~~~∇· u=0. (31) where u is the velocity field, p is the pressure and ReRe is the Reynolds number. In accordance with the conventional framework of linear stability analysis, we can linearize Eq. 31 around a steady base state (b,Pb)( U_b,P_b) to obtain the linearized Navier-Stokes equations governing perturbations (u′,p′)(u ,p ): ∂′∂t+b⋅∇′+∇b⋅′=−∇p′+Re−1∇2′,∇⋅′=0. ∂ u ∂ t+ U_b·∇ u +∇ U_b· u =-∇ p +Re^-1∇^2 u ,~~~~∇· u =0. (32) Based on aforementioned normal mode expansions (′(x,y,t),p′(x,y,t))T=(^(x,y),p^(x,y))Teσt( u (x,y,t),p (x,y,t))^T=( u(x,y), p(x,y))^Te^σ t with σ=λ+iωσ=λ+iω, we substitute it into the (32), which then can be reformulated in matrix form: (^p^)=σ(^p^),with=(−b⋅∇−∇b⋅+Re−1∇2−∇⋅).A pmatrix u\\ p pmatrix=σ pmatrix u\\ p pmatrix, ~~A= pmatrix- U_b·∇-∇ U_b·+Re^-1∇^2&-∇\\ ∇·&0\\ pmatrix. (33) This leads to the eigenvalue problem about matrix A, where the stability of the flow is determined by the sign of the imaginary part of the eigenvalues σ. And the non-trivial solutions of (^,q^)T( u, q)^T constitute global linear modes of the problem. For complex flow configurations (e.g., cylinder flow), it is rarely practical to assemble and store the full matrix representation A explicitly. This global eigenvalue problem can be solved by an iterative approach, and the most popular one is the Arnoldi algorithm (Arnoldi, 1951; Saad, 1980). It is a time-stepping–based Jacobian-free method and has been widely used for linear stability analyses in the fluid mechanics community. C.2 Orr-Sommerfeld/Squire equations For comparatively simple parallel flows, one may manipulate the linearized Navier-Stokes equations ((32)) analytically and solve smaller eigenproblems. For instance, as the base flow of the laminar channel flow case is known as U(y)=1−y2U(y)=1-y^2, we can plug it into (32) and obtain its component form: ∂u′∂t+U∂u′∂x+v′dUdy ∂ u ∂ t+U ∂ u ∂ x+v dUdy =−∂p′∂x+Re−1∇2u′ =- ∂ p ∂ x+Re^-1∇^2u (34) ∂v′∂t+U∂v′∂x ∂ v ∂ t+U ∂ v ∂ x =−∂p′∂y+Re−1∇2v′ =- ∂ p ∂ y+Re^-1∇^2v (35) ∂w′∂t+U∂w′∂x ∂ w ∂ t+U ∂ w ∂ x =−∂p′∂z+Re−1∇2w′ =- ∂ p ∂ z+Re^-1∇^2w (36) ∂u′∂x+∂v′∂y+∂w′∂z ∂ u ∂ x+ ∂ v ∂ y+ ∂ w ∂ z =0 =0 (37) We introduce the wall-normal velocity v′(x,y,z,t)v (x,y,z,t) and the wall-normal vorticity η′(x,y,z,t)=∂u′∂z−∂w′∂xη (x,y,z,t)= ∂ u ∂ z- ∂ w ∂ x. By combining (34)-(37) to eliminate the pressure and tangential velocity components, one obtains a coupled system for v′v and η′η . The equation governing v′v is the three-dimensional Orr-Sommerfeld equation, ∂t(∇2v′)+U∂x(∇2v′)−d2Udy2∂v′∂x=Re−1∇4v′, ∂ t(∇^2v )+U ∂ x(∇^2v )- d^2Udy^2 ∂ v ∂ x=Re^-1∇^4v , (38) while the wall-normal vorticity satisfies the Squire equation, ∂η′∂t+U∂η′∂x+dUdy∂v′∂z=Re−1∇2η′. ∂η ∂ t+U ∂η ∂ x+ dUdy ∂ v ∂ z=Re^-1∇^2η . (39) Here, ∇2=∂xx+∂yy+∂zz∇^2= _x+ _y+ _z and ∇4=∇2∇2∇^4=∇^2∇^2. We can further assume perturbations in normal modes exhibit both streamwise and spanwise periodic dependence, i.e. v′,η′(x,y,z,t)=v^(y),η^(y)ei(kxx+kzz)eσt=v^(y),η^(y)eikx(x−ct)eikzz,\v ,η \(x,y,z,t)=\ v(y), η(y)\e^i(k_xx+k_zz)e^σ t=\ v(y), η(y)\e^ik_x(x-ct)e^ik_zz, (40) where v^(y) v(y) and η^(y) η(y) are the wall-normal amplitude functions, kxk_x and kzk_z are the real streamwise and spanwise wavenumbers, and c=cr+icic=c_r+ic_i is the complex wave speed. In other words, the perturbation is considered in wavenumber space to facilitate simplification. Substituting this normal mode form into (38) and (39) transforms the partial derivatives as follows: ∂t→−ikxc,∂x→ikx,∂z→ikz,∇2→D2−k2, ∂ t→-ik_xc, ∂ x→ ik_x, ∂ z→ ik_z, ∇^2→ D^2-k^2, where D=dyD= ddy and k2=kx2+kz2k^2=k_x^2+k_z^2. After simplification, we arrive at the Orr-Sommerfeld/Squire system: (U−c)(D2−k2)v^−d2Udy2v (U-c)(D^2-k^2) v- d^2Udy^2 v =1ikxRe(D2−k2)2v^, = 1ik_xRe(D^2-k^2)^2 v, (41) (U−c)η^+kzkxdUdyv (U-c) η+ k_zk_x dUdy v =1ikxRe(D2−k2)η^. = 1ik_xRe(D^2-k^2) η. (42) Since the no-slip boundary conditions at the walls (y=±1y=± 1) require u′=v′=w′=0u =v =w =0, the corresponding boundary conditions become v^(±1)=0,Dv^(±1)=0,η^(±1)=0. v(± 1)=0, D v(± 1)=0, η(± 1)=0. The Orr-Sommerfeld/Squire equations together with these homogeneous boundary conditions, again form an eigenvalue problem. It is worth noting that the two-dimensional case is recovered by setting kz=0k_z=0, in which case the Squire equation decouples and the above system reduces to the classical Orr-Sommerfeld eigenvalue problem. Appendix D Ablation study In this section, we present a series of ablation and robustness studies based on the Lorenz system to further improve our understanding of neural emulators. Before doing so, we first introduce the metrics used to evaluate emulator performance. To quantify the emulator prediction accuracy, we adopt the metrics introduced in our previous work(Koehler et al., 2024), namely the normalized root mean squared error (nRMSE) and its temporally aggregated version. The mean nRMSE over M samples (30 by default) is defined as LnRMSE=1M∑i=1M‖^i−i‖2‖i‖2.L_nRMSE= 1M _i=1^M \| q_i- q_i \|^2 \| q_i \|^2. (43) This metric accounts for scale differences between the predicted state q and the reference state q, which is especially useful when evaluating rollout errors for dynamical systems whose state amplitudes vary substantially in time, such as decaying channel flows. Since we are particularly interested in how the error develops over time, we further aggregate the nRMSE over 100 rollout steps using the geometric mean (Agg. nRMSE) LAgg. nRMSE=exp(1100∑t=1100log(LnRMSE[t])),L_Agg. nRMSE= ( 1100 _t=1^100 (L_nRMSE^[t]) ), (44) where the superscript [t][t] denotes the nRMSE after t time steps. For the Lorenz system, whose Lyapunov time is of order one time unit, this 100-step aggregation provides a meaningful measure of temporal generalization. The geometric mean removes the requirement to hand-tune upper limits for temporal aggregation in cases where error metrics surpass 1. The relative Frobenius-norm error is used to quantify the prediction accuracy of element-wise Jacobian evaluated at state q as follows: ε()=‖N()−ana()‖F‖ana()‖F. ( q)= \|J_N( q)-J_ana( q) \|_F \|J_ana( q) \|_F. (45) D.1 Dataset range First, we ablate both the number of trajectories and the length of the training temporal horizon. In this study, we kept the network architecture identical to that used in the main text and evaluated combinations of temporal horizons of 100, 400, and 800 time steps together with training sets containing 1, 2, 3, 5, 10, 15, 20, 30, 40, 50, 70, and 100 trajectories. For each combination, we trained five neural networks and plotted a shaded line curve showing the mean aggregated nRMSE as a function of the number of dataset samples, with the shaded region indicating the corresponding standard deviation. The results are presented in Fig. S1. Across all combinations of these two parameters, the performance improves consistently as the amount of training data increases, up to a certain number of training samples beyond which gains become minor. This indicates that the emulator’s temporal generalization capability can converge with substantially less data than that used in the main text. The 80,000-snapshot dataset was chosen for robustness rather than as a minimal requirement. Owing to the chaotic nature of the Lorenz system, extending the temporal horizon effectively provides additional independent samples, which is reflected by the convergence of the three curves to a similar performance level. However, this does not suggest that the temporal horizon plays an insignificant role. For systems with multiple dynamical stages, an excessively short horizon may fail to include important physical regimes. For example, in the Burgers equation, a horizon that is too short may miss the shock propagation stage. This is relevant because data-driven emulators may struggle to recover such behavior unless it is sufficiently represented in the training data. More generally, the amount of data required by our approach is not a universal number, but depends on several factors: (i) the dimensionality of the system, (i) the degree of nonlinearity, (i) whether the system has multiple stages, and (iv) how well the data cover the neighborhood of the base state at which the Jacobian is queried. D.2 Dataset regime To assess how the learned neural emulator generalizes beyond the region covered by the training data, we perform an ablation study by deliberately restricting the training trajectories to a single lobe of the attractor. The purpose of this experiment is to examine how the accuracy of the surrogate-based Jacobian depends on the local support of the training data in state space. Fig. S2(a) illustrates the data construction. The 100 initial conditions denoted by blue circles are used to generate the training data, and the red curves denote the corresponding trajectories used for training. Each training trajectory is evolved for 88 time units so that the total number of samples is the same as in the setting of the main text. To restrict the training data to the right lobe, a trajectory is accepted into the training set only if the proportion of points satisfying x>0x>0 is less than 1%1\%. Using this criterion, the associated trajectories all concentrated on a single attractor lobe. The gray background trajectory in Fig. S2(a) is a long reference trajectory generated from three initial conditions located near the three Lorenz equilibria. Each trajectory is integrated with a time step Δt=0.01 t=0.01 for a total duration of 5050 time units, so that the resulting reference set spans both lobes of the attractor as well as the transitions between them. After training, the neural emulator is evaluated via the relative Frobenius-norm error on these long reference trajectories. We again train five neural networks with the same architecture as before. The color shown in Fig. S2(b) corresponds to the mean relative Frobenius error evaluated along the reference trajectory. The results reveal a clear spatial dependence of the Jacobian accuracy. In the region of state space well represented by the training trajectories, the neural emulator recovers the Jacobian with substantially smaller error. As the evaluation points move away from that region, particularly toward the opposite lobe and the connecting transition region, the error increases progressively. This behavior indicates that the learned surrogate does not provide uniformly accurate local linearizations throughout the full attractor when its training support is restricted to only one portion of the state space. This experiment highlights an important aspect of the proposed framework. Although the neural emulator can exhibit a certain degree of extrapolative capability, its Jacobian is most reliable in regions that are sufficiently sampled during training. In other words, the method should be interpreted as a tool for extracting local linear characteristics within a well-covered dynamical regime, rather than as a model expected to maintain quantitative Jacobian accuracy across dynamically distinct regimes that are absent from the training data. D.3 Noise level In this section, we performed an additional ablation study to assess the noise robustness of the neural emulator. This experiment probes the extent to which the learned one-step map and its induced local linearization remain reliable when the training trajectories are contaminated by measurement-like perturbations. We perturbed the training data with additive random noise with noise levels of 1%1\%, 2%2\%, 5%5\%, and 10%10\% and then constructed the corresponding one-step state-transition pairs from the noisy data. For each noise level, to evaluate the fidelity of the learned local dynamics, we compared the eigenvalues of the Jacobians predicted by the neural emulator at the three Lorenz equilibria and one non-equilibrium state against the analytical reference values; the results are summarized in Tab. S3. The results show that the emulator remains remarkably robust to moderate noise contamination in the training data. As the noise level increases from 1%1\% to 10%10\%, the predicted eigenvalues indeed deviate from the ground truth progressively, indicating the expected degradation in local linear characteristic accuracy. However, this deviation is acceptable within low to moderate noise levels. This suggests that, for the present Lorenz setting, the emulator is able to retain a stable approximation of the underlying dynamics even when the training observations are corrupted to some extent. Moreover, we can also notice that the spectral characteristics of two fixed points associated with the attractor centers are consistently more difficult to predict than the fixed point at the origin and the non-equilibrium state. This increased difficulty is likely related to the strongly chaotic dynamics near the attractor centers in the Lorenz system, where trajectories initiated from nearby points can diverge at different rates and along different directions. Overall, this added noise study shows that the proposed framework is not purely limited to noise-free benchmark data for the Lorenz case. Appendix E Additional examples E.1 2D cylinder flow As a unsteady example for our framework, we consider a flow with an immersed obstacle, the two-dimensional cylinder flow at Re = 100 (based on the cylinder diameter and inflow velocity). This Reynolds number setting is greater than the critical threshold, and the cylinder wake flow experiences a Hopf bifurcation and evolves to be a time-periodic non-symmetric state. In contrast to channel flow, the linearized operator for such an intricate configuration cannot be derived analytically and must instead be computed numerically. The spectral element method (Fischer et al., 2007; Loiseau et al., 2014) is used for steady-state solution, numerical linear stability analysis and fully nonlinear simulations. The cylinder with a diameter D is placed at (0,0)(0,0) in a grid extending from −15D-15D to 35D35D and −15D-15D to 15D15D in x and y directions, respectively. We employ the same boundary conditions as in (Barkley, 2006) for linearized and nonlinear NS equations, respectively. The nonlinear simulations are initialized with the superposition of the base flow and different infinitesimal perturbations for faster saturation. The dataset comprises 30 trajectories, each sampled for m=200m=200 snapshots every 0.30.3 time units. A third-order temporal scheme with a time step of 0.0050.005 is chosen for nonlinear numerical simulation. The flow snapshots are uniformly resampled within [−8,12]×[−3,3][-8,12]×[-3,3] using an 80×2480× 24 Cartesian grid. A masking approach is applied by assigning zero values to the network inputs and outputs inside the cylinder. The network architecture is similar to that of the 2D channel flow case, except that the input tensor now comprises three features: pressure, streamwise velocity, and transverse velocity (see Fig. S9). The projection subspace is truncated to r=30r=30 modes, and the model is trained using a twenty-step unrolling strategy to ensure coverage of a full period of vortex shedding. The eigenspectra obtained from numerical direct stability analysis and N-based Jacobian are presented in Fig. S3(a), along with the leading eigenmodes shown in Fig. S3(b). As expected, the cylinder flow at Re=100Re=100 exhibits one unstable leading eigenvalue located in the upper half of the complex plane. The leading eigenvalue predicted by the N emulator at frequency f=0.12f=0.12 aligns well with the numerical ground truth, and its eigenmode in the downstream area of the cylinder also shows an excellent agreement. This observation again demonstrates that the N emulator can not only advance the flow field in an autoregressive fashion, but also capture spatial perturbation structures consistent with those in the numerical modal analysis, thereby reflecting the underlying dynamics of the system. Note that the frequency of the leading mode differs from the vortex shedding frequency, i.e., Strouhal number (0.1660.166) at Re=100Re=100. Previous literature (Barkley, 2006) has reported that only linear stability analysis around the mean flow can predict the vortex shedding frequency with good accuracy, but we reiterate that our current method is not applicable to time-averaged states, as this would invalidate the derivation. Since the steady-state solution of cylinder flow at Re=100Re=100 represents an artificial pseudo-equilibrium that is not physically realizable, any perturbation introduced around this base state will grow exponentially, ultimately leading to vortex shedding. Consequently, there is no reference for resolvent analysis based on this steady-state solution. In prior resolvent studies of cylinder wakes, the time-averaged flow is typically adopted as the base state, with the unknown nonlinear terms treated as effective forcing to uncover coherent structures. Nevertheless, we can follow the procedure of previous examples to conduct the resolvent analysis around the steady-state flow based on the obtained N-based Jacobian. The resulting gain curves and resolvent pressure modes at the optimal forcing frequency are presented in Fig. S4. Despite the lack of a ground truth solution, the results offer some clear physical insights: the optimal forcing frequency again coincides with the resonance frequency in Fig. S3, and the corresponding response mode forms symmetric lobe-like patterns downstream of the cylinder, similar to the leading eigenmode. The forcing mode also displays lobe-like structures in the downstream region much closer to the cylinder. This is likely associated with the wavemaker region, the area immediately downstream of the cylinder that is most sensitive to external forcing (Giannetti and Luchini, 2007; Marquet et al., 2008). E.2 Reduced complex Ginzburg-Landau system To further verify the spectral consistency between the POD-reduced system and the original physical-space system, we additionally consider a POD-based reduced-order experiment for the complex Ginzburg-Landau system. The sample trajectories and all parameter settings are kept identical to those of the full-order case, except that here we consider only the nonlinear dataset. Since the original spatial discretization yields a state vector of dimension 220, we retain only the first 20 POD modes in the reduced representation. A four-layer complex-valued perceptron with 256 neurons in each hidden layer is employed, and only a ten-step rollout is used during training. As shown in Fig. S5, we compare the heatmaps of the N-based Jacobian learned from the reduced data with the analytical Jacobian NPODN_POD obtained by projecting the full Jacobian onto the POD subspace. The improvement is particularly evident from the prediction difference in Fig. S5, where the error magnitude is significantly smaller than that of the corresponding full-order case shown in the main text. This improvement is expected since POD substantially reduces the input and output dimensions of the neural emulator, thereby simplifying the learning problem. As we mentioned in the section “Stability in reduced order models”, if one is only interested in the leading eigenvalues, it is not necessary to retain all POD modes to construct this subspace. In Fig. S6(a), we show the consistency between the eigenvalues of the full Jacobian N and those of PODN_POD when different numbers of POD modes are retained. It can be seen that retaining at most 100 POD modes is already sufficient to recover the first seven eigenvalues with high accuracy. Even with only 20 POD modes, the reduced system still captures the unique leading unstable mode of the system. Figure S6(b) further shows the eigenvalues predicted by the neural emulator trained on the reduced POD coefficients. In addition to the leading unstable mode, many higher-order modes are also found to approach the analytical spectrum. This observation is consistent with the very small Jacobian error seen in Fig. S5. By contrast, the N-based eigenvalues obtained in the main text from training on the full-order data show much poorer accuracy for the higher-order modes (Fig. 4(b) in the main text). This behavior is physically reasonable for the complex Ginzburg-Landau system. For most snapshots, the state amplitude is concentrated only in a localized region near the origin, while the far-field state remains close to zero. Consequently, when training directly on the full-order data, a large fraction of the input-output mapping learned by the neural network is spent on dynamically uninformative regions of the domain. By projecting the state onto a compact subspace spanned by the dominant POD modes, the neural emulator is relieved from learning spatially redundant information and can focus on the temporal‑coefficient dynamics associated with structures at different energy levels. As a consequence, the resulting reduced-order model yields a more accurate and significantly improved prediction of the Jacobian compared with direct full-state learning. Figure S1: Effect of training trajectory number and temporal horizon on emulator generalization in the Lorenz system. The blue line indicates the aggregated nRMSE error across five trained networks, while the shaded region denotes the standard deviation. Figure S2: Effect of training regimes on emulator generalization in the Lorenz system. (a) Phase-space visualization of the Lorenz attractor, with the training trajectories (red curves) confined to one lobe and initial points marked (blue circles). The full attractors (gray curves) are generated by integrating trajectories initialized near the three equilibria. (b) Frobenius norm error between N-based and operator-based Jacobians along the full attractors. Base state b q_b Eigenvalues 1%1\% 2%2\% 5%5\% 10%10\% Ground Truth [000] bmatrix0\\ 0\\ 0 bmatrix [1.12570.97350.7939] bmatrix1.1257\\ 0.9735\\ 0.7939 bmatrix [1.12560.97350.7979] bmatrix1.1256\\ 0.9735\\ 0.7979 bmatrix [1.12540.98110.7600] bmatrix1.1254\\ 0.9811\\ 0.7600 bmatrix [1.12431.01350.6988] bmatrix1.1243\\ 1.0135\\ 0.6988 bmatrix [1.12560.97370.7959] bmatrix1.1256\\ 0.9737\\ 0.7959 bmatrix [626227] bmatrix6 2\\ 6 2\\ 27 bmatrix [0.86850.9957±0.1018i] bmatrix0.8685\\ 0.9957± 0.1018i bmatrix [0.85900.9947±0.1018i] bmatrix0.8590\\ 0.9947± 0.1018i bmatrix [0.81900.9910±0.0997i] bmatrix0.8190\\ 0.9910± 0.0997i bmatrix [0.62060.9553±0.0958i] bmatrix0.6206\\ 0.9553± 0.0958i bmatrix [0.87060.9957±0.1019i] bmatrix0.8706\\ 0.9957± 0.1019i bmatrix [−62−6227] bmatrix-6 2\\ -6 2\\ 27 bmatrix [0.86780.9958±0.1014i] bmatrix0.8678\\ 0.9958± 0.1014i bmatrix [0.86150.9948±0.1014i] bmatrix0.8615\\ 0.9948± 0.1014i bmatrix [0.81900.9877±0.0994i] bmatrix0.8190\\ 0.9877± 0.0994i bmatrix [0.63130.9481±0.0868i] bmatrix0.6313\\ 0.9481± 0.0868i bmatrix [0.87060.9957±0.1019i] bmatrix0.8706\\ 0.9957± 0.1019i bmatrix [1.46898.7045−7.6905] bmatrix1.4689\\ 8.7045\\ -7.6905 bmatrix [1.14430.97930.7775] bmatrix1.1443\\ 0.9793\\ 0.7775 bmatrix [1.14740.95710.7881] bmatrix1.1474\\ 0.9571\\ 0.7881 bmatrix [1.14520.98200.7661] bmatrix1.1452\\ 0.9820\\ 0.7661 bmatrix [1.16090.93410.7780] bmatrix1.1609\\ 0.9341\\ 0.7780 bmatrix [1.14610.97750.7786] bmatrix1.1461\\ 0.9775\\ 0.7786 bmatrix Table S3: Comparison of Jacobian eigenspectra obtained from neural emulators trained on the Lorenz dataset with different levels of noise: 1%1\%, 2%2\%, 5%5\%, and 10%10\%. Figure S3: Linear stability analysis results of the 2D cylinder flow system. (a) Eigenspectra of the numerical direct stability analysis (operator-based) and N-based Jacobian for unsteady two-dimensional cylinder flow. (b) Eigenmodes of the unstable eigenvalue from the operator-based and N-based Jacobians. From top to bottom: pressure, streamwise velocity, and transverse velocity. Figure S4: Resolvent analysis results of the 2D cylinder flow system. (a) Resolvent gain as a function of frequency in the interval [0, 3]. (b) Optimal pressure forcing and response modes evaluated at the frequency corresponding to the peak gain. Figure S5: Comparison of the real and imaginary parts of the reduced Jacobians between the N-based operator and the operator-based ground truth. Figure S6: Linear stability analysis results of the POD-reduced complex Ginzburg-Landau systems. (a) Operator-based eigenvalues obtained by projecting the full Jacobian onto POD subspaces with different truncation levels, compared with the full-order spectrum. (b) Comparison between the operator-based and N-based eigenspectra in the 20-mode POD subspace. Figure S7: Sample trajectories from datasets of (a) linear and (b) nonlinear complex Ginzburg-Landau systems. Figure S8: Snapshots of the wall-normal velocity perturbation field during the transient evolution stage of two-dimensional channel flow. The initial perturbation profiles are the same but with different magnitude scaling factors ϵε: (a) weakly nonlinear scenario with ϵ=10−6ε=10^-6; (b) nonlinear scenario with ϵ=10−1ε=10^-1. The decaying perturbations in (a) remain symmetric structures when evolving, while those in (b) encounter symmetry-breaking phenomena. Figure S9: Neural emulator architecture for all examples used in this paper. (a) A three-layer perceptron with heterogeneous activation functions (log-sigmoid, radial basis, and linear), used for the Lorenz system. (b) A six-layer complex-valued perceptron with complex GELU activation function (Hendrycks and Gimpel, 2016) for each layer, used for the full-order complex Ginzburg-Landau system. (c) A classical U-net used for the 2D channel flow syste,, where the two-channel inputs and outputs represent the perturbation fields of streamwise and transverse velocity around the base flow, respectively. (d) A four-layer perceptron used for the POD-reduced 2D channel flow system. (e) A classical U-net used for the two-dimensional cylinder flow system, where the three-channel inputs and outputs represent the pressure, streamwise velocity, and transverse velocity fields, respectively. (d) A four-layer complex-valued perceptron used for the POD-reduced complex Ginzburg-Landau system system. Data Availability The source code for the current study is available at Github (https://github.com/tum-pbs/NonlinearRA). Data available upon request from chengyun.wang@tum.de. Previously published platform (Koehler et al., 2024) used for this work is available at GitHub (https://github.com/tum-pbs/apebench). Acknowledgements C.W. is supported by the China Scholarship Council (No. 202406340048) and L.C. is supported by the National Natural Science Foundation of China (Grant No. 92470120). We are grateful for many fruitful discussions with Jean-Christophe Loiseau and Eduardo Martini. We would also like to thank Benjamin Herrmann for kindly providing the dataset used in his paper. Competing Interests The authors declare no competing interests. References [1] K. Allen, T. Lopez-Guevara, K. L. Stachenfeld, A. Sanchez Gonzalez, P. Battaglia, J. B. Hamrick, and T. Pfaff (2022) Inverse design for fluid-structure interactions using graph network simulators. Advances in Neural Information Processing Systems 35, p. 13759–13774. Cited by: §3. [2] W. E. Arnoldi (1951) The principle of minimized iterations in the solution of the matrix eigenvalue problem. Quarterly of applied mathematics 9 (1), p. 17–29. Cited by: §C.1. [3] P. J. Baddoo, B. Herrmann, B. J. McKeon, and S. L. Brunton (2022) Kernel learning for robust dynamic mode decomposition: linear and nonlinear disambiguation optimization. Proceedings of the Royal Society A: Mathematical, Physical and Engineering Sciences 478 (2260). Cited by: §1. [4] P. J. Baddoo, B. Herrmann, B. J. McKeon, J. Nathan Kutz, and S. L. Brunton (2023) Physics-informed dynamic mode decomposition. Proceedings of the Royal Society A 479 (2271), p. 20220576. Cited by: §1. [5] S. Bagheri, D. Henningson, J. Hœpffner, and P. J. Schmid (2009) Input-output analysis and control design applied to a linear model of spatially developing flows. Applied Mechanics Reviews 62 (2), p. 020803. Cited by: §B.2, §3, §3, §3. [6] Y. Bar-Sinai, S. Hoyer, J. Hickey, and M. P. Brenner (2019) Learning data-driven discretizations for partial differential equations. Proceedings of the National Academy of Sciences 116 (31), p. 15344–15349. Cited by: §1. [7] D. Barkley and R. D. Henderson (1996) Three-dimensional floquet stability analysis of the wake of a circular cylinder. Journal of Fluid Mechanics 322, p. 215–241. Cited by: §2. [8] D. Barkley (2006) Linear analysis of the cylinder wake mean flow. Europhysics Letters 75 (5), p. 750. Cited by: §E.1, §E.1. [9] S. L. Brunton and J. N. Kutz (2022) Data-driven science and engineering: machine learning, dynamical systems, and control. Cambridge University Press. Cited by: §3, §3. [10] S. L. Brunton, J. L. Proctor, and J. N. Kutz (2016) Discovering governing equations from data by sparse identification of nonlinear dynamical systems. Proceedings of the national academy of sciences 113 (15), p. 3932–3937. Cited by: §1. [11] K. K. Chen and C. W. Rowley (2011) H2 optimal actuator and sensor placement in the linearised complex ginzburg–landau system. Journal of Fluid Mechanics 681, p. 241–260. Cited by: §3. [12] L. Chen and N. Thuerey (2023) Towards high-accuracy deep learning inference of compressible flows over aerofoils. Computers & Fluids 250, p. 105707. Cited by: §1. [13] L. Chen and N. Thuerey (2024) Deep learning-based predictive modeling of transonic flow over an airfoil. Physics of Fluids 36 (12). Cited by: §1. [14] B. Chu (1965) On the energy transfer to small disturbances in fluid flow (part i). Acta Mechanica 1 (3), p. 215–234. Cited by: §A.2. [15] K. Chua, R. Calandra, R. McAllister, and S. Levine (2018) Deep reinforcement learning in a handful of trials using probabilistic dynamics models. Advances in neural information processing systems 31. Cited by: §3. [16] K. Cohen, S. Siegel, T. McLaughlin, E. Gillies, and J. Myatt (2005) Closed-loop approaches to control of a wake flow modeled by the ginzburg–landau equation. Computers & Fluids 34 (8), p. 927–949. Cited by: §B.1. [17] A. M. Edstrand, T. B. Davis, P. J. Schmid, K. Taira, and L. N. Cattafesta I (2016) On the mechanism of trailing vortex wandering. Journal of Fluid Mechanics 801, p. R1. Cited by: §2. [18] U. Ehrenstein and F. Gallaire (2005) On two-dimensional temporal modes in spatially evolving open flows: the flat-plate boundary layer. Journal of Fluid Mechanics 536, p. 209–218. Cited by: §B.1. [19] P. Fischer, J. Lottes, and H. Tufo (2007-05) Nek5000. Note: Nek5000 is high-order spectral element code used for the simulation of incompressible fluid flow, MHD, and heat transfer in general three-dimensional domains that has been specifically developed for deployment on large-scale distributed memory parallel architectures. The code features state of the art matrix free multilevel iterative solvers coupled with stable second and third order operator-split stepping schemes. Scalable performance has been established onup to 32,000 processors and the overall quality of the algorithmic approach, implementation, and performance were recognized with the award of the Gordon Bell Prize in 1999. External Links: Document, Link Cited by: §E.1. [20] D. Floryan and M. D. Graham (2022) Data-driven discovery of intrinsic dynamics. Nature Machine Intelligence 4 (12), p. 1113–1120. Cited by: §1. [21] F. Giannetti and P. Luchini (2007) Structural sensitivity of the first instability of the cylinder wake. Journal of Fluid Mechanics 581, p. 167–197. Cited by: §E.1. [22] J. Guckenheimer and P. Holmes (2013) Nonlinear oscillations, dynamical systems, and bifurcations of vector fields. Vol. 42, Springer Science & Business Media. Cited by: §1. [23] F. Guéniat, L. Mathelin, and L. R. Pastur (2015) A dynamic mode decomposition approach for large and arbitrarily sampled systems. Physics of Fluids 27 (2). Cited by: §4. [24] D. Hendrycks and K. Gimpel (2016) Gaussian error linear units (gelus). arXiv preprint arXiv:1606.08415. Cited by: Figure S9. [25] C. G. Hernández, K. Cao, B. Herrmann, S. L. Brunton, and B. J. McKeon (2023) Toward data-driven resolvent analysis of nonlinear flows. CTR Annual Research Briefs, p. 33–42. Cited by: §1. [26] B. Herrmann, P. J. Baddoo, R. Semaan, S. L. Brunton, and B. J. McKeon (2021) Data-driven resolvent analysis. Journal of Fluid Mechanics 918, p. A10. Cited by: §1, §2, §2, §3, §3, §3. [27] C. Ho and P. Huerre (1984) Perturbed free shear layers. Annual review of fluid mechanics 16, p. 365–424. Cited by: §B.1. [28] P. Holmes (2012) Turbulence, coherent structures, dynamical systems and symmetry. Cambridge university press. Cited by: §1. [29] M. R. Jovanović and B. Bamieh (2005) Componentwise energy amplification in channel flows. Journal of Fluid Mechanics 534, p. 145–183. Cited by: §3. [30] M. R. Jovanović, P. J. Schmid, and J. W. Nichols (2014) Sparsity-promoting dynamic mode decomposition. Physics of Fluids 26 (2). Cited by: §4. [31] M. R. Jovanović (2021) From bypass transition to flow control and data-driven turbulence modeling: an input–output viewpoint. Annual Review of Fluid Mechanics 53 (1), p. 311–345. Cited by: §1. [32] D. Kochkov, J. A. Smith, A. Alieva, Q. Wang, M. P. Brenner, and S. Hoyer (2021) Machine learning–accelerated computational fluid dynamics. Proceedings of the National Academy of Sciences 118 (21), p. e2101784118. Cited by: §1. [33] F. Koehler, S. Niedermayr, N. Thuerey, et al. (2024) APEBench: a benchmark for autoregressive neural emulators of pdes. Advances in Neural Information Processing Systems 37, p. 120252–120310. Cited by: Appendix D, Data Availability, §1, §3. [34] J. N. Kutz, S. L. Brunton, B. W. Brunton, and J. L. Proctor (2016) Dynamic mode decomposition: data-driven modeling of complex systems. SIAM. Cited by: §1. [35] E. Lauga and T. R. Bewley (2004) Performance of a linear robust control strategy on a nonlinear model of spatially developing flows. Journal of Fluid Mechanics 512, p. 343–374. Cited by: §B.1. [36] L. Lesshafft, P. Huerre, P. Sagaut, and M. Terracol (2006) Nonlinear global modes in hot jets. Journal of Fluid Mechanics 554, p. 393–409. Cited by: §B.1. [37] Z. Li, N. Kovachki, K. Azizzadenesheli, B. Liu, K. Bhattacharya, A. Stuart, and A. Anandkumar (2020) Fourier neural operator for parametric partial differential equations. arXiv preprint arXiv:2010.08895. Cited by: §1. [38] B. List, L. Chen, and N. Thuerey (2022) Learned turbulence modelling with differentiable fluid solvers: physics-based loss functions and optimisation horizons. Journal of Fluid Mechanics 949, p. A25. Cited by: §1. [39] J.-Ch. Loiseau, J.-Ch. Robinet, S. Cherubini, and E. Leriche (2014) Investigation of the roughness-induced transition: global stability analyses and direct numerical simulations. J. Fluid Mech. 760, p. 175–211. Cited by: §E.1. [40] J. L. Lumley (1967) The structure of inhomogeneous turbulent flows. Atmospheric turbulence and radio wave propagation, p. 166–178. Cited by: §1. [41] L. M. Mack (1976) A numerical study of the temporal eigenvalue spectrum of the blasius boundary layer. Journal of Fluid Mechanics 73 (3), p. 497–520. Cited by: §3. [42] O. Marquet, D. Sipp, and L. Jacquin (2008) Sensitivity analysis and passive control of cylinder flow. Journal of Fluid Mechanics 615, p. 221–252. Cited by: §E.1. [43] E. Martini, D. Rodríguez, A. Towne, and A. V. Cavalieri (2021) Efficient computation of global resolvent modes. Journal of Fluid Mechanics 919, p. A3. Cited by: §3. [44] M. McCabe, P. Harrington, S. Subramanian, and J. Brown (2023) Towards stability of autoregressive neural operators. Transactions on Machine Learning Research. Note: External Links: ISSN 2835-8856, Link Cited by: §1. [45] B. J. McKeon and A. S. Sharma (2010) A critical-layer framework for turbulent pipe flow. Journal of Fluid Mechanics 658, p. 336–382. Cited by: §1. [46] M. Mortensen (2017) Shenfun - automating the spectral galerkin method. In MekIT’17 - Ninth national conference on Computational Mechanics, B. H. Skallerud and H. I. Andersson (Eds.), p. 273–298. External Links: ISBN 978-84-947311-1-2, 1708.03188 Cited by: §3. [47] M. Mortensen (2018) Shenfun: high performance spectral galerkin computing platform. Journal of Open Source Software 3 (31), p. 1071. External Links: Document Cited by: §3. [48] J. Morton, A. Jameson, M. J. Kochenderfer, and F. Witherden (2018) Deep dynamical modeling and control of unsteady fluid flows. In Advances in Neural Information Processing Systems, Cited by: §1. [49] J. Page, P. Norgaard, M. P. Brenner, and R. R. Kerswell (2024) Recurrent flow patterns as a basis for two-dimensional turbulence: predicting statistics from structures. Proceedings of the National Academy of Sciences 121 (23), p. e2320007121. Cited by: §1. [50] J. L. Proctor, S. L. Brunton, and J. N. Kutz (2016) Dynamic mode decomposition with control. SIAM Journal on Applied Dynamical Systems 15 (1), p. 142–161. Cited by: §1. [51] M. Raissi, P. Perdikaris, and G. E. Karniadakis (2019) Physics-informed neural networks: a deep learning framework for solving forward and inverse problems involving nonlinear partial differential equations. Journal of Computational physics 378, p. 686–707. Cited by: §1. [52] S. C. Reddy and D. S. Henningson (1993) Energy growth in viscous channel flows. Journal of Fluid Mechanics 252, p. 209–238. Cited by: §2. [53] L. V. Rolandi, J. H. M. Ribeiro, C. Yeh, and K. Taira (2024) An invitation to resolvent analysis. Theoretical and Computational Fluid Dynamics 38 (5), p. 603–639. Cited by: §1. [54] O. Ronneberger, P. Fischer, and T. Brox (2015) U-net: convolutional networks for biomedical image segmentation. In International Conference on Medical image computing and computer-assisted intervention, p. 234–241. Cited by: §3. [55] C. W. Rowley and S. T. Dawson (2017) Model reduction for flow analysis and control. Annual Review of Fluid Mechanics 49 (1), p. 387–417. Cited by: §1. [56] C. W. Rowley, I. Mezić, S. Bagheri, P. Schlatter, and D. S. Henningson (2009) Spectral analysis of nonlinear flows. Journal of fluid mechanics 641, p. 115–127. Cited by: §1. [57] Y. Saad (1980) Variations on arnoldi’s method for computing eigenelements of large unsymmetric matrices. Linear algebra and its applications 34, p. 269–295. Cited by: §C.1. [58] P. J. Schmid and D. S. Henningson (2012) Stability and transition in shear flows. Vol. 142, Springer Science & Business Media. Cited by: §2, §3, §3. [59] P. J. Schmid (2007) Nonmodal stability theory. Annu. Rev. Fluid Mech. 39 (1), p. 129–162. Cited by: §1. [60] P. J. Schmid (2010) Dynamic mode decomposition of numerical and experimental data. Journal of fluid mechanics 656, p. 5–28. Cited by: §1, §1. [61] S. C. Shadden, F. Lekien, and J. E. Marsden (2005) Definition and properties of lagrangian coherent structures from finite-time lyapunov exponents in two-dimensional aperiodic flows. Physica D: Nonlinear Phenomena 212 (3-4), p. 271–304. Cited by: §3. [62] J. D. Simon and S. K. Mitter (1968) A theory of modal control. Information and control 13 (4), p. 316–353. Cited by: §1. [63] L. Sirovich (1987) Turbulence and the dynamics of coherent structures. i. coherent structures. Quarterly of applied mathematics 45 (3), p. 561–571. Cited by: §2. [64] S. Symon, K. Rosenberg, S. T. Dawson, and B. J. McKeon (2018) Non-normality and classification of amplification mechanisms in stability and resolvent analysis. Physical Review Fluids 3 (5), p. 053902. Cited by: §3. [65] K. Taira, S. L. Brunton, S. T. Dawson, C. W. Rowley, T. Colonius, B. J. McKeon, O. T. Schmidt, S. Gordeyev, V. Theofilis, and L. S. Ukeiley (2017) Modal analysis of fluid flows: an overview. AIAA journal 55 (12), p. 4013–4041. Cited by: §1. [66] V. Theofilis (2003) Advances in global linear instability analysis of nonparallel and three-dimensional flows. Progress in aerospace sciences 39 (4), p. 249–315. Cited by: §1. [67] V. Theofilis (2011) Global linear instability. Annual Review of Fluid Mechanics 43 (1), p. 319–352. Cited by: §1. [68] N. Thuerey, K. Weißenow, L. Prantl, and X. Hu (2020) Deep learning methods for reynolds-averaged navier–stokes simulations of airfoil flows. AIAA Journal 58 (1), p. 25–36. Cited by: §1. [69] A. Tran, A. Mathews, L. Xie, and C. S. Ong (2021) Factorized fourier neural operators. arXiv preprint arXiv:2111.13802. Cited by: §2. [70] L. N. Trefethen, A. E. Trefethen, S. C. Reddy, and T. A. Driscoll (1993) Hydrodynamic stability without eigenvalues. Science 261 (5121), p. 578–584. Cited by: §1. [71] K. Um, R. Brand, Y. R. Fei, P. Holl, and N. Thuerey (2020) Solver-in-the-loop: learning from differentiable physics to interact with iterative pde-solvers. Advances in neural information processing systems 33, p. 6111–6122. Cited by: §2. [72] E. Van der Avoird, H. Dijkstra, J. Nauw, and C. Schuurmans (2002) Nonlinearly induced low-frequency variability in a midlatitude coupled ocean–atmosphere model of intermediate complexity. Climate dynamics 19 (3), p. 303–320. Cited by: §4. [73] J. A. Weideman and S. C. Reddy (2000) A matlab differentiation matrix suite. ACM transactions on mathematical software (TOMS) 26 (4), p. 465–519. Cited by: §B.2. [74] R. Yu and R. Wang (2024) Learning dynamical systems from data: an introduction to physics-guided deep learning. Proceedings of the National Academy of Sciences 121 (27), p. e2311808121. Cited by: §1. [75] H. Yuan, J. Kou, C. Gao, and W. Zhang (2023) Resolvent and dynamic mode analysis of flow past a square cylinder at subcritical reynolds numbers. Physics of Fluids 35 (7). Cited by: §A.2. [76] Y. Zhao, D. B. Grayden, M. Boley, Y. Liu, P. J. Karoly, M. J. Cook, and L. Kuhlmann (2025) Cortical stability and chaos during focal seizures: insights from inference-based modeling. Journal of Neural Engineering 22 (3), p. 036021. Cited by: §4.