Conceptio › Archive › arXiv CS

arXiv CSopen access

Hypothesis-driven construction of mesoscopic dynamics

2026 · arxiv_cs

arXiv CS · Papers · License: Open Access · 2026

Open Source ↗Direct PDF ↓

neural-networks

machine learning, deep learning, neural networks

Hypothesis-driven construction of mesoscopic dynamics ∗ Zhuoyuan Li† §

Aiqing Zhu‡ §

Qianxiao Li† ‡ ¶

arXiv:2605.16211v1 [cs.LG] 15 May 2026

Abstract Traditional scientific modeling typically begins with fixed, instance-wise effective equations and then carries out equation-specific analysis and computation, a procedure that becomes exceptionally challenging in complex applications such as multiscale systems. We propose an alternative paradigm by learning mesoscopic dynamics within a mathematically constrained hypothesis class. Building upon a generalized Onsager principle, we introduce a unified framework encompassing both dissipative and conservative mesoscopic dynamics. We establish uniform and a priori theoretical guarantees, including global well-posedness, asymptotic stability, unique factorization identifiability, and discrete energy dissipation, applicable to all spatio-temporal evolution equations within this hypothesis class prior to all learning stages. Data from each problem instance is then used to guide the identification of members within our hypothesis class, giving rise to accurate, robust and interpretable dynamical models. We empirically validate this framework on both data from continuum PDE models as a check, and on data arising from microscopic chain models for which exact meso-scale models are unknown. The proposed approach not only acts as an effective dynamics learner, but also offers vital interpretable diagnostics of the underlying physics. Keywords: scientific computing, multiscale modeling, Onsager principle, interpretability

Introduction

Mathematical modeling of dynamical processes is traditionally organized around a well-established pipeline [16, 33, 36]: one starts from a physical law in the form of a mathematical equation with a finite set of parameters, uses experiments or observations to determine these parameters, thereby obtains a fixed governing equation; One then studies the well-posedness of the resulting model, designs stable and accurate discretizations, and analyzes the resulting simulations. This framework has been highly successful, and scientific modeling and computing have accordingly centered on two tasks: the identification of an effective governing equation [16, 24] and, once such an equation is available, its instance-wise mathematical and numerical analysis [5, 17]. In many problems of interest, the effective equation itself emerges from underlying microscopic descriptions through various limiting procedures, such as kinetic or hydrodynamic limits, homogenization or averaging arguments [16, 25]. However, it is often unclear how to derive the limiting equation in a systematic way, and even when such a derivation is available, favorable analytical properties of the approximating models need not pass automatically to the limit. This phenomenon is typical in multiscale dynamics, where the microscopic dynamics are often directly specified and analytically more accessible (yet expensive to simulate), whereas the corresponding mesoscopic spatio-temporal evolution equations are much more subtle to derive and analyze. Representative examples include ∗ This project is supported by the National Research Foundation, Singapore, under its AI Singapore Programme

(AISG Award No.: AISG3-RP-2022-028). Z.L. is supported by the Ministry of Education, Singapore, under its Research Centre of Excellence award to the Institute for Functional Intelligent Materials. (Project No. EDUNC-3318-279-V12). † Institute for Functional Intelligent Materials, National University of Singapore ([email protected], [email protected]) ‡ Department of Mathematics, National University of Singapore ([email protected]). § Zhuoyuan Li and Aiqing Zhu contributed equally to this work. ¶ Corresponding author.

Classical Modeling Framework model

Hypothesis-Driven Framework for Mesoscopic Dynamics Analysis

Hypothesis Space

Computation

model

Analysis

model model

model

Computation model

model model

Computation Interpretation

Analysis Computation

Analysis

Figure 1: Comparison of our framework and classical modeling framework. Classical approaches treat modeling and analysis case by case. Our approach instead first builds a hypothesis space with a unified analytical theory and then uses data to identify an admissible model from this space. the passage from particle chains to KdV-type equations [1, 41], from molecular dynamics to kinetic equations [16, 19], and from the asymmetric exclusion process to the KPZ equation [3, 14]. These considerations motivate us to leverage structured approaches and deep learning to circumvent these difficulties. Under suitable hypotheses, one can construct a class whose elements are admissible equations involving operators, functionals, or functions as unknown components. The class is chosen to balance two requirements: it must be sufficiently expressive for some elements of it to represent the underlying system with adequate fidelity, and sufficiently structured for well-posedness, stability, and related analytical properties to be studied uniformly across the class. Training is then used to identify a specific representative within such a class. In this way, one works within a subspace of possible equations that is both close enough to the target system for modeling purposes and sufficiently structured to support a family-level analysis. This change of framework is summarized in Figure 1. In this paper, we demonstrate this program for the specific class of systems admitting a possibly infinite-dimensional Onsager-type dynamics, formulated as ∂t u = − [M(u) + W(u)]

δV , δu

(1)

which will be explained later. In addition, we choose some structural hypotheses under which the operator constraints admit a pointwise spectral characterization in terms of Fourier multipliers. This defines a hypothesis class broad enough to include both dissipative and conservative thermodynamic systems, yet structured enough to admit a unified analysis. Within this setting, we prove global well-posedness and asymptotic stability. By securing these properties, we establish an a priori analytical foundation for the entire hypothesis class, ensuring that our models are endowed with physical stability independent of, and prior to, any data-driven learning process. We then establish identifiability of the Onsager factorization eq. (1) up to a one-dimensional rescaling, and derive explicit conditions under which the forward Euler scheme preserves the energy-dissipation inequality. For numerical validation, we introduce Spectral OnsagerNet, a neural architecture that enforces the structural constraints by construction. Numerical experiments on data generated from the KdV and Allen–Cahn equations confirm accuracy, stability, and recovery of the expected physical structure. Experiments on microscopic particle chains, including an FPUT chain [18] with a known KdV-type continuum limit [41] and an FENE chain [38] without a known closed-form continuum limit, further show that the method can recover interpretable mesoscopic dynamics directly from microscopic trajectory data. The remainder of this paper is organized as follows. After discussing our approach in the context of the existing literature in section 2, section 3 formally introduces the generalized Onsager principle on a Gelfand triple and establishes its spectral representation. Section 4 then provides the theoretical foundation for our framework, including rigorous analytical guarantees for global well-posedness, 2

asymptotic stability, factorization identifiability, and discrete energy dissipation. Transitioning from theory to computational methodology, section 5 details the neural architecture of Spectral OnsagerNet. Section 6 validates the proposed method through numerical experiments on both continuum PDE models and microscopic particle chains, and section 7 concludes the paper. For clarity of exposition, we only provide proof sketches for some of the propositions, and detailed versions are deferred to section C.

Related works

The identification of governing equations from observational data has been studied extensively [7, 12, 10, 6, 15]. Most existing works are formulated for finite-dimensional systems or for finite-dimensional state representations inferred from data. However, many problems of interest in scientific computing, including fluid flows, phase transitions, and pattern formation, are governed by evolution equations on infinite-dimensional state spaces. Accordingly, the object to be learned is not merely a vector field on Rd , but an evolution system posed on a function space. Recent data-driven approaches address such problems in two main ways. One line of work models the dynamics through finite-dimensional state-space or latent representations, as in ODE-Net [12, 11]. Another line learns the solution map or evolution operator more directly, as in neural operator methods such as the Fourier Neural Operator [28] and DeepONet [30]. While these approaches provide expressive approximations of nonlinear dynamics, they typically treat the learned vector field or operator as an unconstrained map. As a result, they generally do not identify structural features of the dynamics, such as an underlying energy, a dissipative–conservative splitting, or a built-in stability mechanism, which are often central to the analysis of infinite-dimensional evolution systems. A natural strategy to impose such structure is through physics-informed model classes rooted in variational or thermodynamic principles. Representative examples include Hamiltonian models for conservative systems [2, 9, 20, 21, 39] and Onsager- or GENERIC-type architectures for dissipative dynamics [34, 35, 31, 32, 40, 27, 42]. While related concepts have been successfully adapted for stochastic systems and model reduction [13, 43], these frameworks remain predominantly confined to finite-dimensional settings. A notable recent exception is Stat-PINNs [23], which incorporates thermodynamic constraints into infinite-dimensional evolution systems; however, it is strictly limited to equilibrium systems. Consequently, there remains a critical need for structure-preserving frameworks capable of learning non-equilibrium dynamics in infinite-dimensional spaces.

Generalized Onsager principle

The objective of this paper is to model the evolution of physical fields of our interest over a spatial domain, specifically targeting mesoscopic dynamics. We characterize this mesoscopic regime as the continuous spatio-temporal evolution that bridges the critical gap between discrete, microscopic particle interactions and bulk, macroscopic thermodynamics. By capturing these intricate spatial and temporal dependencies, we seek to develop a mesoscopic model that is globally wellposed, thermodynamically consistent, and stable for long-time evolution. This section establishes the continuous mathematical framework necessary to achieve this goal. We first extend the classical finite-dimensional Onsager principle to infinite-dimensional spaces. After that, we introduce the core structural hypotheses, including the rigorous formulation of the underlying function spaces, that render the hypothesis class analytically and computationally tractable. Finally, we derive a spectral representation of these dynamics, bridging the continuous thermodynamic theory with the practical neural network implementation.

3.1

From finite to infinite dimensions

In finite dimensions, learning from trajectories is often formulated as learning an ODE vector field żt = f (zt ). The Neural ODE approach parameterizes f as an unconstrained neural network [12], 3

which offers expressive power but does not generally preserve the structure underlying physical systems. OnsagerNet provides a structured alternative for finite-dimensional dissipative dynamics [40, 13]. Rather than treating the vector field as an unconstrained black box, it decomposes the dynamics into a dissipative part and a conservative part, thereby extending the classical Onsager principle [34] to strongly non-equilibrium regimes. Specifically, the evolution for the time-dependent state variable zt is modeled as żt = −[M (zt ) + W (zt )]∇V (zt ),

(2)

where M (zt ) is a symmetric and positive semi-definite matrix, W (zt ) is an anti-symmetric matrix, and V (zt ) is a scalar function modeling a potential energy. Each component is parameterized by a neural network and inferred from data. By embedding the dissipative–conservative decomposition into the parameterization, OnsagerNet enforces thermodynamic consistency and improves long-term stability, which is generally not guaranteed by unconstrained Neural ODE approaches. These features make it suitable for learning finite-dimensional dissipative dynamics, with applications including Rayleigh–Bénard convection and polymer-chain stretching dynamics [12, 40]. We next extend this principle to mesoscopic dynamics, where the evolving state is a spatial field rather than a finite-dimensional vector. Typical examples include strain fields in nonlinear particle chains, density fields in interacting particle systems and phase fields in phase-separation models. Such field-level descriptions provide spatially resolved effective closures of microscopic dynamics, thereby revealing how local microscopic laws generate wave propagation, energy transfer, pattern formation, phase separation, and relaxation at larger scales. Consider a free-energy functional V : X → R and its Fréchet derivative µ(u) = δV /δu ∈ X ∗ , which plays the role of the generalized thermodynamic driving force [37]. We define the linear operators M(u) and W(u) mapping from the dual space (forces) to the primal space (fluxes/rates): M(u), W(u) : X ∗ → X .

(3)

The generalized Onsager principle in infinite dimensions then takes the form ∂t u = − [M(u) + W(u)] µ(u),

µ(u) =

δV ∈ X ∗. δu

(4)

However, unlike the finite-dimensional setting, the infinite-dimensional formulation involving statedependent functionals and operators poses significant computational and theoretical challenges. A direct implementation via spatial discretization leads to prohibitively high-dimensional models, making naive parameterization and training costly. Enforcing symmetry or anti-symmetry and positive semi-definiteness adds additional constraints that further increase the computational burden. From a theoretical perspective, it is also nontrivial to ensure that the learned model is globally well-posed and remains stable over long-time evolution. These challenges motivate the structural hypotheses we introduce next, which preserve the Onsager structure at the operator level, admit efficient spectral representations, and provide principled guarantees of energy dissipation and long-time stability.

3.2

Structural hypotheses

We aim to study a broad dynamic class of the following form ∂t u(t) = G(u(t)),

u(0) = uinit ∈ X

(5)

with the solution space W (0, T ; X ) := u ∈ L2 (0, T ; X ) ∂t u ∈ L2 (0, T ; X ) .

(6)

In this subsection, we formulate the structural hypotheses used in the subsequent analysis, which keep the hypothesis class analytically tractable while preserving sufficient approximation flexibility. We start by stating several mild hypotheses on the function space X , where the field of the physical mesoscopic state of our interest stays. 4

Hypothesis 3.1. X is a uniformly smooth and uniformly convex Banach space compactly embedded within H = L2 (Ω), and the solution domain Ω is bounded in Rd . Throughout this work, we use (·, ·) to denote the inner product on the Hilbert space H = L2 (Ω) and ⟨·, ·⟩X ∗ ,X to denote the duality pairing on X . For a Banach space X compactly embedded within H = L2 (Ω), it is clear that we have the Gelfand triple X ,→ H ∼ = H∗ ,→ X ∗ .

(7)

Moreover, uniformly convexity implies the reflexivity of X by Milman–Pettis theorem. Sobolev spaces offer a concrete class of Banach spaces satisfying Hypothesis 3.1. A Sobolev space W s,p (Ω) is uniformly smooth and uniformly convex for 1 < p < +∞, so it suffices to ensure that W s,p (Ω) is compactly embedded within L2 (Ω), which is a direct result of the Sobolev embedding theorem. In particular, we need s ≥ d(1/p − 1/2)+ for the Sobolev space W s,p (Ω) to satisfy our hypothesis. Having established the functional spaces, we now introduce the fundamental structural hypotheses governing the dynamics. These conditions are specifically designed to enforce thermodynamic consistency directly at the operator level. Hypothesis 3.2 (Dissipation structure). The nonlinear operator G : X → X admits a decomposition G(u) = −L(u)µ(u)

and

L(u) = M(u) + W(u).

(8)

Here, for each u ∈ X , µ(u) ∈ X ∗ , and L(u) is a linear operator from X ∗ to X . The operators M, W and µ satisfy 1. M(u) : X ∗ → X is symmetric and positive-semidefinite: ⟨ξ, M(u)η⟩X ∗ ,X = ⟨η, M(u)ξ⟩X ∗ ,X

∀ξ, η ∈ X ∗

(9)

and ⟨ξ, M(u)ξ⟩X ∗ ,X ≥ 0

∀ξ ∈ X ∗ .

(10)

2. W(u) : X ∗ → X is skew-symmetric: ⟨ξ, W(u)η⟩X ∗ ,X = − ⟨η, W(u)ξ⟩X ∗ ,X .

(11)

3. µ(u) ∈ X ∗ is the Fréchet derivative of a coercive functional: there exists a Fréchet differentiable functional V : X → R such that • V (un ) < +∞ implies ∥un ∥X < +∞ for any sequence {un }n ⊆ X ; • µ(u) = δV /δu as the Fréchet derivative of V . These structural hypotheses are a direct infinite-dimensional generalization of the finite-dimensional Onsager principle eq. (2). Physically, the functional V acts as a generalized free energy, while the operators M and W capture the dissipative and conservative transport mechanisms, respectively. A vital consequence of this imposed structure is the continuous dissipation of free energy, which guarantees thermodynamically consistent macroscopic behavior. Proposition 3.3 (Energy dissipation). With Hypothesis 3.2, the free energy does not increase along any solution trajectory as d V (u) = ⟨µ(u), ∂t u⟩X ∗ ,X = − ⟨µ(u), M(u)µ(u)⟩X ∗ ,X ≤ 0. dt

(12)

The skew-symmetric part W does not contribute to the energy evolution. Hence, any solution trajectory, if it exists, is uniformly bounded in X due to the coercivity of V .

While Proposition 3.3 guarantees that trajectories remain uniformly bounded due to the coercivity of V , establishing the global well-posedness of the dynamics class eq. (5) requires two additional hypotheses. The first assumption is physically motivated by translational invariance on periodic domains, though formulated here as a general property of linear operators. Crucially, this convolution structure not only makes the continuous analysis tractable but also admits a diagonalized representation in the frequency domain. This explicitly bridges the continuous theory with our computational method, directly inspiring the spectral training and evolution of Spectral OnsagerNet. Meanwhile, following standard well-posedness arguments, we impose local Lipschitz continuity to ensure the existence and uniqueness of the solution flow. Hypothesis 3.4 (Convolution structure). The operators M and W act as convolutions in the physical domain Ω = Td . Namely, there exist KM (·; u), KW (·; u) ∈ L2 (Ω) such that for any η ∈ X ∗ , M(u)(η)(x) = ⟨η, KM (x − ·; u)⟩X ∗ ,X , W(u)(η)(x) = ⟨η, KW (x − ·; u)⟩X ∗ ,X .

(13)

Hypothesis 3.5. The operators L and µ are X -Lipschitz continuous and the operator G is L2 Lipschitz continuous on bounded sets. Formally, there exist L1 (R), L2 (R), and L3 (R) such that ∥µ(v1 ) − µ(v2 )∥X ∗ ≤ L1 (R)∥v1 − v2 ∥X , ∥L(v1 ) − L(v2 )∥X ∗ →X ≤ L2 (R)∥v1 − v2 ∥X ,

∀ v1 , v2 ∈ BX (R),

(14)

∥G(v1 ) − G(v2 )∥2 ≤ L3 (R)∥v1 − v2 ∥2 , where BX (R) := {v ∈ X | ∥v∥X ≤ R}. The generalized Onsager principle eq. (4) provides a unifying framework for a diverse array of infinite-dimensional dynamics. Its scope encompasses many classical PDE dynamics, including purely dissipative systems (e.g., the Allen–Cahn and Cahn–Hilliard equations), purely conservative systems (e.g., the KdV and nonlinear Schrödinger equations), and mixed dynamics typified by the incompressible Navier–Stokes equations. Since the generalized Onsager principle is formulated at the abstract operator level, it extends naturally to non-local integro-differential operators and the effective mesoscopic limits of discrete microscopic systems, such as the nonlinear particle chains explored later in our numerical experiments. More details of these examples are exhibited in section B.

3.3

Spectral representation

A key consequence of Hypothesis 3.4 is that the operators M(u) and W(u) diagonalize in Fourier space. Theorem 3.6. Under Hypothesis 3.4, we define the Fourier series c[k] (u) = (KM (·; u), ek ) M

c[k] (u) = (KW (·; u), ek ) , W

and

∀k ∈ Z

(15)

for the kernels KM (·; u) and KW (·; u), respectively, where {ek }k∈Z denotes the Fourier basis of L2 (Ω). Then it holds for M (and W, respectively) that X c[k] (u)η̂ [k] ek → M(u) M (16) |k|≤N

in X as N → +∞ and ⟨ξ, M(u)η⟩X ∗ ,X =

c[k] (u)η̂ [k] ξˆ[−k] M

(17)

k∈Z

for any ξ and η ∈ X , where we set ξˆ[k] = ⟨ξ, ek ⟩X ∗ ,X

η̂ [k] = ⟨η, ek ⟩X ∗ ,X ,

and

∀k ∈ Z.

(18)

Furthermore, Hypothesis 3.2 holds if and only if c[k] (u) = M c[−k] (u) ≥ 0 M

and

for all k ∈ Z. 6

c[k] (u) = W c[−k] (u) ∈ iR W

(19)

Proof. By definition, Z ⟨η, KM (x − ·; u)⟩X ∗ ,X ek (x) dx Ω Z = η, KM (x − ·; u)ek (x) dx

(M(u)(η), ek ) =

X ∗ ,X

(20)

Z η, ek (·)

KM (x − ·; u)ek (x − ·) dx X ∗ ,X

c[k]

[k]

= ⟨η, (KM (u), ek ) ek ⟩X ∗ ,X = M (u)η̂ . The second equality holds because the operator µ is continuously linear and the mapping x 7→ ek (x)KM (x − ·; u) is Bochner integrable due to the boundness of Ω = Td and the L2 -integrability of KM . Similar deduction can be applied to W as well, and all the conclusions follow immediately by direct evaluation. This theorem lies at the heart of our proposed framework. By establishing that the stringent infinite-dimensional constraints of symmetry, positive semi-definiteness, and anti-symmetry (Hypothesis 3.2) are strictly equivalent to pointwise conditions on Fourier multipliers, the problem is rendered both analytically and computationally tractable. The spectral characterization also serves as the foundational bridge connecting the theoretical guarantees derived in section 4 with the practical network parameterization detailed in section 5.

Theoretical analysis

In this section, we establish a rigorous theoretical foundation for the generalized Onsager dynamics (4). Throughout our analysis, all mathematical proofs inherently rely on the structural hypotheses formalized in Hypotheses 3.2, 3.4, and 3.5. We begin by proving the global well-posedness of the system and characterizing the long-time asymptotic behavior of its solution trajectories. To guarantee that the learned operators carry genuine physical meaning rather than acting as unconstrained fitting devices, we then establish the identifiability of the Onsager factorization when M and W are independent of the state u. Finally, bridging the continuous formulation with its practical implementation, we derive explicit stability conditions that ensure the time-discrete approximation strictly preserves the energy dissipation property.

4.1

Global well-posedness

We begin with the most basic analytical requirement for the proposed hypothesis class by establishing global well-posedness under the structural hypotheses. Theorem 4.1 (Global existence). For any initial condition u(0) = uinit ∈ X , the Cauchy problem eq. (5) admits a global solution u defined on [0, +∞). Proof sketch. To construct a global solution, we first define a sequence of finite-dimensional approximate solutions {un }n , establishing their local existence via Peano’s theorem. We then leverage the energy dissipation property of the operator L according to Hypotheses 3.2 and 3.4, combined with the coercivity of the functional V , to derive a uniform bound for this sequence. Applying the Aubin– Lions lemma alongside a diagonal extraction process yields a weakly convergent subsequence. By invoking the Lipschitz continuity stipulated in Hypothesis 3.5, we demonstrate that this weak limit is indeed a global solution to the Cauchy problem. A detailed proof can be found in section C.1. Let J : X → X ∗ denote the duality mapping, then the uniform smoothness of X guarantees that J is single-valued and uniformly continuous on bounded sets. The mapping J can also be treated as the Fréchet derivative of 21 ∥ · ∥2X .

Lemma 4.2. For any u ∈ X , the functional J(u) ∈ X ∗ satisfies ⟨J(u), u⟩X ∗ ,X = ∥u∥2X ,

∥J(u)∥X ∗ = ∥u∥X .

(21)

The proof is deferred to section C.2, leading to the following uniqueness result. Theorem 4.3 (Global uniqueness and continuous dependence). The global solution described in Theorem 4.1 is unique. Moreover, if ũ(t) and ū(t) are two solutions with initial conditions ũinit and ūinit respectively, then 2

∥ũ(t) − ū(t)∥X ≤ e2L4 (R)t ∥ũinit − ūinit ∥X ,

∀t > 0,

(22)

where R is a uniform bound on ∥ũ(t)∥X and ∥ū(t)∥X , and L4 (R) is the X -Lipschitz constant for G. Proof. By Proposition 3.3 and the coercivity of V , both trajectories are uniformly bounded in X . Let R denote this common bound. Then d ∥ũ(t) − ū(t)∥2X = 2 ⟨J(ũ(t) − ū(t)), ∂t ũ(t) − ∂t ū(t)⟩X ∗ ,X dt = 2 ⟨J(ũ(t) − ū(t)), G(ũ(t)) − G(ū(t))⟩X ∗ ,X (23) ≤ 2 ∥J(ũ(t) − ū(t))∥X ∗ ∥G(ũ(t)) − G(ū(t))∥X ≤ L4 (R)∥ũ(t) − ū(t)∥2X , where the last inequality follows from Lemma 4.2 and the X -Lipschitz condition of G. The conclusion follows by Grönwall’s inequality. Uniqueness is the special case ũinit = ūinit .

4.2

Asymptotic behavior

After establishing global well-posedness, we turn our attention to the long-time asymptotic behavior of the system. According to Proposition 3.3, the free-energy functional V acts strictly as a Lyapunov function for the dynamics. Consequently, the asymptotic limits of the trajectories are linked to the kernel structure of the dissipation operator M. To formalize this relationship, we distinguish between two critical invariant sets. Definition 4.4 (Equilibrium and zero-dissipation sets). We define the following subsets of X : 1. E as the set of equilibria, consisting of states where the system is stationary. E = {φ ∈ X | G(φ) = −[M(φ) + W(φ)]µ(φ) = 0}.

(24)

2. Z as the zero-dissipation set, consisting of states where the free energy dissipation vanishes. n o Z = φ ∈ X ⟨µ(φ), M(φ)µ(φ)⟩X ∗ ,X = 0 . (25) Clearly, E ⊆ Z. By virtue of the global well-posedness, the generalized Onsager dynamics induces a well-defined continuous flow S(t) : uinit 7→ u(t). The following theorem demonstrates that every trajectory asymptotically approaches the zero-dissipation set, which, under specific strict positivity conditions, collapses exactly to the set of stationary equilibria. Theorem 4.5. By global well-posedness in section 4.1, the solution induces a continuous flow S(t) : u(0) 7→ u(t). For any trajectory u(t), the ω-limit set \ ω(u(0)) := {u(t) : t ≥ T } (26) T >0

is a non-empty and compact subset of Z, which is positively invariant under S(t). Moreover, if the dissipation operator M(u) satisfies the elliptic condition, which means there exists c > 0 such that 2

⟨ξ, M(u)ξ⟩X ∗ ,X ≥ c ∥ξ∥X ∗ ,

∀ξ ∈ X ∗ ,

then Z = E, and any trajectory converges asymptotically to the set of static equilibria. 8

(27)

Proof. A detailed proof can be found in section C.3. According to Proposition 3.3, any trajectory {u(t)}t∈[0,+∞) is uniformly bounded in X . By the Rellich–Kondrachov Theorem, the embedding X = H 1 (Ω) ,→ L2 (Ω) is compact. Therefore, we have a non-empty ω-limit set ω(u(0)). We claim that the energy V is a constant on the ω-limit set. Otherwise there exist two sequence of snapshots with energy converging to distinct values, which contradict the monotonicity of V along the trajectory. Consider a trajectory lying on the ω-limit set. By taking the time derivative of V , we can easily verify that the free energy dissipation vanishes, and consequently we have ω(u(0)) ⊆ Z. Furthermore, if the elliptic condition eq. (27) holds, then for any φ ∈ Z, we have that ⟨µ(φ), M(φ)µ(φ)⟩X ∗ ,X = 0 implies µ(φ) = 0, so φ ∈ E. In the general case where M(u) has a nontrivial kernel, the zero-dissipation set Z may be strictly larger than the equilibrium set E. Physically, this means that rather than settling into a static equilibrium, the system may converge to a persistent dynamic regime such as motion along a Hamiltonian orbit on a constant-energy surface driven entirely by the conservative operator W. The significance of the convergence result lies in its strictly a priori nature. We do not rely on the data-driven optimization process to discover, approximate, or enforce long-term stability; instead, up to fitting and discretization errors, any learned model is mathematically ensured to converge to a physically meaningful orbit. Consequently, the discovered mesoscopic dynamics remain physically interpretable and structurally sound over infinite time horizons. Definition 4.6. In fact, we have two physically distinct regimes of the generalized Onsager dynamics eq. (5): 1. Elliptic dissipative systems (e.g., Allen–Cahn): M satisfies the elliptic condition as in eq. (27), the free energy decays strictly, and trajectories converge to static equilibria. 2. Conservative systems (e.g., KdV and other Hamiltonian PDEs): M ≡ 0, the free energy is exactly conserved, and trajectories may exhibit persistent dynamics on constant-energy surfaces. While intermediate regimes naturally exist where M is non-zero but possesses a non-trivial kernel, violating the elliptic condition, a detailed theoretical treatment of such mixed dynamics falls outside the scope of the present study. Instead, we focus on the two principal extremes defined above. Crucially, as we demonstrate in the subsequent section, these regimes dictate entirely different identifiability properties for the learned Onsager factorization.

4.3

Uniqueness of the Onsager factorization

Recall that a central objective of our framework is to recover interpretable mesoscopic mechanisms. Hence, it would be desirable if the operators (M and W) and free energy functionals V could be uniquely identified. While in general there may be a gauge dependence as shown in [43], we demonstrate that when the operators M and W are state-independent, the generalized Onsager factorization eq. (4) is uniquely determined under mild spectral irreducibility conditions. Definition 4.7 (spectral irreducibility). Consider a functional V : X → R. We say that the functional V is irreducible on the spectral space if it is second-order Fréchet differentiable, and there exists no partition Z = S1 ⊔ S2 such that δ2 V d2 (e , e ) = V (u + τ1 ep + τ2 eq ) ≡ 0, p q δu2 dτ1 dτ2 τ1 =τ2 =0

∀p ∈ S1 , q ∈ S2 ,

(28)

where {ek }k∈Z is the spectral basis of H = L2 (Ω). To build intuition for the spectral irreducibility condition, it is sufficient to note that nonlinear interaction within the integrand of the free energy potential V will probably mix the spectral modes. For instance, the potential V for the KdV dynamics has a cubic component u3 within the integrand, whose computation mixes the spectral modes and thus the functional naturally satisfies the irreducibility requirement. 9

Theorem 4.8. Let (M1 , W1 ) and (M2 , W2 ) be independent of u and satisfy (M1 + W1 )µ1 (u) = (M2 + W2 )µ2 (u),

µ1 =

δV1 , δu

µ2 =

δV2 δu

(29)

with our previous hypotheses. Supposing that V2 is irreducible on the spectral space, we have • if M1 is positive definite and V2 is not a constant functional, then there exists a constant λ > 0 such that M2 = λM1 , W2 = λW1 , µ1 = λµ2 ; (30) [k]

• if M2 = 0 and µ̂2 ̸≡ 0 for all k ∈ Z, then M1 µ1 (u) ≡ 0. Furthermore, if ker W1 = 0, then µ1 = λ′ µ2 for a constant λ′ ∈ R. The detailed proof can be found in section C.5. Remark 4.9. Crucially, the irreducibility requirement is a natural and mild assumption rather than a restrictive limitation. If a functional V is reducible, it physically implies that the state variable can be partitioned into two or more mutually non-interacting subsystems. In such cases, the dynamics can simply be decomposed into their respective irreducible components. The uniqueness result of Theorem 4.8 can then be applied independently to each block. Consequently, the generalized Onsager factorization remains unique up to a set of independent, one-dimensional scaling factors, with exactly one scaling constant corresponding to each decoupled physical subsystem. Theorem 4.8 offers identifiability guarantees for the generalized Onsager principle when M and W are state-independent. Interpreting the left-hand side of eq. (29) as the learned dynamics and the right-hand side as the reference, two physically distinct regimes emerge. 1. Elliptic dissipative systems: Every learned factorization reproducing the reference dynamics recovers M, W, and µ simultaneously, each up to a common positive scalar. The dissipative and conservative operators are therefore individually identifiable, and the learned energy is an affine transformation of the true one. 2. Conservative systems: Theorem 4.8 forces the learned dissipation contribution M1 µ1 to vanish identically, regardless of M1 itself. If additionally ker W1 = 0, the learned functional derivative coincides with the true one up to a real scalar, so the learned energy is conserved along every trajectory. In particular, the framework cannot spuriously inject dissipation into a genuinely conservative system. In both regimes, the Onsager factorization is unique modulo a one-dimensional rescaling. For conservative systems an additional freedom remains in M1 , but it is physically inert since its contribution to the dynamics vanishes. These guarantees ensure that the learned energy functional and transport operators carry genuine physical meaning rather than serving as unconstrained fitting devices, thereby providing a rigorous foundation for interpreting the outputs of Spectral OnsagerNet as discovered governing laws. Some experimental evidence for the identifiability can be found in section 6.2.2.

4.4

Discrete energy dissipation

While the continuous formulation guarantees thermodynamic consistency, practical implementation requires evaluating the system on a computational grid and advancing it via a numerical integrator. In our framework, we employ the explicit forward Euler method, where continuous stability does not automatically ensure numerical stability due to the discrepancies introduced by spatial and temporal discretization. Let {usp }p∈G denote the discrete state at time level s, where usp = up (s∆t) and p indexes the grid points in a fixed grid G. We consider the time-discrete scheme us+1 = us − ∆t[M(us ) + W(us )]µ(us ).

(31)

Our goal is to identify conditions on the time step ∆t under which the discrete evolution inherits an energy-dissipation property. The estimate is obtained by applying the descent lemma to the discrete increment from us to us+1 . Theorem 4.10 (Discrete energy dissipation). The discrete evolution satisfies the energy dissipation property V (us+1 ) ≤ V (us ) if the time step ∆t ≤ 2κL1 (R0 )−1 , where κ=

2 ⟨µs , Ms µs ⟩X ∗ ,X

(32)

∥(Ms + W s )µs ∥X

and R0 depends only on V . Here, we use the abbreviations Ms ≡ M(us ), W s ≡ W(us ), and µs ≡ µ(us ). Moreover, • for elliptic dissipative dynamics where eq. (27) holds, the parameter κ ≥ κ0 > 0, and the constant κ0 depends only on M, W, and V . • for conservative dynamics where M = 0, the one-step variation of the potential V (us+1 ) − V (us ) ≤ C(R̃)(∆t)2 = O(∆t2 ),

(33)

where the parameter R̃ depends on W, V , and the maximal X -norm across the discrete trajectory. It follows that the potential changes V (us ) − V (u0 ) on the whole trajectory is O(∆t) for any s∆t ≤ T if the trajectory does not blow up within finite time T . Proof sketch. By the Lipschitz continuity of µ1 , applying the descent lemma we have V (us+1 ) − V (us ) ≤ ⟨µs , Ms µs ⟩X ∗ ,X ∆t + O(∆t2 )

(34)

max ∥us ∥X , us+1 X ≤ R.

(35)

as long as Then, we can prove by induction that either µs = 0, or ∆t ≤

2 ⟨µs , Ms µs ⟩X ∗ ,X 2

∥(Ms + W s )µs ∥X L1 (Rs )

and

Rs = D +

2 ⟨µs , Ms µs ⟩X ∗ ,X ∥(Ms + W s )µs ∥X L1 (D)

(36)

for the sth step will ensure that V (us+1 ) ≤ V (us ) and ∥us ∥X ≤ D for all s. The constant D is chosen such that V (w) ≤ V (u0 ) =⇒ ∥w∥X ≤ D (37) for any w ∈ X , and the existence is guaranteed by the coercivity of V . The conclusion for elliptic dissipative dynamics and conservative dynamics follows respectively by direct deduction. A detailed proof can be found in section C.6. Consequently, for elliptic dissipative systems, the explicit forward Euler scheme strictly preserves the monotonic decay of the free energy when the time step is sufficiently small. For conservative systems, while exact energy preservation is relaxed due to time discretization, the discrete energy fluctuates by at most an O(∆t2 ) for each update. This result reflects a fundamental philosophical departure from the traditional scientific computing and numerical analysis workflows. In classical settings, ensuring long-term stability and thermodynamic consistency typically requires the rigorous design of complex, equation-specific discretizations. In contrast, our framework reverses this by implicitly imposing these stability constraints by defining a rigorous continuous hypothesis class prior to any data-driven identification of equations. Consequently, even the simplest explicit integrator automatically inherits these robust physical guarantees. This decoupled approach not only provides critical flexibility for neural network model design, but also enables a unified, computationally efficient temporal discretization that applies universally across an entire class of unknown, data-driven mesoscopic dynamics without requiring instance-specific numerical engineering. 11

Spectral OnsagerNet

Starting from this section, we focus our discussion on a specific neural network implementation of the framework, where we set X = H 1 (Ω). We name the architecture Spectral OnsagerNet, abbreviated as SpecOnsNet. To construct a network for learning the generalized Onsager principle, we take the ansatz Z Z h i β α 2 (38) V (u) = |u| + F (u(x)) dx + |∇u|2 dx +V2 (u) Ω 2 Ω 2 {z } | {z } | V0 (u)

V1 (u)

for the free energy V with positive scalars α and β. Here, V0 and V1 account for pointwise potentials and local interactions, and V2 : X → R is a functional capturing the residuals of the integrands. The motivation for this explicit decomposition comes from physics. In fact, the ansatz mirrors the standard structure of thermodynamic free energies: the β-term penalizes steep spatial gradients to account for interfacial energy, capillarity, or dispersion when coupled with a conservative transport operator; α provides a global quadratic stabilizing potential that ensures the coercivity of the functional; the pointwise function F (u) models local energy contributions; and V2 (u) serves as a flexible residual functional designed to capture complex, non-local, or unknown mesoscopic interactions. We then obtain the functional derivative µ(u) =

δV2 δV = αu − β∆u + F ′ (u) + . δu δu

(39)

To ensure that the preceding theoretical analysis applies to our ansatz eq. (38), we must impose specific constraints on the functional V . More specifically, we let α, β > 0 be positive scalars and define Z Q(u) = F (u(x)) dx + V2 (u). (40) Ω

We further enforce that the functional Q : X → R is bounded from below and that its Fréchet derivative δQ/δu is X -Lipschitz continuous on bounded set. These conditions guarantee that both the coercivity requirement in Hypothesis 3.2 and the X -Lipschitz continuity in Hypothesis 3.5 are satisfied. For network implementation, our network is designed to follow the update eq. (31) in spectral space. More specifically, we consider the explicit Euler step h i c[k] (us ) + W c[k] (us ) µ̂s,[k] ûs+1,[k] = ûs,[k] − ∆t M (41) for each spectral mode, where we define µ̂s,[k] as the kth mode for µ(us ). Since δV2 µ̂[k] := ⟨µ(u), ek ⟩X ∗ ,X = ⟨αu − β∆u + F ′ (u), ek ⟩X ∗ ,X + , ek δu X ∗ ,X [k] [k] ∂v2 ′ (u) \ + = α + (2πk)2 β û[k] + F , ∂ û

(42)

where v2 is the spectral representation of V2 , we have to specify appropriate parameterizations for the terms M, W, F , and v2 . Convolutions M and W To enforce the constraints in Hypothesis 3.2, it suffices, by Theorem 3.6, cψ (us ) to be a non-negative real number and each entry of W cψ (us ) to be to require each entry of M c and W c as purely imaginary . Consequently, we parameterize M cψ (u) = ℜ (Gψ (û))2 M

cψ (u) = iℑ (Gψ (û)) , and W

where ℜ(·) and ℑ(·) denote the real and imaginary part of the network output Gψ (û).

(43)

Free energy functional V With the ansatz given by eq. (38), we directly parameterize the functional variation µ = δV /δu as δVϕ , δu then according to eq. (42), the spectral modes of µ are written as ′ \ µ̂ϕ (u) = α + (2πk)2 β ⊙ û + F ϕ (u) + ∇vϕ (û), µϕ (u) = αu − β∆u + Fϕ′ (u) +

(44)

(45)

where k is the Fourier mode-index vector, with the entry corresponding to mode k equal to k. As explained previously, to make our theoretical analysis applicable for our network, we only need to force α, β > 0 with a softplus layer and control the boundedness of the networks since the remaining continuity assumptions are automatically satisfied for common network architectures. In particular, we adopt a sinusoidal activation function for the final layers. The Fourier spectral modes are computed with the discrete Fourier transform, and the gradient is evaluated by autodifferentiation. In summary, we use the update h i 2 ûs+1 ≈ Tθ (ûs ) := ûs − ∆t ℜ (Gψ (û)) + iℑ (Gψ (û)) ⊙ µ̂ϕ (us ), (46) where ′ s s µ̂ϕ (us ) = α + (2πk)2 β ⊙ ûs + F\ ϕ (u ) + ∇vϕ (û ).

(47)

The parameter family θ contains all the trainable parameters ψ, ϕ, α, and β.

Numerical experiments

We empirically validate the proposed framework in two complementary settings: (i) learning from data generated from PDE models, where the underlying dynamics are known and serve as ground truth; and (ii) learning effective mesoscopic dynamics directly from microscopic simulation data, where no closed-form macroscopic equation is available a priori. The first setting allows us to rigorously assess predictive accuracy, long-time stability, and the identifiability properties established in section 4, while the second setting demonstrates that our framework can serve as a practical tool for discovering interpretable mesoscopic models from microscopic data.

6.1

Network training

The evolution of the system is solved in the spectral domain using a simple explicit forward Euler scheme. To optimize the network parameters θ, we compute the reconstruction error for each reference trajectory {us }Ss=1 using a K-step accumulated loss function, evaluated as J (θ) =

S−K K XX

Tθk (ûs )∨ − us+k 2 ,

(48)

s=1 k=1

where Tθ represents the Euler update specified in eq. (46). The symbol “∨” denotes the inverse discrete Fourier transform, so that the norm is computed in physical space. We adopt the multistep loss specifically to penalize accumulating integration errors and enhance long-term predictive stability, but importantly, our framework is not limited to this exact learning approach. Any other objective capable of training the dynamics, such as the detaching trick used in [4, 29], can be applied without compromising the rigorous thermodynamic structure.

6.2

Learning from known PDE dynamics

We first evaluate the predictive capabilities and numerical stability of the proposed framework on two known PDE dynamics: the Allen–Cahn model (both 1D and 2D cases), representing purely dissipative phase-field dynamics, and the KdV model, representing conservative dispersive wave propagation. Readers may refer to section D for more details about these PDE models. 13

Figure 2: The general workflow for our Spectral OnsagerNet. 100

Relative RMSE

10 1

10 2

10 3

FNO Allen-Cahn (classical solver) OnsagerNet SpecOnsNet (ours)

Prediction step

10 1

10 2

10 3

100

FNO KdV (classical solver) OnsagerNet SpecOnsNet (ours)

Prediction step

100

(b) KdV dataset

(a) Allen–Cahn (1D) dataset

Figure 3: Comparison of long-time prediction performance 6.2.1

Predictive accuracy and long-term stability

To demonstrate the capability of SpecOnsNet, we benchmark our model with 3 different baselines: (i) A classical one-step numerical scheme analogous to that for generating the reference trajectory, but with a larger (×25) step size for a fair comparison; (ii) A data-driven approach based on the Fourier Neural Operator (FNO) [28] without any structural assumptions; (iii) The original OnsagerNet [40] designed for finite-dimensional cases, where we fix the spatial discretization and restate the learning task in the discrete space. The general workflow of our model is displayed in fig. 2. Detailed description about the baselines as well as the evaluation metric can be found in sections A.1 and A.2. As illustrated in Figure 3, the rigorous enforcement of the generalized Onsager structure of our model serves as a highly effective inductive bias. Compared to the FNO model as a purely black-box learning approach and the OnsagerNet model as a naive generalization to PDE cases, SpecOnsNet demonstrates superior predictive accuracy and robust stability over extended integration horizons. Note that we do not apply the classical OnsagerNet to the 2D Allen–Cahn model due to its impractical memory cost since we would need around 100 times more budgets than those for the 1D case for a 128 × 128 grid. The advantage of our SpecOnsNet is also quantitatively confirmed in table 1, where SpecOnsNet consistently achieves the lowest 5-step relative prediction error with the minimal number of parameters. 6.2.2

Identifiability of the dynamics

Beyond mere trajectory matching, a core objective of our framework is the faithful extraction of the underlying thermodynamic mechanisms. To this end, we investigate the temporal behavior of the learned potentials evaluated on the reference trajectories to validate the identifiability guarantees established in section 4.3. As depicted in fig. 4, for the dissipative Allen–Cahn system, the learned free-energy functional Vθ exhibits a strictly affine relationship (highlighted by the dashed linear fit) with the true physical potential V . This empirical result corroborates the theoretical prediction that,

Task

Classical solver

FNO[28]

OnsagerNet[40]

SpecOnsNet

#params (1D; 2D)

—

(66.3 K; 1.3 M)

(1.3 M; N/A)

(18.2 K; 50.0 K)

Allen–Cahn (1D) KdV Allen–Cahn (2D)

±0.0029

0.0085 0.0541±0.0313 0.0074±0.0008 — —

FPUT chain FENE chain

0.0204 0.1747±0.0820 0.0311±0.0047

±0.0037

0.0260 0.0671±0.0241 Out of Memory

0.0055±0.0022 0.0100±0.0037 0.0037±0.0011

0.0360±0.0096 0.0186±0.0066

0.0242±0.0075 0.0139±0.0076

0.0082±0.0026 0.0012±0.0004

±0.0081

Table 1: 5-step relative prediction error (↓) for various baselines.

Learned potential V (u(t))

Learned potential V

1250

fit: V = 1.85V + 137.97

200 150 100 50 0 50 100 150

1000 750 500 250 0 250 500

150

125

100

Real potential V

(a) Allen–Cahn (1D) dataset

Prediction step

100

(b) KdV dataset

Figure 4: Identifiability of the PDE dynamics. (a) relationship between the learned potential and the real potential, and each point corresponds to one snapshot; (b) evolution of the learned potentials for the KdV dataset. For both the two figures, we use colors to distinguish 10 different trajectories in the test dataset. for elliptic dissipative systems with M and W independent of the state, the Onsager factorization is identifiable up to an affine transformation. Meanwhile, for the purely conservative KdV system, the learned dynamics strictly conserve the learned potential. This confirms our theoretical assertion that the framework cannot spuriously inject artificial dissipation into a genuinely conservative system. To sum up, the experimental results on the PDE cases have numerically confirmed the identifiability of the proposed generalized Onsager principle shown in Theorem 4.8. 6.2.3

Empirical validation of a priori discrete stability

In section 4, we have established rigorous a priori guarantees that apply universally to any dynamics within our hypothesis class. An important consequence of this framework, formalized in Theorem 4.10, is the preservation of discrete energy dissipation. Specifically, for elliptic dissipative systems, the explicit numerical integration strictly enforces a monotonically non-increasing free energy, guaranteeing trajectory stability provided the inference step size remains below a critical threshold. To empirically validate this theoretical bound, we evaluate the inference dynamics of both 1D and 2D Allen–Cahn models. As illustrated in fig. 5, the observed behavior aligns with our theoretical predictions. When the inference time step exceeds the characteristic scale used during training (3 × 10−3 compared to 1 × 10−3 ), the predicted trajectory diverges, and the strict monotonic decay of the learned potential is subsequently violated. Conversely, selecting an inference step smaller than 1 × 10−3 safely preserves the monotonic decay of the free energy, leading to robust, long-term trajectory stability without any degradation in relative error. These observations empirically confirm the discrete stability of our SpecOnsNet, consistent with the theoretical guarantee induced by the structurally constrained hypothesis class.

10 1

100

t = 10 5 t = 3 × 10 5 t = 10 4 t = 3 × 10 4 t = 10 3 (train) t = 3 × 10 3

V (u(t))

Relative RMSE

10 2

10 3

t = 10 5 t = 3 × 10 5 t = 10 4 t = 3 × 10 4 t = 10 3 (train) t = 3 × 10 3 0.00

0.02

0.04

0.06

Prediction lead time t

0.08

40 20 0 20 40 60

0.10

0.00

(a) prediction error along rollout (1D)

0.02

0.04

0.06

Prediction lead time t

t = 10 5 t = 3 × 10 5 t = 10 4 t = 3 × 10 4 t = 10 3 (train) t = 3 × 10 3

200

V (u(t))

Relative RMSE

400

10 2 t = 10 5 t = 3 × 10 5 t = 10 4 t = 3 × 10 4 t = 10 3 (train) t = 3 × 10 3 0.00

0.02

0.04

0.06

Prediction lead time t

0.08

0.10

(b) learned potential along rollout (1D)

10 1

10 3

0.08

0 200 400

0.10

0.00

0.02

0.04

0.06

Prediction lead time t

0.08

(d) learned potential along rollout (2D)

Figure 5: Inference step size scaling for the Allen–Cahn PDE models.

0.10

6.3

Learning mesoscopic dynamics from microscopic data

Whereas the continuum examples validate the framework against known reference dynamics, we next turn to a more practical multiscale setting: inferring effective mesoscopic spatio-temporal dynamics directly from microscopic simulation data without assuming a closed-form PDE a priori. In particular, we demonstrate our approach on nonlinear particle chains, a classical setting in which the mesoscopic energy structure reveals how local interactions give rise to energy transfer, nonlinear wave propagation, recurrence, and continuum-scale behavior. The Fermi–Pasta–Ulam–Tsingou (FPUT) chain [18] represents the classical setting in which a weak nonlinear perturbation of a harmonic interaction, under the long-wave and small-amplitude scaling, produces the KdV-type mesoscopic potential that balance nonlinearity and dispersion [41]. The polymer chain with Finitely Extensible Nonlinear Elastic (FENE) [38], by contrast, introduces a finite-extensibility constraint and the corresponding mesoscopic description is not available. Existing data-driven methods may reproduce coarse trajectories, but they do not identify the key mesoscopic potential that determines whether coarse-graining preserves or reshapes the KdV-type potential. In this subsection, we show that our method recovers this learned mesoscopic potential as a physically interpretable quantity, enabling a direct comparison between the established FPUT energy mechanism and the finite-extensibilityinduced effective potential of the polymer chain. 6.3.1

Microscopic chain models

We consider a one-dimensional nearest-neighbor chain with displacement qn and strain rn = qn+1 − qn . For an interaction potential V , the microscopic dynamics are q̈n = ∇V (rn ) − ∇V (rn−1 ),

yielding r̈n = ∇V (rn+1 ) − 2∇V (rn ) + ∇V (rn−1 ).

(49)

For the FPUT chain, V (r) = c2 r2 /2 + αr3 /3. In this case, the mesoscopic limit is analytically understood. Under the long-wave, small-amplitude scaling rn (t) = ε2 u(ξ, τ ),

ξ = ε(n − ct),

τ = ε3 t,

(50)

Taylor expansion and dominant balance give −2cuξτ =

c2 uξξξξ + α(u2 )ξξ . 12

(51)

With periodic or decaying boundary conditions, integration in ξ yields the KdV equation ut + akdv uux + bkdv uxxx = 0,

(52)

where akdv = α/c and bkdv = c/24. Figure 6 illustrates how microscopic particle-chain models are coarse-grained into mesoscopic variables, As a contrasting microscopic model, we use the FENE potential and force law HR2 r2 Hr V (r) = − log 1 − 2 , F (r) = ∇V (r) = . (53) 2 R 1 − (r/R)2 Here H and R are fixed. Unlike the polynomial FPUT force, this finite-extensibility law cannot be faithfully approximated by the Taylor expansion and therefore does not reduce to a simple standard continuum model under the same coarse-graining procedure eq. (50), making the FENE chain a more demanding test case for learning mesoscopic dynamics directly from microscopic simulations. First, we evaluate the prediction accuracy of our model compared with other data-driven approaches, and the quantitative performance is exhibited in table 1. The trends of superior performance of SpecOnsNet on the microscopic datasets mirror those of the continuum examples. Our model consistently yields the highest overall accuracy and stability among all tested methodologies, which provide a reliable learned dynamics for our further study on the interpretability.

trajectory #1 ( =1.07, R 2=0.992) trajectory #2 ( =1.13, R 2=0.978) trajectory #3 ( =1.00, R 2=0.999) trajectory #4 ( =1.07, R 2=0.997) trajectory #5 ( =1.27, R 2=0.988) time step t used for training: 10 2

medians |V (us + 1) V (us)|

Figure 6: Illustration of microscopic chain models and their coarse-graining.

10 1

10 3

5 × 10 3

Inference step size t

trajectory #1 ( =1.12, R 2=0.998) trajectory #2 ( =1.08, R 2=0.985) trajectory #3 ( =1.12, R 2=0.998) trajectory #4 ( =1.16, R 2=0.991) trajectory #5 ( =1.19, R 2=0.981) time step t used for training: 5 × 10 3

10 2

5 × 10 4

10 2

(a) FPUT chain

10 3

Inference step size t

5 × 10 3

(b) FENE chain

Figure 7: One-step variation of the learned potentials Vθ for the microscopic chain models. For each test trajectory we run a log-log linear fit; the resulting scaling exponent α and the the coefficients of determination R2 are appended in the legend. 6.3.2

Empirical validation of a priori theoretical analysis

Recall that in section 4 we have provided a priori theoretical analysis, which applies to any dynamics in our hypothesis class. In particular, we validate the discrete stability described in Theorem 4.10, which states that the one-step variation for the learned potential should be of order O(∆t2 ). To this end, we plot the variation under different time-step sizes (from 0.1× to 1.0× the step size used for training) for both the two microscopic models in fig. 7. We compute the median of the variation for robustness since for some rare cases, the predicted trajectory diverges because the network has never been trained with different step sizes. Empirical scaling indicates a variation of approximately O(∆t1.15 ). While this deviates from the theoretical quadratic law, the discrepancy is physically reasonable. By the one-step variation of the discrete energy eq. (34), if the network perfectly identifies the conservative dynamics, then only the quadratic term remains. However, in practice, the learned operators are approximations, resulting in a small but non-zero residual dissipation. The observed scaling of 1.15 thus reflects a superposition of a small O(∆t) contribution from the fitting error dominating the high-order term O(∆t2 ). 6.3.3

Exploring the potential evolution

As shown in fig. 8, the learned dynamics approximately conserve both the learned potential. These results indicate that the proposed framework can recover the conservative property of the underlying mesoscopic dynamics from microscopic trajectory data alone, without assuming any prior knowledge of the mesoscopic evolution equation. 6.3.4

Interpretability of the learned potentials

By revisiting our ansatz of V defined in eq. (38), intuitively we may want to know the learned component of the potential, or more specifically, the pointwise F (u) term and the V2 (u) designed for the residuals.

Learned potential V (u(t))

825

110 100 90 80 70

820 815 810 805 800 795 790 785

Prediction step

100

(a) FPUT chain

Prediction step

100

(b) FENE chain

Figure 8: Evolution of the learned potentials Vθ for the microscopic chain models. We use colors to distinguish different trajectories. Firstly, we evaluate the learned components V0 and V2 on a randomly sampled initial state u0 scaled by a varying amplitude a in fig. 9. Two critical observations emerge from this scaling. For both cases, V0 (au0 ) carries the dominant amplitude dependence, whereas V2 (au0 ) varies on a much smaller scale. Thus V2 represents a weaker correction beyond the pointwise constitutive ansatz, capturing residual effects such as nonlocal or higher-order interactions. Besides, the form of V0 is strongly dependent on the microscopic physics. For the FPUT chain, V0 (au0 ) displays a polynomial-type growth, consistent with the KdV-type Hamiltonian mechanism in which the local energy generates nonlinear transport. In contrast, the learned pointwise energy V0 for the FENE chain initially follows a similar behavior for small amplitudes but then departs markedly from the FPUT profile, exhibiting substantially stronger, super-polynomial growth. This contrast shows that the mesoscopic dynamics learned from the FENE chain is governed by a different nonlinear response once the finite-extensibility effect becomes active. We also plot the learned F (u) term in fig. 10 by varying the variable u within a range where most of the training data stay to avoid out-of-distribution testing. The learned profiles further confirm the distinction between the two microscopic chains. For the FPUT chain, F (u) is predominantly asymmetric, nearly odd over the sampled range, and displays super-quadratic growth, in qualitative agreement with the form of a KdV-type local Hamiltonian density. In contrast, the learned profile for the FENE chain is similar to the FPUT case only near the origin, but departs substantially as the amplitude increases. Taken together, these observations suggest that the two microscopic interaction laws induce different effective nonlinear constitutive relations at the mesoscopic level. In particular, the learned FENE dynamics are not well described by a KdV-type equation with merely modified coefficients, but instead encode a distinct effective nonlinear constitutive relation at the mesoscopic scale.

Conclusion

In this work, we propose a departure from standard pipelines of modelling dynamic processes via evolution equations. Rather than deriving a fixed, instance-wise effective equation and subsequently analyzing its properties, we introduce a framework centered on a highly expressive yet mathematically constrained hypothesis class. By lifting the generalized Onsager principle to infinite-dimensional Gelfand triples, we define a unified class encompassing both dissipative and conservative mesoscopic dynamics. Crucially, we establish theoretical guarantees uniformly across this entire class. This ensures that foundational physical properties, including global well-posedness, asymptotic stability, unique factorization identifiability, and discrete energy dissipation, are guaranteed a priori. Consequently, regardless of the specific data-driven model used to fit the dynamics, the resulting system inherently preserves these structural and analytical certainties.

FPUT

|V2(au0) V2(0)|

10 1

101

u0 #1 ( =2.37, R 2=1.000) u0 #2 ( =2.36, R 2=0.999) u0 #3 ( =2.34, R 2=1.000) u0 #4 ( =2.41, R 2=0.999) u0 #5 ( =2.38, R 2=0.999) fit window: a < 0.2

10 1

10 2 1 10

100

scalar a

100

scalar a FENE V2 scaling u0 #1 u0 #2 u0 #3 u0 #4 u0 #5

100

10 2 1 10

10 1

10 2 1 10

100

scalar a FENE V0 scaling

V2 scaling

u0 #1 u0 #2 u0 #3 u0 #4 u0 #5

100

10 2 1 10

|V0(au0) V0(0)|

V0 scaling

u0 #1 ( =2.51, R 2=0.999) u0 #2 ( =2.10, R 2=1.000) u0 #3 ( =1.98, R 2=1.000) u0 #4 ( =2.26, R 2=0.999) u0 #5 ( =2.46, R 2=0.999)

|V2(au0) V2(0)|

|V0(au0) V0(0)|

101

100

scalar a

Figure 9: Scaling property of the components V0 and V2 for the microscopic chain models. For each test u0 , we run a log-log linear fit of V0 w.r.t. the scalar a; the resulting scaling exponent α and the the coefficients of determination R2 are appended in the legend.

FPUT

pointwise density of V0

FENE

F(u) F(0) V0(u) V0(0) region where u is small

energy density (shifted)

5 0 5 10 15

1.00

0.75

0.50

0.25

0.00

0.25

0.50

0.75

0.2 0.1 0.0 0.1 0.2 0.3 1.00

1.00

pointwise density of V0

F(u) F(0) V0(u) V0(0) region where u is small

0.3

0.75

0.50

0.25

0.00

0.25

0.50

Figure 10: The pointwise component F (u) and V0 (u) of the learned potentials Vθ

0.75

1.00

Empirically, this structured approach has shown effectiveness across both continuum and particlebased settings with a specific implementation, Spectral OnsagerNet. For continuum PDEs such as the KdV and Allen–Cahn equations, Spectral OnsagerNet consistently improves predictive accuracy and long-time numerical stability over unconstrained baselines. Moreover, it faithfully recovers the correct underlying physical regimes, explicitly capturing the exact affine transformations of dissipative free-energy functionals and the conservation of purely Hamiltonian potentials. Furthermore, when learning mesoscopic dynamics directly from microscopic particle chains, the formulation of Spectral OnsagerNet provides vital interpretable diagnostics. By analyzing the amplitude dependence of the learned energy components, we directly observe how microscopic interaction laws are reflected in the effective nonlinear constitutive relation at the mesoscopic scale. This interpretability allows the framework to distinguish KdV-type Hamiltonian behavior from mesoscopic dynamics governed by finite extensibility in the FENE chain. While the proposed approach is a first step towards data-driven construction of mesoscopic descriptions, important theoretical and computational challenges remain. The current reliance on convolution structures to diagonalize the operators M and W restricts the methodology to periodic boundary conditions. Extending this uniform class-level analysis to handle complex spatial geometries, diverse boundary conditions, and fully state-dependent transport operators remains a critical next step. More fundamentally, while the proposed hypothesis class is highly expressive, it does not universally encompass all non-equilibrium mesoscopic phenomena. For complex multiscale systems governed by fundamentally different thermodynamic or kinetic rules, it will be necessary to propose and integrate alternative structure-preserving physical principles within this hypothesis-driven framework. Ultimately, we envision this approach as a rough blueprint for the automated discovery of reliable, physically consistent, and analytically sound multiscale models across a diverse spectrum of governing laws.

Data availability All the codes for data generation, network training, and visualization are publicly available on the GitHub repository: https://github.com/MLDS-NUS/Meso-SpecOnsNet.

Acknowledgments The computational work for this article was partially performed on resources of the National Supercomputing Centre, Singapore (https://www.nscc.sg). During the preparation of this work, the authors used large language models (LLMs) to assist with the code writing and to polish the written text for spelling and grammar. The authors thoroughly reviewed and edited the content and take full responsibility for the final publication.

Experiment details

A.1

Baselines

Classical solver with true operators We use the same equation as that of the reference trajectory to evolve the solution numerically. A larger step size is set to keep consistent with the data-driven approaches. Meanwhile, we adopt the SBDF1 scheme (see section D.3) as a one-step variant of the reference scheme for a fair comparison. This model serves as a baseline for one-step classical methods with a large step size. Fourier Neural Operator (FNO) We use the FNO architecture [28] to directly learn the dynamics in the physical space, which serves as a robust baseline in the field of operator learning. More specifically, we set us+1 ≈ us + γθ (us )∆t, (54)

where the network γθ is implemented with the FNO architecture. OnsagerNet We use the OnsagerNet architecture [40] originally designed for learning Onsager principle. Since the reference snapshots are saved as grid data, we can restate the learning task in the discrete space. Concretely, we treat the Nx -point PDE grid directly as an Nx -dimensional ODE, evolving under du (55) = − Mθ (u) + Wθ (u) ∇Vθ (u), u ∈ RNx . dt The matrix M (u) and W (u) are parameterized as low-rank matrices for computational efficiency, and the potential is implemented as an MLP. Spectral OnsagerNet (SpecOnsNet, our work) We apply the parameterization as described in eqs. (46) and (47). The detailed workflow is displayed in fig. 2.

A.2

Evaluation metric

To quantify the prediction performance, we introduce the relative RMSE at time t, evaluated as v uP 2 u ∥upred (t)(xi ) − uref (t)(xi )∥2 Relative RMSE(t) = t i , (56) P 2 i ∥uref (t)(xi )∥2 where uref (t) and upred (t) are reference and prediction states for time t, respectively.

A.3

Training cost

The network training is performed on a single NVIDIA A100 GPU based on the PyTorch framework, and each experimental configuration requires about 2 hours. Details of the optimizer and learningrate scheduler for each experiment are given in section G.

Examples

In this section, we exhibit some classical PDE dynamics covered by the generalized Onsager principle, or more specifically, those satisfying Hypothesis 3.2, including pure gradient flows (W ≡ 0), purely Hamiltonian flows (M ≡ 0), and mixed dynamics such as the incompressible Navier–Stokes equation, where both M and W do not vanish. Example B.1 ((Hyper/Fractional) diffusion). Let the state space be X = L2 (Ω). Consider the quadratic energy V : X → R: Z 1 2 δV V (u) = u dx, = u. (57) δu Ω 2 This yields the (hyper/fractional) diffusion dynamics: ∂t u = −(−∆)s u = − (−∆)s + |{z} 0 |{z} u . | {z } M(u)

W(u)

(58)

δV δu

Note that for s = 1, this recovers the classical Heat equation. Example B.2 (Reaction–Diffusion). Let the state space be X = H 1 (Ω). Consider the Ginzburg– Landau free energy V : X → R: Z 2 δV ε V (u) = |∇u|2 + F (u) dx, = −ε2 ∆u + F ′ (u), (59) 2 δu Ω where F (u) is a double-well potential. This energy generates two distinct dynamics depending on the choice of operators: 22

1. Allen–Cahn (Non-conserved): By choosing M = I, we obtain the L2 -gradient flow: i h ∂t u = ∆u − F ′ (u) = − |{z} −∆u + F ′ (u) . I + |{z} 0 {z } | M(u)

W(u)

δV δu

2. Cahn–Hilliard (Conserved): By choosing M = −∆, we obtain the H −1 -gradient flow: i h 2 ∂t u = ∆(−ε2 ∆u + F ′ (u)) = − −∆ + |{z} −ε ∆u + F ′ (u) . 0 |{z} {z } | W(u)

M(u)

(60)

(61)

δV δu

Example B.3 (KdV). Let the spatial dimension d = 1. Consider the energy V : X → R defined as Z 1 2 δV V (u) = ux − u3 dx, = −uxx − 3u2 . (62) 2 δu R We have the KdV dynamics: ∂t u = −uxxx − 6uux = −[ |{z} 0 + (−∂x )](−uxx − 3u2 ). {z } | {z } | M(u)

(63)

δV δu

W(u)

Example B.4 (nonlinear Schrödinger). Let u = (p, q) ∈ L2 (Ω; R2 ) and write ψ = p + iq. By setting 0 −1 M(u) ≡ 0, W(u) ≡ , (64) 1 0 and

Z V (u) = Ω

1 λ |∇u|2 + |u|4 dx, 2 4

(65)

we have the nonlinear Schrödinger dynamics: iψt = −∆ψ + λ|ψ|2 ψ.

(66)

Example B.5 (1D conservation laws). Let the spatial dimension d = 1. Consider the energy V : X → R defined as Z δV = G′ (u). (67) V (u) = G(u) dx, δu R We have the 1D conservation dynamics: ∂t u = −G′ (x)x = −( |{z} 0 + ∂x ) G′ (u) . |{z} | {z } M(u)

W(u)

(68)

δV δu

Example B.6 (Incompressible Navier–Stokes). Let Ω ⊂ Rd and X be the L2 -closure of smooth divergence-free vector fields with suitable boundary conditions. Let Z 1 δV V (u) = |u|2 dx, = u ∈ X. (69) 2 Ω δu Let P be the Leray projector onto divergence-free fields, and define M(u)v = νP ∆v,

W(u)v = −P [(u · ∇)v] ,

(70)

then we have the incompressible Navier–Stokes dynamics: ∂t u + (u · ∇)u = −∇p + ν∆u, with the pressure enforced by P .

∇·u=0

(71)

Dynamics diffusion Allen-Cahn Cahn-Hilliard KdV Schrödinger 1D conservation Navier–Stokes

differential equation

0 0 0

ut + (−∆)s u = 0 ∂t u = ∆u − F ′ (u) ∂t u = ∆[−ε2 u + F ′ (u)]

G(u)

0 0 0

−∂x i ∂x

ut + uxxx + 6uux = 0 iψt = −∆ψ + λ|ψ|2 ψ ut + G′ (u)x = 0

1 2 2 |u|

νP ∆

−P [(u · ∇)·]

ut + (u · ∇)u = −∇p + ν∆u

V (u) =

Ω 1 2 u 2 2

· dx

1 2 |∇u| + F (u) ε2 2 2 |∇u| + F (u) 1 2 3 2 ux − u 1 λ 2 4 2 |∇u| + 4 |u|

M (−∆) I −∆

Table 2: Some classical PDE dynamics that fit our setup

Detailed proofs

C.1

Proof of Theorem 4.1

Let Pn be the canonical projection from H = L2 (Ω) to Hn defined as X Pn : u 7→ (u, ek ) ek ,

(72)

|k|≤n

where we set Hn as the finite-dimensional space spanned by the Fourier modes ek for |k| ≤ n. Consider the time-dependent ODE system ∂t un (t) = Pn G(un (t)),

un (0) = Pn uinit ,

(73)

which admits a local solution un on some interval [0, Tn ) by the Peano existence theorem since the subspace Hn is finite-dimensional and the right-hand side is continuous. Next, we show that the functional V decays for any solution un (t). For simplicity, we suppress the dependency on u in the operators. For each n, we have X c[k] ⟨µ, ek ⟩ ∗ ek (z) Pn (Mµ)(z) = (Mµ, ek ) ek (z) = M (74) X ,X |k|≤n

by Theorem 3.6. Meanwhile, defining Pn∗ : X ∗ → X ∗ as the adjoint of Pn via ⟨Pn∗ µ, v⟩X ∗ ,X = ⟨µ, Pn v⟩X ∗ ,X ,

(75)

we have * M (Pn∗ µ) (z) = ⟨Pn∗ µ, KM (z − ·)⟩X ∗ ,X =

+ µ,

(KM (z − ·), ek ) ek

|k|≤n

where

(76)

X ∗ ,X

Z (KM (z − ·), ek ) =

KM (z − y)ek (y) dy Z c[−k] . = ek (z) KM (z − y)ek (z − y) dy = ek (z)M Ω

(77)

Ω

Hence M (Pn∗ µ) (z) =

c[−k] ⟨µ, ek ⟩ ∗ ek (z), M X ,X

(78)

|k|≤n

which is identical with eq. (74) since e−k (z) = ek (−z) = ek (z). The same argument applies to W, so Pn Mµ = MPn∗ µ, Pn Wµ = WPn∗ µ. (79) 24

It follows that d V (un (t)) = − ⟨µ(un ), Pn [M(un ) + W(un )]µ(un )⟩X ∗ ,X dt = − ⟨Pn∗ µ(un ), Pn [M(un ) + W(un )]µ(un )⟩X ∗ ,X

(80)

= − ⟨Pn∗ µ(un ), [M(un ) + W(un )]Pn∗ µ(un )⟩X ∗ ,X ≤ 0. Consequently, the energy dissipation leads to V (un (t)) ≤ V (un (0)) = V (Pn (uinit )),

∀t ∈ [0, Tn ).

(81)

Note that Pn (uinit ) → uinit as n → ∞, so the right-hand side is bounded by a constant independent of n from above. By the coercivity of V , ∥un (t)∥X ≤ C1 ,

∀n,

∀t ∈ [0, Tn ).

(82)

Hence each solution un cannot blow up in finite time, and we may extend it to [0, +∞). Furthermore, ∥∂t un (t)∥X ≤ ∥G(un (t))∥X ≤ C2 ,

∀n,

∀t ∈ [0, +∞),

(83)

since G is Lipschitz continuous resulting from Hypothesis 3.5. Till now, we have shown the uniform boundness for un and ∂t un . According to Hypothesis 3.1, the embedding X ,→ L2 (Ω) is compact and L2 (Ω) ,→ X ∗ is continuous, sothe space W̃ (0, T ; X ) := u ∈ L2 (0, T ; X ) | ∂t u ∈ L2 (0, T ; X ∗ ) (84) is compactly embedded in L2 (0, T ; L2 (Ω)) for any fixed T > 0 by the Aubin–Lions lemma. Let (0) sequence Tm → +∞ as m → +∞. For each m > 0, un = un and take an arbitrary n increasing o (m−1) we can find a subsequence of un with a weak limit in L2 (0, Tm ; X ) by the reflexivity of n

L2 (0, Tom ; X )nand theoboundness according to eqs. (82) and (83). Recursively we find a subsequence n (m) (m−1) un ⊆ un such that n

u(m) ⇀ u(m) ∈ L2 (0, Tm ; X ), n ∂t u(m) ⇀ v (m) ∈ L2 (0, Tm ; X ), n

(85)

u(m) → u(m) ∈ L2 (0, Tm ; L2 (Ω)). n By the uniqueness of strong limits, u(m+1) coincides with u(m) on L2 (0, Tm ; L2 (Ω)) for any m > (m) 0, so by induction there exists a unique limit u defined on t ∈ [0, +∞) such that un → u in 2 2 L (0, Tm ; L (Ω)) for any m > 0. Meanwhile, for any m > 0 and any sufficiently smooth test function ϕ, Z Tm D Z Tm D E E (m) ϕ, ∂t u(m) dt = − ∂ ϕ, u dt. (86) t n n X ∗ ,X

By taking the limit n → ∞, we have that Z Tm D E ϕ, v (m) 0

X ∗ ,X

Z Tm dt = − 0

⟨∂t ϕ, u⟩X ∗ ,X dt,

(87)

which indicates that ∂t u = v (m) in L2 (0, Tm ; X ). To sum up, without loss of generality, we may assume that there exists u defined on t ∈ [0, +∞) such that for any finite time T > 0, un ⇀ u ∈ L2 (0, T ; X ),

∂t un ⇀ ∂t u ∈ L2 (0, T ; X ), and un → u ∈ L2 (0, T ; H).

(88)

Finally, we prove that the limit function u satisfies eq. (4), and we need to handle the convergence of the nonlinear operator G(un ). Consider Z T Xn (w) = (G(un ) − G(w), un − w) dt (89) 0

for an arbitrary w ∈ L2 (0, T ; X ). The term Xn (w) can be decomposed as Z T

Z T (G(un ), un ) dt − {z }|

Xn (w) = 0

Z T (G(un ), w) dt − (G(w), un − w) dt . 0 {z }

(90)

Zn (w)

The Zn term converges as Z T Zn (w) = −

Z T (∂t un , w) dt −

Z T →−

(G(w), un − w) dt

0 Z T

(∂t u, w) dt − 0

(91) (G(w), u − w) dt

as n → ∞ due to the weak convergence of un and ∂t un . Meanwhile, for the Yn term we have Z T 1 1 (∂t un , un ) dt = ∥un (T )∥22 − ∥un (0)∥22 . Yn = (92) 2 2 0 For any test function ϕ ∈ X , Z t ⟨ϕ, un (T )⟩X ∗ ,X − ⟨ϕ, un (0)⟩X ∗ ,X =

⟨ϕ, ∂t un ⟩X ∗ ,X dt,

(93)

which converges to Z t 0

⟨ϕ, ∂t u⟩X ∗ ,X dt = ⟨ϕ, u(T )⟩X ∗ ,X − ⟨ϕ, u(0)⟩X ∗ ,X

(94)

as n → ∞. Recall that by definition, un (0) = Pn (u(0)) → u(0)

in L2 (Ω),

(95)

so we have ⟨ϕ, un (T )⟩X ∗ ,X → ⟨ϕ, u(T )⟩X ∗ ,X

(96)

as n → ∞, and thus un (T ) converges weakly to u(T ) as well. As a result, we have lim inf Xn (w) = lim inf Yn + lim Zn (w) n→∞

n→∞

1 1 ≥ ∥u(T )∥22 − ∥u(0)∥22 + lim Zn (w) n→∞ 2 2 Z T Z T Z T = (∂t u, u) dt − (∂t u, w) dt − (G(w), u − w) dt 0

(97)

Z T (∂t u − G(w), u − w) dt.

= 0

Besides, by the Lipschitz continuity of G in Hypothesis 3.5, Z T lim sup Xn (w) ≤ lim sup ∥G(un ) − G(w)∥2 ∥un − w∥2 dt n→∞

n→∞

Z T ≤ L3 (R) 0

(98) lim sup ∥un − w∥22 dt = L3 (R)∥u − w∥2L2 (0,T ;H) n→∞

if w(t) ∈ BX (R) for all t, where we take R = 1 + C1 . We choose w = u + λψ for a smooth test function ψ and a sufficiently small scalar λ > 0, then Z T (∂t u − G(u + λψ), ψ) dt ≤ lim inf Xn ≤ λL3 (R)∥ψ∥2L2 (0,T ;H) . (99) 0

n→∞

Taking the limit λ → 0+ , we obtain Z T (∂t u − G(u), ψ) dt ≤ 0,

∀ψ ∈ X .

(100)

Since the inequality holds for ±ψ, the integral vanishes, which implies that ∂t u and G(u) coincide a.e. on L2 (0, T ; X ).

C.2

Proof of Lemma 4.2

Since 12 ∥ · ∥2X is Fréchet differentiable by Hypothesis 3.1, the limit J(u)(v) =

1 lim t−1 (∥u + tv∥2X − ∥u∥2X ) 2 t→0

(101)

exists for all v ∈ X . By the triangle inequality, ∥u + tv∥X ≤ ∥u∥X + |t|∥v∥X , so |J(u)(v)| ≤

1 lim sup |t|−1 (∥u∥X + |t|∥v∥X )2 − ∥u∥2X = ∥u∥X ∥v∥X , 2 t→0

(102)

which gives ∥J(u)∥X ∗ ≤ ∥u∥X . Meanwhile, J(u)(u) =

∥u + tu∥2X − ∥u∥2X 1 (1 + t)2 − 1 1 lim = lim ∥u∥2X = ∥u∥2X , 2 t→0 t 2 t→0 t

(103)

and therefore ∥J(u)∥X ∗ = ∥u∥X .

C.3

Proof of Theorem 4.5

According to Proposition 3.3, any trajectory {u(t)}t∈[0,+∞) is uniformly bounded in X . By the Rellich–Kondrachov Theorem, the embedding X = H 1 (Ω) ,→ L2 (Ω) is compact. Therefore, we have a non-empty ω-limit set ω(u(0)). We show that ω(u(0)) is positively invariant under the flow S(t), and the energy V is a constant on the ω-limit set. Consider any sequence tn → ∞ such that u(tn ) → v ∈ ω(u(0)). For any τ > 0, S(τ )v = S(τ ) lim u(tn ) = lim S(τ )u(tn ) = lim u(tn + τ ), (104) n→∞

n→∞

which implies S(τ )v ∈ ω(u0 ). If V is not a constant on the ω-limit set, we may assume that u(tnk ) → v1 and u(tmk ) → v2 in H as k → +∞ and V (v1 ) < V (v2 ). By the definition of convergence, there exists a K large enough such that V (u(tnK )) <

1 1 V (v1 ) + V (v2 ) < V (v2 ). 2 2

(105)

However, since V is non-increasing on the trajectory, the sequence V (u(tmk )) will eventually be smaller than V (u(tnK )), and thus the corresponding limit V (v2 ) is not greater than V (u(tnK )), which is a contradiction. Consider a trajectory lying on the ω-limit set. By taking the time derivative of V it is easy to verify that the free energy dissipation vanishes, and consequently we have ω(u(0)) ⊆ Z. Furthermore, if the elliptic condition eq. (27) holds, then for any φ ∈ Z, we have that ⟨µ(φ), M(φ)µ(φ)⟩X ∗ ,X = 0 implies µ(φ) = 0, so φ ∈ E.

C.4

Proof of Lemma C.1

(106)

C.5

Proof of Theorem 4.8

Lemma C.1. Consider a Fréchet differentiable functional V : X → R. The Fourier transform of the Fréchet derivative µ = δV /δu can be written as µ̂[k] = ∂−k v = ∂k v = µ̂[−k] ,

(107)

where we use the notation ∂k v = ∂v/∂ û[k] for any k ∈ Z. Proof. By definition we have d ∂v δV d = = = ∂−k v. µ̂[k] = , ek = V (u + τ ek ) v(û + τ ed −k ) [−k] δu dτ dτ ∂ û ∗ τ =0 τ =0 X ,X

(108)

Furthermore,

δV , ek + e−k ∈ R, δu X ∗ ,X δV [k] [−k] , ie−k − iek ∈ R, i∂k v − i∂−k v = iµ̂ − iµ̂ = δu X ∗ ,X ∂k v + ∂−k v = µ̂[k] + µ̂[−k] =

(109)

which gives µ̂[−k] = µ̂[k] and ∂−k v = ∂k v for all k. Lemma C.2. Let (M1 , W1 ) and (M2 , W2 ) be independent of u and satisfy (M1 + W1 )µ1 (u) = (M2 + W2 )µ2 (u),

µ1 =

δV1 , δu

µ2 =

δV2 δu

(110)

with Hypotheses 3.2 and 3.4. Assume that V2 is irreducible on the spectral space. Define n o c[k] + W c[k] ̸= 0 , Z1 = k ∈ Z | M 1 1

then

λk :=

c[k] + W c[k] M 2 2 c[k] + W c[k] M 1

(111)

[k]

is a constant λ for k ∈ Z1 . Moreover, either µ̂2 ≡ 0 for all k ∈ Z1 or λ ∈ R holds. Proof. By transforming eq. (110) from the physical space to the spectral space, we have c[k] + W c[k] µ̂[k] , ∀k ∈ Z. c[k] + W c[k] µ̂[k] = M M 2 2 2 1 1 1

(112)

By the definition of λk we have [k]

[k]

µ̂1 = λk µ̂2 ,

∀k ∈ Z1 .

(113)

[k] Since that µ̂i = ∂vi /∂ û[k] for i ∈ {1, 2} by Lemma C.1, it follows that [k]

[l]

∂ µ̂i ∂ µ̂i = , [l] ∂ û ∂ û[k]

i ∈ {1, 2}

(114)

for any indices k and l. Then, plugging eq. (113) into the above equation gives [k]

[l]

∂λk [k] ∂ µ̂ ∂λl [l] ∂ µ̂2 µ̂ + λk 2[l] = µ̂ + λl [k] , ∂ û[l] 2 ∂ û ∂ û[k] 2 ∂ û

(115)

which implies ∂ 2 v2 =0 (116) ∂ û[k] ∂ û[l] as the derivatives of λk and λl vanish. It follows immediately by the irreducibility of V2 that [k ] λk = λ ∈ C for any k ∈ Z1 . Suppose that µ̂2 0 ̸≡ 0 for some k0 ∈ Z1 . It is clear that −k0 ∈ Z1 by Lemma C.1, then by substituting k with ±k0 in eq. (113), we have (λk − λl )

[k ]

[−k0 ]

λk0 µ̂2 0 = µ̂1 0 = µ̂1

[−k0 ]

= λ̄−k0 µ̂2

which indicates that λ = λk0 = λ̄−k0 = λ̄ as a real number. 28

[k ]

= λ̄−k0 µ̂2 0 ,

(117)

Proposition C.3. With the same assumptions as in Lemma C.2, if M1 is positive definite and V2 is not a constant functional, then there exists a constant λ > 0 such that M2 = λM1 ,

W2 = λW1 ,

µ1 = λµ2 .

(118)

Proof. It follows by the positive definiteness of M1 that Z1 = Z. By Lemma C.2, there exists a real constant λ such that c[k] + W c[k] = λM c[k] + λW c[k] , M 2 2 1 1

[k]

µ1 = λµ2 ,

∀k ∈ Z.

(119)

Due to the positive (semi-)definiteness of M1 and M2 , the positiveness of λ follows, and the proof is completed by separating the real and imaginary parts and taking the inverse Fourier transform. [k]

Proposition C.4. With the same assumptions as in Lemma C.2, if M2 = 0 and µ̂2 ̸≡ 0 for all k ∈ Z, then M1 µ1 (u) ≡ 0, W2 = λ′ W1 (120) for a constant λ′ ∈ R. Furthermore, we have µ1 = λ′ µ2 if ker W1 = 0. c[k0 ] ̸= 0 for some k0 ∈ Z. It follows that k0 ∈ Z1 by definition, then we have Proof. Suppose that M 1 c[k0 ] + W c[k0 ] = λM c[k0 ] + λW c[k0 ] , M 2 2 1 1

[k ]

µ1 0 = λµ2 0

(121)

for some λ ∈ R according to Lemma C.2. Comparing the real parts on both sides leads to λ = 0, [k ] c[k] µ[k] ≡ 0 for all k ∈ Z. Back to the which implies µ̂1 0 ≡ 0. Therefore, it can be concluded that M 1 1 spectral representation of eq. (110), we have c[k] µ̂[k] = W c[k] µ̂[k] , ∀k ∈ Z. W (122) 1 1 2 2 n o c[k] ̸= 0 and let λ′ := W c[k] /W c[k] ∈ R for any k ∈ Z ′ . Analogously, we may define Z1′ = k ∈ Z | W 1 1 2 1 k With the same deduction in Lemma C.2, we can show that λ′k remains a constant for k ∈ Z1′ because [k] [k] of the irreducibility of V2 . Hence µ̂1 = λ′ µ̂2 with a constant λ′ ∈ R for any k ∈ Z1′ . Meanwhile, c[k] for any k c[k] = 0 implies W c[k] = 0 according to eq. (122). As a result, W c[k] = λ′ W for k ∈ / Z1′ , W 1 1 2 2 ′ ′ ′ and thus W2 = λ W1 . Additionally, ker W1 = 0 implies Z1 = Z and thus µ1 = λ µ2 .

C.6

Proof of Theorem 4.10

By the Lipschitz continuity of µ1 , applying the Descent Lemma to the step from us to us+1 , we have that Z 1

V (us+1 ) − V (us ) =

µ(us + τ (us+1 − us )), us+1 − us X ∗ ,X dτ

Z 1 = 0

µ(us + τ (us+1 − us )) − µ(us ), us+1 − us X ∗ ,X dτ

+ ⟨µ(us ), us+1 − us ⟩X ∗ ,X Z 1 ≤ ⟨µ(us ), us+1 − us ⟩X ∗ ,X + L1 (R) τ ∥us+1 − us ∥2X dτ

(123)

L1 (R) s+1 2 ≤ µs , us+1 − us X ∗ ,X + u − us X 2 as long as max ∥us ∥X , us+1 X ≤ R.

(124)

Combining with the discrete update eq. (31), we may conclude that it suffices to ensure that 2

∥[M(us ) + W(us )]µ(us )∥X L1 (R)∆t ≤ 2 ⟨µ(us ), M(us )µ(us )⟩X ∗ ,X , 29

(125)

for some R satisfying eq. (124). Note that in eq. (124), R is partially determined by ∥un+1 ∥X , and we have to dismiss the dependency. Now, we claim that either µs = 0, or ∆t ≤

2 ⟨µs , Ms µs ⟩X ∗ ,X 2

∥(Ms + W s )µs ∥X L1 (Rs )

Rs = D +

2 ⟨µs , Ms µs ⟩X ∗ ,X ∥(Ms + W s )µs ∥X L1 (D)

(126)

for the sth step will ensure that V (us+1 ) ≤ V (us ) and ∥us ∥X ≤ D for all s, and we will prove it by induction. The constant D is chosen such that V (w) ≤ V (u0 ) =⇒ ∥w∥X ≤ D

(127)

for any w ∈ X , and the existence is guaranteed by the coercivity of V . When µs = 0, the claim is trivial. Otherwise, by triangle inequality, we can decompose the update as us+1 X ≤ ∥us ∥X + ∆t ∥(Ms + W s )µs ∥X 2 ⟨µs , Ms µs ⟩X ∗ ,X

≤ ∥us ∥X +

(128)

∥(Ms + W s )µs ∥X L1 (Rs ) 2 ⟨µs , Ms µs ⟩X ∗ ,X ≤D+ = Rs , ∥(Ms + W s )µs ∥X L1 (D)

which means that R = Rs makes eq. (123) hold. It follows that L1 (Rs ) s+1 2 u − us X 2 L1 (Rs ) 2 s s s = −∆t ⟨µ , M µ ⟩X ∗ ,X + ∥(Ms + W s )µs ∥X (∆t)2 2 ≤ 0.

V (us+1 ) − V (us ) ≤ µs , us+1 − us X ∗ ,X +

(129)

Hence we obtain the non-increasing behavior of V (us ), and by induction we have V (us+1 ) ≤ V (us ) ≤ · · · ≤ V (u0 ), which implies that us+1 X ≤ D as well. After proving the claim, by eq. (126) we have Rs ≤ D +

2 ⟨µs , (Ms + W s )µs ⟩X ∗ ,X

≤D+

∥(Ms + W s )µs ∥X L1 (D) 2(∥µ(0)∥X + L1 (∥us ∥X ) ∥us ∥X ≤D+ L1 (D)

2 ∥µs ∥X L1 (D) (130)

≤ 3D + 2 ∥µ(0)∥X L1 (D)−1 =: R0 . Therefore, the upper bound for the time step 2 ⟨µs , Ms µs ⟩X ∗ ,X 2

∥(Ms + W s )µs ∥X L1 (Rs )

= 2κL1 (Rs )−1 ≥ 2κL1 (R0 )−1

(131)

Moreover, by Hypothesis 3.5 we have κ=

⟨µs , Ms µs ⟩X ∗ ,X 2 ∥(Ms + W s )µs ∥X

≥

c 2 ∥Ms + W s ∥

≥

[∥M(0) + W(0)∥ + L2 (∥us ∥X ) ∥us ∥X ] c ≥ 2 =: κ0 [∥M(0) + W(0)∥ + L2 (D)D]

is always positive for elliptic dynamics when eq. (27) holds. 30

(132)

For conservative systems where M = 0, a similar deduction as eq. (123) gives V (us+1 ) − V (us ) ≤

1 1 2 2 L1 (R) us+1 − us X ≤ L1 (R)(∆t)2 ∥W s µs ∥X . 2 2

(133)

whenever eq. (124) holds. By assuming that the discrete trajectory does not blow up in finite time T , we set R̃ as the maximal X -norm on the trajectory, and the term ∥W s µs ∥X can be bounded by R̃ as well due to their Lipschitz continuity given by Hypothesis 3.5. The boundness of the potential changes across the whole trajectory follows immediately by decomposing the increment into multiple steps.

PDE models

We describe the detailed setup for the PDE models involved in this work, including the Allen–Cahn and KdV dynamics.

D.1

Allen–Cahn dataset

We consider the following Allen–Cahn dynamics ∂t u = ∆u + ε−2 (u − u3 ) = − |{z} I + |{z} 0 M(u)

W(u)

−∆u − ε−2 (u − u3 ) {z } |

(134)

δV δu

with ε = 0.1. The free energy functional V can be defined as Z 1 2u2 − u4 VAC = |∇u|2 − dx. 2 4ε2 Ω

D.2

(135)

KdV dataset

We consider the following KdV dynamics ∂t u = −uxxx − 6uux = − |{z} 0 + (−∂x ) | {z } M(u)

W(u)

i −uxx − 3u2 . | {z }

The free energy functional V can be defined as Z 1 2 VKdV = ux − u3 dx. 2 Ω

D.3

(136)

δV δu

(137)

Numerical steppers

Consider solving PDE on the spectral domain h i ∂t ût = − M̂(ût ) + Ŵ(ût ) µ̂t .

(138)

Suppose that the term µ̂t can be decomposed as a linear function plus a nonlinear function of ût , formulated as µ̂t = L̂ût + N̂ (ût ), (139) then we may use the time discretization h ih i ût+∆t − ût = − M̂(ût ) + Ŵ(ût ) L̂ût+∆t + N̂ (ût ) ∆t,

(140)

which is equivalent to h i ût − M̂(ût ) + Ŵ(ût ) N̂ (ût )∆t h i ût+∆t = . 1 + M̂(ût ) + Ŵ(ût ) L̂∆t

(141)

Such an update is usually referred to as the semi-implicit backward differentiation stepper of order 1 (SBDF1). Analogous to the leap-frog scheme, we may also use the time discretization h ih i ût+∆t − ût−∆t = −2 M̂(ût ) + Ŵ(ût ) L̂ût+∆t + N̂ (ût ) ∆t, (142) and the resulted scheme is abbreviated as SBDF2.

Microscopic chain models

Consider a one-dimensional chain where the displacement of the nth node is denoted by qn . The strain between the nth and the (n + 1)th nodes is defined as rn = qn+1 − qn .

E.1

(143)

FPUT chain model

The Fermi–Pasta–Ulam–Tsingou (FPUT) chain model employs a truncated quadratic force law given by F (r) = c2 r + αr2 . (144) Consequently, the equations of motion for the node displacements can be written as q̈n = F (rn ) − F (rn−1 ). This system can be cast in Hamiltonian form, with the Hamiltonian defined as X 1 vn2 + V (rn ) , F (r) = V ′ (r), H= 2 n yielding the canonical equations  ∂H   q̇n = ∂v = vn , n  ∂H  v̇n = − = F (rn ) − F (rn−1 ). ∂qn

(145)

(146)

(147)

Taking the second time derivative of the strain yields the governing equation for rn 2 2 r̈n = q̈n+1 − q̈n = c2 (rn+1 − 2rn + rn−1 ) + α(rn+1 − 2rn2 + rn−1 ).

(148)

To capture the continuum limit, we adopt a long-wave, small-amplitude ansatz describing a right-moving wave rn (t) = ε2 u(ξ, τ ), ξ = ε(n − ct), τ = εp t, (149) where the temporal scaling exponent p > 1 will be determined via dominant balance. Applying the chain rule, the differential time operator becomes ∂t = −cε∂ξ + εp ∂τ .

(150)

Consequently, the time derivatives of the strain are ṙn = −cε3 uξ + εp+2 uτ 32

(151)

and r̈n = ε2 c2 ε2 uξξ − 2cεp+1 uξτ + ε2p uτ τ = c2 ε4 uξξ − 2cεp+3 uξτ + O(ε2p+2 ).

(152)

Expanding the spatial shift terms rn±1 = ε u(ξ ± ε, τ ) in ε using a Taylor series gives rn+1 − 2rn + rn−1 = ε4 uξξ +

ε6 uξξξξ + O(ε8 ), 12

(153)

and 2 2 rn+1 − 2rn2 + rn−1 = ε6 (u2 )ξξ + O(ε8 ).

Substituting these expansions back into eq. (148) yields 2 c r̈n = c2 ε4 uξξ + ε6 uξξξξ + α(u2 )ξξ + O(ε8 ). 12

(154)

(155)

Equating the terms in eq. (152) and eq. (155), dominant balance dictates a cubic slow-time scale (p = 3). Collecting the O(ε6 ) terms leaves −2cuξτ =

c2 uξξξξ + α(u2 )ξξ . 12

(156)

Assuming periodic or rapidly decaying boundary conditions, integrating with respect to ξ yields uτ +

α 2 c (u )ξ + uξξξ = 0, 2c 24

(157)

which is the standard form of the Korteweg–de Vries (KdV) equation ut + akdv uux + bkdv uxxx = 0,

(158)

with akdv = α/c and bkdv = c/24.

E.2

FENE chain model

Alternatively, the Finitely Extensible Nonlinear Elastic (FENE) model utilizes the potential V (r) = −

HR2 log(1 − (r/R)2 ), 2

(159)

Hr , 1 − (r/R)2

(160)

which yields the force law F (r) =

where the parameters H and R are fixed. Because the macroscopic continuum limit of this specific force law does not readily reduce to a known, standard PDE, it provides a nontrivial testbed for evaluating the capabilities of our framework.

E.3

Numerical simulation

In this section, we detail the numerical simulation of the microscopic chain. We illustrate the procedure using the chain-KdV model, and an analogous approach applies to the FENE model. Recall that the microscopic chain operates in laboratory time t, whereas the macroscopic KdV equation evolves in slow time τ = ε3 t for a fixed scaling factor 0 < ε ≪ 1. Given an initial macroscopic profile u(ξ, 0), the corresponding initial conditions for the microscopic strain and strain rate are rn (0) = ε2 u(εn, 0), ṙn (0) = −cε3 uξ (εn, 0). (161) We integrate the system in time using the velocity Verlet algorithm. Introducing a micro timestep δt such that the macroscopic saving interval is ∆t = mδt for m substeps, the update rules for a single substep are as follows. 33

1/2

1. Compute the half-step velocity ṙn the current strain rn .

= ṙn + 21 r̈n δt, where r̈n is evaluated via eq. (148) using

1/2

2. Update the strain rn+ = rn + ṙn δt. 1/2

3. Compute the full-step velocity ṙn+ = ṙn the updated strain rn+ .

+ 12 r̈n+ δt, where r̈n+ is evaluated via eq. (148) using

After advancing the system by m substeps, we reconstruct the macroscopic state u at time τ + ∆τ (where ∆τ = ε3 ∆t) via the inverse relation u(ε(n − ct), ε3 t) = ε−2 rn (t)

(162)

for each nodal index n. In practice, we enforce periodic boundary conditions by applying a modulo operation to the spatial coordinates relative to the domain length, and subsequently interpolate the scattered u values onto a uniform spatial grid.

Data generation

Our numerical experiments encompass four distinct test cases: two continuum PDE models (the KdV and Allen–Cahn equations) and two microscopic particle systems (the FPUT and FENE chain models). We set the recording step size to ∆t = 0.001, resolving the internal dynamics with 25 intermediate integration substeps per recorded step. The total simulation horizon is T = 0.1, yielding exactly 100 temporally equidistant snapshots per trajectory. For 1D cases, we set the grid size as 256, and for 2D cases, we set the spatial resolution as 128 × 128. The data generation process for the PDE models are implemented via the Python package dedalus1 [8] with the SBDF2 numerical stepper. The generation for the microscopic chain models follows section E.3. For the microscopic systems, the FPUT chain parameters are selected to correspond to a canonical continuum KdV equation with akdv = 1 and bkdv = 1/24, which fixes the microscopic parameters at c = 1 and α = 1. For the FENE chain model, we set the parameters to H = 1 and R = 50ε2 . The spatial scaling parameter ε is fixed at 0.05 for the FPUT chain and 0.03 for the FENE chain. The initial conditions across all four datasets are generated using a consistent methodology based on random superpositions of sinusoidal waves. For the PDE datasets, these superpositions are used to directly initialize the continuum solution variable u. For the particle chain models, we first generate a continuum state u, which is subsequently projected onto the discrete microscopic strain variables r.

Networks and training details

G.1

Network architectures

We describe the detailed architecture for each model below. FNO We employ the FNO module officially implemented in the Python package neuralop2 [26]. We set the number of modes as 24, and the number of hidden channels as 32. All the other parameters remain unchanged. 1 https://dedalus-project.org/ 2 https://github.com/neuraloperator/neuraloperator

OnsagerNet We use the classical OnsagerNet formulation that treats the Nx -point PDE grid directly as an Nx -dimensional ODE, evolving under du = − M (u) + W (u) ∇V (u), dt

(163)

where V is a coercive MLP potential with 3 hidden layers of width 256 and tanh activations; M (u) = ϵI + B(u)B(u)⊤ is a low-r ank SPD dissipation matrix of rank 8 with ϵ = 10−3 , and B parameterized by a 2-hidden-layer MLP of width 256; and W (u) = B(u)C ⊤ − CB(u)⊤ is a low-rank anti-symmetric conservation matrix of rank 4, similarly parameterized. SpecOnsNet We use a pointwise MLP applied independently at each grid point to implement Fϕ . For the gradient term ∇vϕ , we represent vϕ as a complex-valued MLP acting on the first 16 Fourier modes with 3 hidden layers of width 64, and the graidient is computed via automatic differentiation. All the activations in Fϕ and vϕ are set as SiLU [22], except for the activations for the final nonlinear layer where we set them as sine to ensure the boundness of the network output.

G.2

Training details

We adopt the Adam optimizer with learning rate 10−3 . No weight decay is used. To manage the learning rate, we decrease the learning rate to half if the validation loss does not decrease for 300 epochs. We validate the networks every 100 epochs, and terminate the training process when the validation loss does not decrease for 5 validation steps. For the loss function, we fix the number of gradient accumulation steps as K = 5 for the FPUT model and K = 3 for the other cases.

Ablation study

We have conducted an ablation study for the modeling choice of our Spectral OnsagerNet. Apart from the baselines (classical solver and FNO), we introduce the following variants of our model. Residual OnsagerNet (abbreviated as Res-OnsagerNet) Instead of using eq. (47), we introduce s 2 s s µ̂Res (164) θ (u ) = α + (2πk) β ⊙ û + Hϕ (û ). where we directly parameterize all the gradient terms Hϕ as a complex-valued network. This is the most naive approach to learning the increment in the spectral space. Note that in such circumstances, the parameterization for the term µt is not guaranteed to be a functional derivative. Spectral OnsagerNet with real V (abbreviated as SpecOnsNet-V) In the update eq. (47) of SpecOnsNet, we set the potential V as the real potential for the PDE models. More specifically, we set µ̂θ (us ) as VKdV and VAC for the KdV equation and the Allen–Cahn equation, respectively.

Model

#params

KdV

Allen–Cahn

Chain–KdV

FENE

Classical solver SpecOnsNet-V

— 258

0.0541±0.0313 0.0096±0.0039

0.0085±0.0029 0.0059±0.0020

— —

FNO[28] OnsagerNet[40] Res-OnsagerNet SpecOnsNet

66.3 K 1.3 M 790 K 18.2 K

0.1747±0.0820 0.0671±0.0241 0.0129±0.0027 0.0100±0.0037

0.0204±0.0081 0.0260±0.0037 0.1225±0.0277 0.0055±0.0022

0.0360±0.0096 0.0242±0.0075 0.0156±0.0043 0.0082±0.0026

0.0186±0.0066 0.0139±0.0076 0.0100±0.0058 0.0012±0.0004

Table 3: 5-step relative prediction error (↓) for various baselines. Some of the entries have been exhibited in table 1 from the main text.

References [1] D. Bambusi and A. Ponno, On metastability in FPU, Commun. Math. Phys., 264 (2006), pp. 539–561, https://doi.org/10.1007/s00220-005-1488-1. [2] T. Bertalan, F. Dietrich, I. Mezić, and I. G. Kevrekidis, On learning Hamiltonian systems from data, Chaos, 29 (2019), p. 121107, https://doi.org/10.1063/1.5128231. [3] L. Bertini and G. Giacomin, Stochastic Burgers and KPZ equations from particle systems, Comm. Math. Phys., 183 (1997), pp. 571–607, https://doi.org/10.1007/s002200050044. [4] J. Brandstetter, D. Worrall, and M. Welling, Message passing neural PDE solvers, in International Conference on Learning Representations, 2022. [5] S. C. Brenner and L. R. Scott, The Mathematical Theory of Finite Element Methods, Springer New York, NY, 2008, https://doi.org/10.1007/978-0-387-75934-0. [6] S. L. Brunton and J. N. Kutz, Data-driven science and engineering: Machine learning, dynamical systems, and control, Cambridge University Press, 2019. [7] S. L. Brunton, J. L. Proctor, and J. N. Kutz, Discovering governing equations from data by sparse identification of nonlinear dynamical systems, Proc. Natl. Acad. Sci. USA, 113 (2016), pp. 3932–3937, https://doi.org/10.1073/pnas.1517384113. [8] K. J. Burns, G. M. Vasil, J. S. Oishi, D. Lecoanet, and B. P. Brown, Dedalus: A flexible framework for numerical simulations with spectral methods, Phys. Rev. Res., 2 (2020), p. 023068, https://doi.org/10.1103/PhysRevResearch.2.023068. [9] E. Celledoni, A. Leone, D. Murari, and B. Owren, Learning hamiltonians of constrained mechanical systems, Journal of Computational and Applied Mathematics, 417 (2023), p. 114608. [10] K. Champion, B. Lusch, J. N. Kutz, and S. L. Brunton, Data-driven discovery of coordinates and governing equations, Proc. Natl. Acad. Sci. USA, 116 (2019), pp. 22445–22451, https://doi.org/10.1073/pnas.1906995116. [11] J. Chen, K. Wu, and D. Xiu, DUE: A deep learning framework and library for modeling unknown equations, SIAM Rev., 67 (2025), pp. 873–902, https://doi.org/10.1137/24M1671827. [12] R. T. Q. Chen, Y. Rubanova, J. Bettencourt, and D. Duvenaud, Neural ordinary differential equations, in Proceedings of the 32nd International Conference on Neural Information Processing Systems, Curran Associates Inc., 2018, p. 6572–6583. [13] X. Chen, B. W. Soh, Z.-E. Ooi, E. Vissol-Gaudin, H. Yu, K. S. Novoselov, K. Hippalgaonkar, and Q. Li, Constructing custom thermodynamics using deep learning, Nat. Comput. Sci., 4 (2024), pp. 66–85, https://doi.org/10.1038/s43588-023-00581-5. [14] I. Corwin, The Kardar–Parisi–Zhang equation and universality class, Random Matrices Theory Appl., 1 (2012), p. 1130001, https://doi.org/10.1142/S2010326311300014. [15] Q. Du, Y. Gu, H. Yang, and C. Zhou, The discovery of dynamics via linear multistep methods and deep learning: error estimation, SIAM J. Numer. Anal., 60 (2022), pp. 2014–2045, https://doi.org/10.1137/21M140691X. [16] W. E, Principles of Multiscale Modeling, Cambridge University Press, 2011. [17] L. C. Evans, Partial Differential Equations, vol. 19, American Mathematical Society, 2010. [18] E. Fermi, P. Pasta, S. Ulam, and M. Tsingou, Studies of the nonlinear problems, tech. report, Los Alamos Scientific Lab., N. Mex., 1955.

[19] I. Gallagher, L. Saint-Raymond, and B. Texier, From Newton to Boltzmann: Hard Spheres and Short-range Potentials, European Mathematical Society Zürich, 2013, https:// doi.org/10.4171/129. [20] S. Greydanus, M. Dzamba, and J. Yosinski, Hamiltonian neural networks, in Proceedings of the 33rd International Conference on Neural Information Processing Systems, Curran Associates Inc., 2019, pp. 15379–15389. [21] M. D. Hansen, E. Celledoni, and B. K. Tapley, Learning mechanical systems from realworld data using discrete forced lagrangian dynamics, arXiv preprint arXiv:2505.20370, (2025). [22] D. Hendrycks and K. Gimpel, Gaussian error linear units (GELUs), arXiv preprint arXiv:1606.08415, (2016). [23] S. Huang, Z. He, N. Dirr, J. Zimmer, and C. Reina, Statistical-Physics-Informed Neural Networks (Stat-PINNs): A machine learning strategy for coarse-graining dissipative dynamics, J. Mech. Phys. Solids, 194 (2025), p. 105908, https://doi.org/10.1016/j.jmps.2024. 105908. [24] M. C. Kennedy and A. O’Hagan, Bayesian calibration of computer models, J. R. Stat. Soc. Ser. B. Stat. Methodol., 63 (2001), pp. 425–464, https://doi.org/10.1111/1467-9868. 00294. [25] I. G. Kevrekidis and G. Samaey, Equation-free multiscale computation: Algorithms and applications, Annu. Rev. Phys. Chem., 60 (2009), pp. 321–344, https://doi.org/10.1146/ annurev.physchem.59.032607.093610. [26] J. Kossaifi, N. Kovachki, Z. Li, D. Pitt, M. Liu-Schiaffini, R. J. George, B. Bonev, K. Azizzadenesheli, J. Berner, V. Duruisseaux, and A. Anandkumar, A library for learning neural operators, arXiv preprint arXiv:2412.10354, (2025). [27] K. Lee, N. Trask, and P. Stinis, Machine learning structure preserving brackets for forecasting irreversible processes, in Proceedings of the 35th International Conference on Neural Information Processing Systems, Curran Associates Inc., 2021, pp. 5696–5707. [28] Z. Li, N. B. Kovachki, K. Azizzadenesheli, B. liu, K. Bhattacharya, A. Stuart, and A. Anandkumar, Fourier neural operator for parametric partial differential equations, in International Conference on Learning Representations, 2021. [29] B. List, L.-W. Chen, K. Bali, and N. Thuerey, Differentiability in unrolled training of neural physics simulators on transient dynamics, Comput. Methods Appl. Mech. Engrg., 433 (2025), p. 117441, https://doi.org/10.1016/j.cma.2024.117441. [30] L. Lu, P. Jin, G. Pang, Z. Zhang, and G. E. Karniadakis, Learning nonlinear operators via DeepONet based on the universal approximation theorem of operators, Nat. Mach. Intell., 3 (2021), pp. 218–229, https://doi.org/10.1038/s42256-021-00302-5. [31] D. Masao, Onsager’s variational principle in soft matter, J. Phys. Condens. Matter, 23 (2011), p. 284118, https://doi.org/10.1088/0953-8984/23/28/284118. [32] D. Masao, Onsager principle as a tool for approximation, Chinese Phys. B, 24 (2015), p. 020505, https://doi.org/10.1088/1674-1056/24/2/020505. [33] W. L. Oberkampf and C. J. Roy, Verification and Validation in Scientific Computing, Cambridge university press, 2010, https://doi.org/10.1017/CBO9780511760396. [34] L. Onsager, Reciprocal relations in irreversible processes. I., Phys. Rev., 37 (1931), pp. 405– 426, https://doi.org/10.1103/PhysRev.37.405.

[35] L. Onsager, Reciprocal relations in irreversible processes. II., Phys. Rev., 38 (1931), pp. 2265– 2279, https://doi.org/10.1103/PhysRev.38.2265. [36] G. Pavliotis and A. Stuart, Multiscale Methods: Averaging and Homogenization, Springer New York, NY, 2008, https://doi.org/10.1007/978-0-387-73829-1. [37] N. Provatas and K. Elder, Phase-Field Methods in Materials Science and Engineering, John Wiley & Sons, Ltd, 2010, https://doi.org/10.1002/9783527631520. [38] H. R. Warner Jr, Kinetic theory and rheology of dilute suspensions of finitely extendible dumbbells, Industrial & Engineering Chemistry Fundamentals, 11 (1972), pp. 379–387, https: //doi.org/10.1021/i160043a017. [39] K. Wu, T. Qin, and D. Xiu, Structure-preserving method for reconstructing unknown hamiltonian systems from trajectory data, SIAM J. Sci. Comput., 42 (2020), pp. A3704–A3729, https://doi.org/10.1137/19M1264011. [40] H. Yu, X. Tian, W. E, and Q. Li, OnsagerNet: Learning stable and interpretable dynamics using a generalized Onsager principle, Phys. Rev. Fluids, 6 (2021), p. 114402, https://doi. org/10.1103/PhysRevFluids.6.114402. [41] N. J. Zabusky and M. D. Kruskal, Interaction of “solitons” in a collisionless plasma and the recurrence of initial states, Phys. Rev. Lett., 15 (1965), p. 240, https://doi.org/10.1103/ PhysRevLett.15.240. [42] Z. Zhang, Y. Shin, and G. Em Karniadakis, GFINNs: GENERIC formalism informed neural networks for deterministic and stochastic dynamical systems, Philos. Trans. Roy. Soc. A, 380 (2022), p. 20210207, https://doi.org/10.1098/rsta.2021.0207. [43] A. Zhu, B. W. Soh, G. A. Pavliotis, and Q. Li, Identifiable learning of dissipative dynamics, arXiv preprint arXiv:2510.24160, (2025).

Record · ID 196442 · SHA-256 a7d20b42307abc8d
Retrieved via Conceptio — every document is proof-bundled with source, license, and retrieval metadata.