From Langevin dynamics to macroscopic thermodynamic models: a general framework valid far from equilibrium

Travis Leadbetter; Prashant K. Purohit; Celia Reina

doi:10.1515/jnet-2025-0071

Article Open Access

From Langevin dynamics to macroscopic thermodynamic models: a general framework valid far from equilibrium

Travis Leadbetter , Prashant K. Purohit and Celia Reina

Published/Copyright: October 16, 2025

Published by

Become an author with De Gruyter Brill

Submit Manuscript Author Information

From the journal Journal of Non-Equilibrium Thermodynamics

Abstract

Given a particle system obeying overdamped Langevin dynamics, we demonstrate that it is always possible to construct a thermodynamically consistent macroscopic model which obeys a gradient flow with respect to its non-equilibrium free energy. To do so, we significantly extend the recent Stochastic Thermodynamics with Internal Variables (STIV) framework, a method for producing macroscopic thermodynamic models far-from-equilibrium from the underlying mesoscopic dynamics and an approximate probability density of states parameterized with so-called internal variables. Though originally explored for Gaussian probability distributions, we here allow for an arbitrary choice of the approximate probability density while retaining a gradient flow dynamics. This greatly extends its range of applicability and automatically ensures consistency with the second law of thermodynamics, without the need for secondary verification. We demonstrate numerical convergence, in the limit of increasing internal variables, to the true probability density of states for both a multi-modal relaxation problem, a protein diffusing on a strand of DNA, and for an externally driven particle in a periodic landscape. Finally, we provide a reformulation of STIV with the quasi-equilibrium approximations in terms of the averages of observables of the mesostate, and show that these, too, obey a gradient flow.

Keywords: stochastic thermodynamics; thermodynamics with internal variables; gradient flow

1 Introduction

The field of non-equilibrium thermodynamics boldly seeks to describe the macroscopic behavior of a wide range of physical systems. Using various formulations, such as irreversible thermodynamics [1], [2], rational thermodynamics [3], [4], thermodynamics with internal variables [5], or the General Equation for Non-Equilibrium Reversible-Irreversible Coupling (GENERIC) formalism [6], [7], [8], complex phenomena like creep [9], damage [10], [11], [12], thermoplatiscity [7], and fluid rheology [13], to name a few, have been modeled effectively and systematically. However, central to our understanding, and at the heart of Hilbert’s sixth problem, is the ability to formally and practically connect these macroscopic models to more fundamental microscopic physics. Perhaps the two most important contributions on this front are the Mori–Zwanzig formalism and other projection operator approaches which seek to project out noisy or irrelevant aspects of a Hamiltonian system, leaving only quantities which vary over long time scales, and the utilization of large deviation principles to coarse-grain microscopic stochastic processes [14], [15], [16], [17]. Although these two theories provide the cornerstone of our theoretical understanding, significant challenges limit their utility in producing usable models.

In a previous work [18], the authors leveraged the theory of stochastic thermodynamics [19], [20] and a variational principle due to Eyink [21] to derive thermodynamically consistent macroscopic models directly from microscopic particle systems obeying an overdamped Langevin dynamics of the form

(1) d x t = − 1 η ∂ x e d t + 2 η β d w t ,

where the interaction energy e(x, λ) is potentially dependent on an external driving protocol λ, η and β are the viscosity and inverse temperature respectively (both assumed to be constant), and w _t is a standard Brownian motion. Stochastic thermodynamics provided the definitions for thermodynamic quantities both at the level of individual trajectories and at the macroscopic scale via ensemble averages. The system was then modeled using a parametrized approximate probability density of states, where the dynamics of the parameters (called internal variables) were determined using the variational method of Eyink. One obtained an approximate macroscopic thermodynamic model with internal variables by replacing the true probability density of states with its parameterized approximation in the definitions of macroscopic quantities from stochastic thermodynamics. The resulting macroscopic models exhibited the thermodynamics with internal variable structure for an arbitrary choice of approximation to the probability density of states, and hence the framework was named Stochastic Thermodynamics with Internal Variables (STIV). For the particular choice of a Gaussian approximation to the probability density of states, the dynamics of the internal variables obeyed a gradient flow with respect to the non-equilibrium free energy, A ̂ neq , thereby guaranteeing a non-negative rate of total entropy production. The STIV framework with a Gaussian approximation to the probability density of states has proven useful so far in providing the first ab-initio derivation of phase field models of front propagation in elastic media [22] as well as capturing the evolution and spreading of orientations of cell populations [23].

In this work, we expand the utility of the STIV framework by showing that the Gaussian approximation is not the only way to ensure the gradient flow structure, and hence a non-negative rate of total entropy production. In fact, the gradient flow structure can be achieved for any choice of approximate probability density of states. This is fitting, as it is well known that the Fokker–Planck equation associated with Langevin dynamics can be formulated as a gradient flow of the free energy A neq = ∫ e + β − 1 ⁡ log ⁡ p p d x with respect to the Wasserstein-2 metric [24]

(2) ∂ t p = 1 η ∂ x i p ∂ x i δ A neq δ p .

Preserving the gradient flow structure offers several benefits. As already mentioned, within the STIV framework it guarantees that the approximate rate of total entropy production is non-negative and ensures thermodynamic consistency with the second law of thermodynamics. Typically, the second law must be enforced in an additional step, such as through the Coleman–Noll procedure [25]. Structure preservation has also been found both to improve the numerical accuracy of an approximation, and lead to more accurate qualitative behaviors [26]. This has proved vital to the numerical study of complex systems in many fields, including celestial mechanics [27], [28], molecular dynamics [29], control theory [30], and non-equilibrium mechanics [31].

2 Derivation

The dynamical equations for the STIV framework are derived from the variational method of Eyink [21], which provides a method for approximating the solution to differential equations arising in non-equilibrium statistical mechanics. Starting from the Fokker–Planck equation associated with the microscopic stochastic differential equation in Eq. (1) ^[1]

∂ t p ( x , t ) = L p ( x , t ) = ∂ x i 1 η p ( x , t ) ∂ x i e + 1 η β ∂ x i p ( x , t )

one can recast the problem in a variational form by introducing a test function ψ, and the action

Γ [ ψ , p ] = ∫ 0 ∞ ∫ X ψ ∂ t − L p d x d t .

Stationary points (ψ*, p*) of this action are solutions to the adjoint and standard forms of the original Fokker–Planck equation respectively, and solve an infinite dimensional Hamilton equation

∂ t p = δ δ ψ H [ ψ , p ] ∂ t ψ = − δ δ p H [ ψ , p ]

where H = ∫ X ψ L p d x is the Hamiltonian. The adjoint equation is trivially solvable with solution ψ* ≡ 1. The equation for p is intractable in general, so to find an approximate solution one proposes a parameterized approximate probability density of states, p ̂ ( x , α ) , and a parameterized test function, ψ ̂ ( x , α ) . Writing u ̂ = log ( p ̂ ) and f p ̂ = ∫ f ( x ) p ̂ ( x , α ) d x , the dynamical equations for the parameters α reduce to

[ ∂ ψ ̂ , ∂ u ̂ ] i j p ̂ α ̇ j = ∂ α i L † ψ ̂ p ̂

where [ ∂ ψ ̂ , ∂ u ̂ ] i j = ∂ α i ψ ̂ ∂ α j u ̂ − ∂ α j ψ ̂ ∂ α i u ̂ and L † is the adjoint of L . The matrix [ ∂ ψ ̂ , ∂ u ̂ ] i j p ̂ is anti-symmetric and obeys the Jacobi identity, so these dynamics retain the Hamiltonian form.

In the original framework [18], the parameters were split into two sets: the variables parameterizing p denoted (with a slight abuse of notation) α and referred to as internal variables, and dummy variables γ of the same dimension as α. The approximate density is assumed to only depend on α while ψ ̂ was defined as ψ ̂ = 1 + γ i ∂ α i s ̂ , where s ̂ = − k B u ̂ = − k B ⁡ log ( p ̂ ) is the approximate microscopic entropy, and k _B is the Boltzmann constant. With this particular choice, the dynamical equations reduce to

∂ α i s ̂ ∂ α j s ̂ p ̂ α ̇ j = − k B L † ∂ α i s ̂ p ̂ , γ ( t ) ≡ 0 .

This “linear ansatz” using the α-gradient of the microscopic entropy in the test function ψ ̂ produces a specific case of the moment-closure equations called “entropic matching” [32], and ultimately leads to irreversible dynamics for α.

Once the dynamics for the internal variables are obtained, one can approximate macroscopic thermodynamic quantities defined via stochastic thermodynamics by substituting the approximate probability density of states for the true solution. For example, the macroscopic non-equilibrium free energy defined as A neq = e p − T s p becomes A ̂ neq = e p ̂ − T s ̂ p ̂ , from which one can approximate both the external force F ex = ∂ λ e p by F ̂ ex = ∂ λ A ̂ neq and the rate of total entropy production d d t S tot by d d t S ̂ tot = 1 T F ̂ ex λ ̇ − d d t A ̂ neq = − 1 T ∂ α i A ̂ neq α ̇ i .

Previously, the STIV framework was effectively used to study the dynamics and thermodynamics of phase front propagation in an elastic medium, where one could derive both the continuum limit [18], and an equivalent phase field model directly from microscopic physics [22]. This previous implementation, however, was limited to a multivariate Gaussian approximation to the probability density of states p ̂ as this was the only form in which the gradient flow structure, and hence non-decreasing total entropy production, could be ensured at the time. The key findings of this work are a number of options for preserving the gradient flow structure of the underlying Fokker–Planck equation which are flexible enough to approximate any probability density to arbitrary accuracy in the limit of increasing internal variables. The first method, which we call the “operator method”, utilizes an alternative form for the approximate test function that allows for a completely unconstrained choice of approximate probability density of states. The second uses the original approximate test functions and entropic matching moment-closure equations, but requires one to restrict the approximate probability density of states using a “quasi-equilibrium” distributions [33]. The final method forgoes the need to include the Boltzmann factor exp(−βe(x, λ)) in the approximate probability density of states, which could reduce computational costs in some cases, but the internal variable dynamics given by entropic matching only approximate a gradient flow. We label this method as “exponential family” throughout. A detailed analysis of each of these three methods is provided next, and the key elements and results for each of them are summarized in Table 1 as reference.

Table 1:

Key features of the closure methods.^a

Method name	Operator method	Quasi-equilibrium	Exponential family
Probability density of states p ̂	Arbitrary	exp(−βe + αϕ − f)	exp(αϕ − f)
Test functions ψ ̂	1 + γ L u ̂ † − 1 ∂ α u ̂	1 + γ ∂ α u ̂	1 + γ ∂ α u ̂
Dynamical equations	F α ̇ = − 1 η ∂ α A ̂ neq	C α ̇ = − 1 η M C − 1 ∂ α A ̂ neq	C α ̇ = − 1 η M C − 1 ∂ α A ̂ neq ^c
Dynamics with external driving	F α ̇ = − 1 η ∂ α A ̂ neq	C α ̇ + ϕ ̄ ∂ λ u ̂ p ̂ λ ̇ = − 1 η M C − 1 ∂ α A ̂ neq	C α ̇ = − 1 η M C − 1 ∂ α A ̂ neq
Dual dynamics	Not studied here	d d t Φ ̂ = − 1 η M ∂ Φ A ̂ neq ^b	d d t Φ ̂ = − 1 η M ∂ Φ A ̂ neq
Entropy production T d d t S ̂ tot	η F i j α ̇ i α ̇ j	1 η β 2 M i j α i α j	η M i j − 1 d d t Φ ̂ i d d t Φ ̂ j

^aDefinition of terms for Table 1: C i j = ∂ α i u ̂ ∂ α j u ̂ p ̂ = Cov p ̂ ( ϕ i , ϕ j ) , M i j = ∂ x k ∂ α i u ̂ ∂ x k ∂ α j u ̂ p ̂ = ∂ x k ϕ i ∂ x k ϕ j p ̂ , L u ̂ † g = ∂ x i u ̂ ∂ x i g + ∂ x i ∂ x i g , F i j = ∂ α i u ̂ L u ̂ † − 1 ∂ α j u ̂ p ̂ , ϕ ̄ = ϕ − Φ ̂ , Φ ̂ = ϕ p ̂ . ^bThis equation holds with and without external driving. ^cThese dynamics are not given by moment closure, but approximate it for large number of observables.

2.1 The operator method

We now motivate an alternative choice for the approximate test function ψ ̂ to show how the STIV method can produce gradient flow equations for an arbitrary approximate probability density of states. Again, we fix an approximate density p ̂ ( x , α ) , and let u ̂ = log ( p ̂ ) . We then assume that χ(x, α) is a vector valued function with the same dimension as α, and define ψ ̂ ( x , α , γ ) = 1 + γ i χ i ( x , α ) . Following the same steps as above, the dynamical equations induced by the variational method become

χ i ∂ α j u ̂ p ̂ α ̇ j = L † χ i p ̂ , γ ( t ) = 0 .

For convenience, we define F i j ≡ χ i ∂ α j u ̂ p ̂ and C i j ≡ ∂ α i u ̂ ∂ α j u ̂ p ̂ .

We have not specified new test functions yet, but to do so we examine the approximate rate of total entropy production, T S ̂ ̇ tot = − α ̇ i ∂ α i A ̂ neq . We recall that A ̂ neq = E ̂ − T S ̂ = e p ̂ + u ̂ / β p ̂ , and since 0 = ∂ α i 1 p ̂ = ∂ α i u ̂ p ̂ and e does not depend on α, we have ∂ α i A ̂ neq = ( e + u ̂ / β ) ∂ α i u ̂ p ̂ . On the other hand, repeated integration by parts reveals that

L † χ i p ̂ = − 1 η ∂ x k e ∂ x k χ i + d ∂ x l ∂ x l χ i p ̂ = − 1 η ∂ x k e + u ̂ / β ∂ x k χ i p ̂ = 1 η ( e + u ̂ / β ) ∂ x k u ̂ ∂ x k χ i + ∂ x k ∂ x k χ i p ̂ ≡ 1 η ( e + u ̂ / β ) L u ̂ † χ i p ̂ .

This new operator L u ̂ † g = ∂ x k u ̂ ∂ x k g + ∂ x l ∂ x l g is also the adjoint of a Fokker Planck operator which only depends on u ̂ . The final line resembles the equation for ∂ α i A ̂ neq , and suggests that one should choose test functions which obey

(3) L u ̂ † χ i = − ∂ α i u ̂ .

Assuming one can solve this differential equation, the dynamical equations become a gradient flow

(4) F i j α ̇ j = − 1 η ∂ α i A ̂ neq ,

since F i j is positive semi-definite as

F i j = χ i ∂ α j u ̂ p ̂ = − χ i L u ̂ † χ j p ̂ = − ∫ χ i ∂ x k p ̂ ∂ x k χ j d x = ∂ x k χ i ∂ x k χ j p ̂ .

Between the second and third lines, we have made use of the fact that p ̂ L u ̂ † χ j = ∂ x k p ̂ ∂ x k χ j . With the gradient flow structure comes the guarantee of non-negative rate of total entropy production

T S ̂ ̇ tot = η F i j α ̇ i α ̇ j ≥ 0 .

This strategy allows for a completely general choice of approximate probability density of states p ̂ ( x , α ) . Solving Eq. (3), though, is challenging in general. However, following a well known trick [34], one can define a “Hamiltonian” conjugation of L u ̂ † by

− H u ̂ g = p ̂ L u ̂ † g p ̂ = ∂ x 2 g − v ( x , α ) g

where v = ∂ x 2 u ̂ / 2 + | ∂ x u ̂ | 2 / 4 . Unlike L u ̂ † , H u ̂ is self-adjoint with respect to the Lebesgue measure. Note, also that H u ̂ defines a Schrödinger operator, and it is immediately clear that if (λ*, g*) is an eigenvalue-eigenvector pair for H u ̂ , then ( − λ * , g * / p ̂ ) is an eigenvalue-eigenvector pair for L u ̂ † . Hence, one can take advantage of existing eigenvalue solvers for quantum systems in order to compute the necessary test functions.

2.2 Quasi-equilibrium distribution

The primary advantage of the previous method for producing a gradient flow dynamics is that it allows for an arbitrary choice of approximate probability density of states. The disadvantage is that the test functions must be computed as solutions to the differential equation in Eq. (3). Alternatively, we show now that if one restricts the approximate probability density of states to a specific class known as quasi-equilibrium densities (a specific form of exponential family) [33], then one may use the original test functions χ i = ∂ α i u ̂ and still achieve a gradient flow (here, we have replaced s ̂ by u ̂ for notational simplicity). This is only a minor restriction as any continuous probability distribution can be approximated by an exponential family or quasi-equilibrium distribution to arbitrary accuracy (see Appendix C for one such proof).

Assume that the approximate probability density of states takes a quasi-equilibrium form. That is

(5) p ̂ ( x , α , β ) = exp ( − β e ( x ) + α i ϕ i ( x ) − f ( α , β ) ) ,

for some set of “observables” ϕ _i(x) and internal variables α _i(t). Here, β (taken as a constant) and the internal variables α, form the natural parameters of the exponential family, while the energy e and observables ϕ form the “sufficient statistics” [35]. It will be helpful to define Φ ̂ i as the p ̂ averages of the observables, i.e., Φ ̂ i ≡ ϕ i p ̂ .

To see the gradient flow structure, note first that in this case

χ i ∂ α j u ̂ p ̂ = ( ϕ i − Φ ̂ i ) ( ϕ j − Φ ̂ j ) p ̂ = C i j .

Second, since e + u ̂ / β = ( α i ϕ i − f ( α , β ) ) / β , we have

∂ α i A ̂ neq = ( e + u ̂ / β ) ∂ α i u ̂ p ̂ = C i j α j β .

Finally, again by integrating by parts

L † χ i p ̂ = − 1 η ∂ x k ( e + u ̂ / β ) ∂ x k χ i p ̂ = − ∂ x k ϕ i ∂ x k ϕ j p ̂ α j η β ≡ − 1 η β M i j α j .

As with F in the earlier discussion, M i j = ∂ x k ϕ i ∂ x k ϕ j p ̂ is positive semi-definite. Putting the pieces together, we get

(6) C i j α ̇ j = − 1 η β M i j α j

or, by inserting the identity as δ j k = C j l − 1 C l k between M and α (we assume C is invertible),

(7) α ̇ i = − 1 η C i j − 1 M j k C k l − 1 ∂ α l A ̂ neq .

Since M is positive semi-definite, and C is symmetric, C − 1 M C − 1 is positive semi-definite and the dynamics obey a gradient flow.

2.3 Exponential family

Using a quasi-equilibrium distribution which includes the total energy leads to gradient flow dynamics exactly. However, one may not always wish to include the energy in the approximate density of states, for example, if it would make computing the normalizing function more challenging. If one removes the energy and uses an exponential family of the form

(8) p ̂ ( x , α ) = exp ( α i ϕ i − f ( α ) )

the dynamics are not always a gradient flow. Specifically, they become

(9) C α ̇ = − 1 η M C − 1 ∂ α A ̂ neq − 1 η ∂ x k e ∂ x k ϕ p ̂ − M C − 1 e ( ϕ − Φ ̂ ) p ̂ .

The second term on the right hand side only depends on the (non-constant) part of the energy which is p ̂ -orthogonal to the span of the ϕ. Mathematically, if we can write e ( x ) = γ j ( ϕ j ( x ) − Φ ̂ j ) + e 0 + e ⊥ ( x ) where e ⊥ ( ϕ j − Φ ̂ j ) p ̂ = 0 for all j, then the second term becomes

∂ x k e ∂ x k ϕ p ̂ − M C − 1 e ( ϕ − Φ ̂ ) p ̂ = ∂ x k e ⊥ ∂ x k ϕ p ̂ .

If the observables are chosen to produce a highly flexible approximation, for example by using machine learning techniques or a basis expansion, it is reasonable to assume that this term is negligible and therefore one can set it to zero, giving a gradient flow structure. In this case, ∂ α A ̂ neq = C α / β + e ( ϕ − Φ ̂ ) p ̂ , and again the dynamical equations are given by Eq. (7). To minimize confusion, we shall refer to an approximate probability density of states of the form p ̂ = exp ( α i ϕ i − f ( α ) ) which does not include the energy as an exponential family approximation. When p ̂ = exp ( − β e + α i ϕ i − f ( α , β ) ) explicitly includes the energy, we shall refer to this as a quasi-equilibrium approximation (although it is also an exponential family).

2.4 Time dependent external driving

Now we consider the case when the total energy is time dependent through some external driving protocol, e = e(x, λ(t)). In the operator formulation, with test functions L u χ i = ∂ α i u ̂ , nothing needs change to maintain the original STIV definitions. Since the approximate probability density of states can still be chosen to be independent of λ (for fixed α), the approximate external force remains F ̂ ex = ∂ λ A ̂ neq [18]. In the quasi-equilibrium formulation, p ̂ explicitly depends on λ through e and so it is no longer true that F ̂ ex = ∂ λ A ̂ neq . Moreover, when deriving the dynamical equations from the variational method [21], it was assumed that all of the time dependence of p ̂ was captured through the internal variables α. When p ̂ takes the quasi-equilibrium form this is not the case since the energy is now time dependent. This is fixed through the proper modification of the dynamical equations given in Eq. (6) (which can be found in Appendix A). The resulting dynamical equations are

(10) C i j α ̇ j + ( ϕ i − Φ ̂ i ) ∂ λ u ̂ p ̂ λ ̇ = − 1 η β M i k α k

and the external force becomes

(11) F ̂ ex = ∂ λ A ̂ neq − α i Φ ̂ i β = − 1 β ∂ λ f ( α , β , λ ) .

The equation for the rate of total entropy production is also impacted, but it still remains non-negative

(12) T d d t S ̂ tot = 1 η β 2 α i M i j α j

(see the Appendix B for the derivation of both the approximate external force and approximate rate of total entropy production).

The new dynamical equations appear to break the gradient flow structure. However, in the following section, we show that the gradient flow structure remains intact when the equations are recast in terms of the dynamics of the averages of observables Φ ̂ .

2.5 Dual formulation

So far, the dynamics of the system have been tracked through the internal variables α. However, in traditional internal variable formulations, the internal variables are generally assumed to be averages of microscopic substructures which are not captured in equilibrium. Thus, it seems more natural to consider the averages of the observables, that is Φ ̂ i = ϕ i p ̂ as independent variables, and develop the system’s dynamics in terms of their behavior. Fortunately, this is readily achievable through a simple change of variables. Assuming the map from the internal variables to the averages of the observables, α → Φ ̂ , is invertible, we can define the non-equilibrium free energy as a function of Φ ̂ . Using the chain rule, we find that ∂ α A ̂ neq = C ∂ Φ A ̂ neq . Combining this with

d d t Φ ̂ i = ϕ i ∂ α j u p ̂ α ̇ j + ϕ i ∂ λ u ̂ p ̂ λ ̇ = C i j α ̇ j + ( ϕ i − Φ ̂ i ) ∂ λ u ̂ p ̂ λ ̇ ,

which reproduces the left hand side of Eq. (10) (or Eqs. (6) and (9) when p ̂ doesn’t depend on λ), gives the equivalent dynamics for Φ ̂

d d t Φ ̂ i = − 1 η M i j ∂ Φ j A ̂ neq

which again have the gradient flow structure. This holds true for both the quasi-equilibrium and exponential family approximations. Moreover, when p ̂ is independent of λ (either the exponential family approximation or the quasi-equilibrium approximation with no external driving), the form of the rate of total entropy production is invariant under the change of variables [22]

T d d t S ̂ tot = − ∂ α i A ̂ neq ( α i , β ) α ̇ i = − ∂ Φ i A ̂ neq ( Φ ̂ , β ) d Φ ̂ i d t = η d Φ ̂ i d t M i j − 1 d Φ ̂ j d t .

In this form, it becomes clear that M i j = ∂ x ϕ i ∂ x ϕ j p ̂ is an approximation to the true average dissipation matrix of the observables ϕ _i(x _t) along fluctuating trajectories. Mathematically, since the governing Langevin equation is given by

d x t = − ∂ x e / η d t + 2 / η β d w t ,

we know by Itô’s formula that

d ϕ i ( x t ) = − 1 η ∂ x e ∂ x ϕ i + 1 η β Δ x ϕ i d t + 2 η β ∂ x ϕ i d w t

and hence the diffusion matrix of ϕ(x _t) is M i j ( x t ) ∝ ∂ x ϕ i ∂ x ϕ j . When p ̂ depends on λ, the approximate rate of total entropy production cannot be written as a quadratic form in d d t Φ ̂ but is still given by Eq. (12).

It is also interesting to note that in the case of the quasi-equilibrium approximation, the non-equilibrium free energy is the convex dual of the normalizing function

f ( α , β , λ ) = log ∫ exp ( − β e + α ϕ ) d x .

f is convex in α since ∂ α i ∂ α j f = ( ϕ i − Φ ̂ i ) ( ϕ j − Φ ̂ j ) p ̂ = C i j is positive semi-definite and we know

Φ ̂ i = ϕ i p ̂ = ∂ α i f .

Using the Legendre transform, we can define the convex dual to f as f * ( Φ , β , λ ) = sup α α Φ − f ( α , β , λ ) so that

α i = ∂ Φ i f * ( Φ = Φ ̂ , β , λ ) .

Looking at the definition for A ̂ neq ( α , β , λ ) = e + u ̂ / β p ̂ = ( α i ϕ i − f ( α , β , λ ) p ̂ / β = ( α i Φ ̂ i − f ( α , β , λ ) ) / β , it’s clear that A ̂ neq ( Φ , β , λ ) = f * ( Φ , β , λ ) / β . Since A ̂ neq ( α , β , λ ) ≠ f ( α , β , λ ) , the non-equilibrium free energies driving the dynamics for α and Φ ̂ , that is A ̂ neq ( α ) and A ̂ neq ( Φ ) , are not related through a Legendre transform as is typically the case when one changes between dual variables in classical thermodynamics.

3 An information theoretic view

Before moving on to the example implementations of the STIV framework, we highlight an information theoretic perspective on the dynamical equations for the internal variables based on minimizing the Kullback–Leibler (KL) divergence between the approximate and true probability density of states. Following Bronstein and Koeppl [32], we assume that the true probability density of states if governed by the equation ∂ t p = L p , and try to minimize the KL-divergence

D K L ( α ) ≡ D K L [ p ̂ ‖ p ] = ∫ p ̂ log p ̂ p d x .

Differentiating in α gives

0 = ∫ log p ̂ p ∂ α i p ̂ d x .

We would like this to be true for all time, and hence we differentiate in time to get

0 = ∫ log p ̂ p ∂ α i ∂ α j p ̂ + p ̂ ∂ α i u ̂ ∂ α j u ̂ α ̇ j d x + ∫ log p ̂ p ∂ α i ∂ λ p ̂ + p ̂ ∂ α i u ̂ ∂ λ u ̂ λ ̇ − ∫ ∂ α i u ̂ L p p p ̂ d x

or in terms of expectations with respect to p ̂ ,

∂ α i u ̂ L p p p ̂ = log p ̂ p ∂ α i ∂ α j p ̂ p ̂ p ̂ α ̇ j + C i j α ̇ j + log p ̂ p ∂ α i ∂ λ p ̂ p ̂ p ̂ λ ̇ + ∂ α i u ̂ ∂ λ u ̂ p ̂ λ ̇ .

Since this equation contains p, it is unlikely that one could simplify it without significant approximation. Moreover, the best approximation one has for p at any given time within our framework is p ̂ . If we make this substitution in the previous equation, we achieve

C i j α ̇ j + ∂ α i u ̂ ∂ λ u ̂ p ̂ λ ̇ = L † ∂ α i u ̂ p ̂

which is exactly the original dynamical equations given by the variational method of Eyink. Thus, these dynamics minimize the KL-divergence between the true probability density of states and the approximate probability density of states at each instance in time.

4 Examples

As a demonstration of the enhanced flexibility of the modified STIV framework presented above, we compare the results of STIV using an exponential family and quasi-equilibrium approximation to the true probability density of states for two example systems. The first example is that of a colloidal protein diffusing along a string of DNA [36], [37]. This system is approximately described by a single particle diffusing in a one dimensional potential consisting of the product of a sinusoidal function and a decaying exponential. Mathematically,

(13) u ( x ) = c 0 ⁡ exp ( − c 1 x ) sin ( c 2 x + c 3 ) + bounding box .

The parameters c ₀, …, c ₃ were chosen to fit the potential shown in Figure 6 of [37], and the bounding box is specified to hold the particle within 60 base pairs of deepest minima (the values of the parameters and the bounding box are detailed in Appendix D). This example allows us to study the approximation error during a relaxation towards equilibrium. The second example includes time dependent driving. A colloidal particle is assumed to diffuse in a sinusoidal landscape while coupled to an external driving protocol via a harmonic potential, a prototypical model of an optical trap. The external driving moves in the positive x axis with constant velocity v. Mathematically,

(14) e ( x , λ ) = − cos 2 π x l + k 2 ( x − λ ) 2 ,

with λ(t) = vt (see Appendix D for parameter details). Here, we can evaluate the accumulation of error in time since this system never reaches a steady state.

For these examples, we compare the approximation accuracy of the quasi-equilibrium approximate density as in Eq. (5), and a traditional exponential family without the energy as in Eq. (8). Since both examples are one dimensional, we chose sinusoidal Fourier expansion functions as the observables ϕ _n = sin(nπ(x − x ₀)/L) where x ₀ and L are chosen to fit the bounding box in the diffusing protein example or the spatial extent of the simulation in the external driving example. Due to the computational complexity of solving the differential equation L u ̂ † χ i = ∂ α i u ̂ to find χ _i at each step, we choose not to implement the operator method. It is likely that this method will only become feasible for specific choices of observables for which the differential equation can be inverted exactly (as happens to be the case for a univariate Gaussian).

Figure 1 shows results for the colloidal protein example, and Figure 2 shows those of the time dependent driving example. Three different implementations of the STIV framework are highlighted throughout. The first, labeled “ p ̂ ∝ exp ( − β e + α ϕ ) ” and colored green, uses the quasi-equilibrium probability density of states and the dynamics obtained in Eq. (6). The second, labeled “ p ̂ ∝ exp ( α ϕ ) Moment Closure” and colored yellow, uses just the exponential family approximation and the dynamics obtained from the variational method of Eyink, Eq. (9). Finally, the method labeled “ p ̂ ∝ exp ( α ϕ ) Gradient Flow” and colored pink also uses the exponential family approximation but the second term on the right hand side of Eq. (9) (which is expected to vanish as the number of observables increases) is set to zero to obtain a gradient flow. The convergence of each method to the ground truth probability density of states is shown in panel A. The metric used is the time averaged total variation distance between the ground truth and the approximation

d ( p , p ̂ ) = 1 2 T ∫ 0 T ∫ R | p ( x , t ) − p ̂ ( x , α ( t ) ) | d x d t .

Figure 1:

The convergence of STIV to the true probability density of states in a model of protein diffusion on DNA. Panel A shows the time averaged total variation distance between each STIV model (the quasi-equilibrium approximation in green, the exponential family approximation without the energy in yellow, and the exponential family approximation with modified dynamics to produce a gradient flow in pink) and the true probability density of states (obtained via a traditional solver) as a function of the number of internal variables. Panel B depicts the energy landscape. Panels C and D show the approximate probability density of states for each model for 8 and 32 internal variables respectively and the true probability density of states (in dark blue) at the same fixed time point (midway between the initial time and fully relaxed in equilibrium).

Figure 2:

Investigation of the error accumulation of the STIV models for a single particle in a sinusoidal potential driven by a linear external force. Panel A shows the total variation distance between each model (same colors as Figure 1) and the true probability density of states when 64 internal variables are used. Panel B depicts the sinusoidal landscape and the quadratic potential associated with the external driving from an optical trap (at t = 0.79 t _final). Panels C and D show the approximate probability density of states for each model using 64 internal variables and the true probability density of states (in dark blue) for two different points in time. (t _final = 10).

Although previously we have made use of the KL divergence to interpret the choice of dynamical equations for the internal variables, here we use the total variation distance to measure the error in the approximation as the numerical computation of the KL divergence tended to be poorly behaved. Fortunately, the total variation distance is bounded above by the KL divergence by [38], [39], [40]

d ( p , q ) ≤ min 1 2 D KL ( p ‖ q ) , 1 − exp ( − D KL ( p ‖ q ) ) .

In panel B of Figure 1, we show the potential energy governing the system, along with a cartoon of the model setup. In the remaining panels C and D, snapshots of the densities at an intermediate time for two different numbers of observables are shown. In these frames, the ground truth probability density of states is shown in dark blue. At this time point, the probability density of states has not fully relaxed from the initial Gaussian distribution, and is still far from equilibrium. All methods converge to the ground truth solution, however, the quasi-equilibrium distribution outperforms the two exponential family models by an order of magnitude for 32 and 64 observables. As is expected, the interaction energy is a good observable and provides important information about the probability density of states, particularly in a problem concerning relaxation to the equilibrium state. Interestingly, the exponential family with gradient flow dynamics outperforms the moment-closure dynamics for small numbers of internal variables. As expected, this gap decreases for 64 internal variables once the approximate density has become highly flexible. Computationally, each of the three models require about the same run time (on a standard laptop computer), and run faster than computing the ground truth pde for few observables (4–32), but become slower for large numbers of internal variables (64) due to the need to solve linear systems involving the covariance matrix C i j .

Figure 2 shows the results for the external driving example. Panel A, again, shows the total variation distance from the ground truth for each model, but now as a function of time over the whole course of the simulation, and for the approximations with 16 and 64 internal variables. After an initial increase, the error of each of the methods is relatively stable in time. This mirrors the observation made in previous implementations of the STIV framework, where the approximation seemed to recover despite passing through regions in time where the approximation was quite poor [18]. Panel B shows the sinusoidal potential in blue and the external driving potential due to an optical trap in pink, along with a cartoon of the model setup. The remaining two panels (C and D) compare the approximate densities with 64 internal variables to the ground truth at two representative points in time. Although not shown here, it is worth noting that when fewer than 64 internal variables are used, including the energy greatly improves the approximation in this example. The external driving is slow enough that the density remains close to equilibrium at all times, and hence the energy is a very good observable.

5 Conclusions

We have shown that given a microscopic particle system governed by overdamped Langevin dynamics it is always possible to construct a thermodynamically consistent macroscopic (ensemble averaged) model with a gradient flow evolution for the internal variables. One strategy, which we’ve called the operator method, applicable for an arbitrary choice of approximate probability density of states, is presented but not pursued due to computational challenges of a general implementation. Alternative methods, restricted to quasi-equilibrium or exponential family approximations, can be implemented relatively simply. These latter models can either be formulated in terms of averaged observables or their dual internal variables, and can approximate continuous probability distributions with arbitrary accuracy.

In this work, we have shown the flexibility of this framework by considering Fourier basis functions as observables. Such basis functions provide an essentially system agnostic model, as the only way the approximate probability density of states depends on the system of study is through the system size. Since the Fourier basis functions can approximate any continuous function, in principle, one could use them to model systems of arbitrary complexity. However, this becomes impracticable both because the number of observables needed for a reasonable approximation would scale exponentially in the number of degrees of freedom in the system, and because the Fourier functions provide relatively little insight into the underlying structure of the system. Alternatively, if one can capture the salient features of the system using observables especially chosen for the problem, one can expect the STIV framework to produce accurate predictions with relatively little computational cost, as was the case in the previous works [18], [22]. The true utility of internal variable models, therefore, comes from the intelligent choice of a few highly descriptive observables. Thus, a key task for continuing work is that of finding descriptive observables for a given microscopic system (see [41], [42] for emerging efforts in this direction).

Corresponding author: Celia Reina, Department of Mechanical Engineering and Applied Mechanics, University of Pennsylvania, 220 South 33rd Street, 229 Towne Building, Philadelphia, PA, 19104, USA, E-mail: creina@seas.upenn.edu

Research ethics: Not applicable.
Informed consent: Not applicable.
Author contributions: All authors have accepted responsibility for the entire content of this manuscript and approved its submission. Travis Leadbetter: Conceptualization, Methodology, Software, Formal Analysis, Writing – Original Draft, Writing – Review & Editing. Prashant K. Purohit: Conceptualization, Writing – Review & Editing, Supervision. Celia Reina: Conceptualization, Methodology, Writing – Review & Editing, Supervision.
Use of Large Language Models, AI and Machine Learning Tools: None declared.
Conflict of interest: All authors state no conflict of interest.
Research funding: T. L. acknowledges funding support from the University of Pennsylvania Applied Mathematics and Computational Science graduate group. C. R. gratefully acknowledges support from NSF CAREER Award (CMMI-2047506). P. K. P. acknowledges partial support from NSF Grant DMR 2212162.
Data availability: Source code used to generate the data and figures displayed in this work is freely available via the GitHub repository celiareina/STIV-GradFlow.

Appendix A: The variational method of Eyink with external driving

In the original derivation of the variational method of Eyink [21], when the true solution to the Fokker–Planck equation ∂ t p ( x , t ) = L p ( x , t ) is approximated by a parameterized density of states p ̂ ( x , α ) , it was assumed that all of the time dependence of p ̂ is captured by the parameters α(t). However, if we wish to implement a quasi-equilibrium approximation p ̂ = exp ( − β e + α i ( t ) ϕ i ( x ) − f ( α , β ) ) and the energy depends on a time-dependent external driving protocol e = e(x, λ(t)), we must modify the variational method since the time dependence assumption has been broken.

Assume for fixed external driving protocol λ(t), we approximate the density of states via p ̂ ( x , α , λ ) (here β is assumed to be fixed), and approximate test function ψ ̂ ( x , α ) . The variational action is S ( α ) = ∫ 0 ∞ ∫ X ψ ̂ ∂ t − L p ̂ d x d t . Stationary points of the action are given by

0 = ∂ α i S ( α ) = ∫ 0 ∞ ∫ X ∂ α i ψ ̂ ∂ t p ̂ − ∂ α i p ̂ ∂ t ψ ̂ d x d t − ∂ α i ∫ 0 ∞ ∫ X ψ ̂ L p ̂ d x d t = ∫ 0 ∞ ∫ X ∂ α i ψ ̂ ∂ α j p ̂ α ̇ j + ∂ α i ψ ̂ ∂ λ p ̂ λ ̇ − ∂ α i p ̂ ∂ α j ψ ̂ α ̇ j − ∫ 0 ∞ ∂ α i L † ψ ̂ p ̂ d t = ∫ 0 ∞ [ ∂ ψ ̂ , ∂ u ̂ ] i j p ̂ α ̇ j + ∂ α i ψ ̂ ∂ λ u ̂ p ̂ λ ̇ − ∂ α i L † ψ ̂ p ̂ d t .

The proper dynamics are recovered by setting the integrand equal to zero.

We now assume the “linear ansatz” for the test function by introducing the paired set of parameters (α, γ) and write ψ ̂ = 1 + γ i ∂ α i u ̂ . Plugging this in to the dynamical equation above leads to two coupled equations

γ k ∂ α i L † ∂ α k u ̂ p ̂ = γ k ∂ α i ∂ α k u ̂ ∂ α j u ̂ − ∂ α j ∂ α k u ̂ ∂ α i u ̂ p ̂ α ̇ j − ∂ α i u ̂ ∂ α j u ̂ p ̂ γ ̇ j + γ k ∂ α i ∂ α k u ̂ ∂ λ u ̂ p ̂ λ ̇ L † ∂ α i u ̂ p ̂ = ∂ α i u ̂ ∂ α j u ̂ p ̂ α ̇ j + ∂ α i u ̂ ∂ λ u ̂ p ̂ λ ̇ .

The first equation is trivially solved by γ _i ≡ 0, whereas the second equation is independent of γ _i and gives the resulting dynamics we need.

Appendix B: Approximate thermodynamic quantities for the quasi-equilibrium approximation of a system with external driving

Since the inclusion of the energy in the approximate probability density of states breaks the assumption that p ̂ be independent of λ, the form of the approximate thermodynamic quantities must change as a result. Recall that the quasi-equilibrium approximation takes the form p ̂ = exp ( − β e ( x , λ ) + α i ϕ i ( x ) − f ( α , β , λ ) ) so that ∂ α i u ̂ = ϕ i − Φ ̂ i , a ̂ neq = e + u ̂ / β = ( α i ϕ i − f ) / β and, following the same steps as in the main text, L † ∂ α i u ̂ p ̂ = − 1 η β M i j α j where M i j = ∂ x ϕ i ∂ x ϕ j p ̂ . Plugging this into the dynamical equations derived from the previous section leads to

C i j α ̇ j + ( ϕ i − Φ ̂ i ) ∂ λ u ̂ p ̂ λ ̇ = − 1 η β M i j α j

where C i j = ( ϕ i − Φ ̂ i ) ( ϕ j − Φ ̂ j ) p ̂ . The approximate external force is defined as F ̂ ex ≡ ∂ λ e p ̂ which is no longer simply ∂ λ A ̂ neq . Now

∂ λ A ̂ neq = ∂ λ e p ̂ + ( e + u ̂ / β ) ∂ λ u ̂ p ̂ = F ̂ ex + 1 β α i ϕ i ∂ λ u ̂ p ̂ = F ̂ ex + 1 β ∂ λ ( α i Φ ̂ i )

which gives the first part of Eq. (11) of the main text. Since A ̂ neq = ( α i Φ ̂ i − f ) / β , A ̂ neq − α i Φ i / β = − f / β and we recover the second part.

To identify the form of the entropy production, we return to its definition from stochastic thermodynamics and its approximation via STIV. The incremental entropy production is the difference in the incremental work minus the incremental free energy, Tds ^tot = dw − da ^neq. Approximating these terms at the ensemble average level leads to the required relation T d d t S ̂ tot = d d t W ̂ − d d t A ̂ neq . The work rate is identified at d d t W ̂ = F ̂ ex λ ̇ and the rate of change in the non-equilibrium free energy writes d d t A ̂ neq = ∂ λ A ̂ neq λ ̇ + ∂ α A ̂ neq α ̇ . Inserting ∂ λ A ̂ neq = F ̂ ex + α i ( ϕ i − Φ ̂ i ) ∂ λ u ̂ p ̂ / β into the equation for T d d t S ̂ tot leaves (note ∂ λ u ̂ has mean zero due to normalization)

T d d t S ̂ tot = − α ̇ ∂ α A ̂ neq − α i β ( ϕ i − Φ ̂ i ) ∂ λ u ̂ p ̂ λ ̇ .

Since ∂ α i A ̂ neq = C i j α j / β , we can multiply the whole dynamical equation by −α _i/β to see

T d d t S ̂ tot = α i M i j α j η β 2

giving Eq. (12) in the main text.

Appendix C: Proof that exponential families and quasi-equilibrium distributions can approximate any continuous probability density on R n in total variation distance

Let p be a continuous probability density on R n and let u = log(p). Since p is integrable, for any fixed 0 < ϵ < 1, we can find a compact subset K ⊂ R n such that K c = R n \ K obeys ∫ K c p ( x ) d x < ϵ /2. Let K ′ = K ∩ { p ≥ ϵ 2 Vol ( K ) } . K′ is compact since p is continuous, and we know

∫ K ′ c p ( x ) d x = ∫ K c p ( x ) d x + ∫ K ∩ { p < ϵ / 2 Vol ( K ) } p ( x ) d x < ϵ / 2 + ϵ / 2 = ϵ .

u is bounded on K′, so by the Stone–Weierstrass theorem, we can find a polynomial q(x) such that ‖ u − q ‖ L ∞ ( K ′ ) = sup x ∈ K ′ | u ( x ) − q ( x ) | < ϵ . As e ^q(x) = e ^u(x) e ^q(x)−u(x) ≤ e ^u(x)e^ϵ, Z _q ≡ ∫ _K′ e ^q(x)dx < ∞ and we can define p ̂ = 1 x ∈ K ′ e q ( x ) Z q . In fact, we have e ^−ϵ(1 − ϵ) ≤ Z _q ≤ e^ϵ. Now

d TV ( p , p ̂ ) = 1 2 ∫ | p ( x ) − p ̂ ( x ) | d x = 1 2 ∫ K ′ | p ( x ) − p ̂ ( x ) | d x + 1 2 ∫ K ′ c p ( x ) d x ≤ 1 2 ∫ K ′ e u 1 − e q − u Z q d x + ϵ 2 ≤ 1 2 ‖ 1 − e q − u Z q ‖ L ∞ ( K ′ ) ∫ K ′ p d x + ϵ 2 ≤ 1 2 e 2 ϵ − 1 + ϵ → 0

as ϵ → 0 since ‖ 1 − e q − u Z q ‖ L ∞ ( K ′ ) ≤ max e 2 ϵ 1 − ϵ − 1,1 − e − 2 ϵ = e 2 ϵ 1 − ϵ − 1 .

An analogous proof shows that the same is true for quasi-equilibrium distributions. Assuming that the energy is continuous as well, we simply choose K′ as above, q such that ‖ u + β e − q ‖ L ∞ ( K ′ ) < ϵ , define Z _q = ∫ _K′e^−βe+qdx. And let p ̂ = 1 x ∈ K ′ e − β e + q Z q . The same bounds hold for Z _q so we again have

d TV ( p , p ̂ ) ≤ 1 2 e 2 ϵ − 1 + ϵ → 0 .

Appendix D: Details of the numerical implementation of STIV for examples

As discussed in the main text, the energy used to study relaxation to equilibrium comes from a model of the energy landscape of a protein diffusing on a strand of DNA. The original experimental data was published by [36], and the particular choice of model energy landscape was taken from [37]. The parameters used in Eq. (13) are c ₀ = −0.87817, c ₁ = 0.04736, c ₂ = 0.60955, c ₃ = 1.25408. A bounding box of the form (where x has units of base pairs)

Bounding Box ( x ) = 100 e − x − 40 1 ( x < − 25 ) + ( x + 10 ) 30 100

was used to constrain the simulated proteins to x ∈ [−40, 20] (that is, the initial range of the model energy landscape in [37]) and to insure that the minima at x ≈ −30 was the left most minima. This range is translated to [0,60] in the video in Appendix E. Since the energy has units of k _B T, we used β = 1. We chose η = 0.83 which is artificially large to accelerate all simulations.

For the external driving example, we chose parameters l = 1.0, k = 1.6287, β = 1.0, and η = 0.8298. The external protocol was chosen as λ(t) = 0.5t.

All models (both STIV and the ground truth) were implemented on a discrete grid of 2,000 equally spaced points (between [−39, 21.1] for the relaxation example and [−2, 8] for the external driving example). The ground truth was obtained from discretizing the true pde, ∂ t p = L p , and solving the resulting ode using the Dormand–Prince adaptive Runge–Kutta scheme with variable step size to ensure the relative error remained less than 0.001. For the STIV models, the energy and observables were computed on the grid points e i = e ( x i ) ∈ R 2000 and ϕ i j = ϕ j ( x i ) ∈ R 2000 × d where d is the number of internal variables used. p ̂ i = p ̂ ( x i , α ) was then approximated as

p ̂ i = 1 Δ x exp ( − β e i + ϕ i j α j − f ( α , β ) )

with

f ( α , β ) = log ∑ i = 1 2,000 exp ( − β e i + ϕ i j α j ) .

Any expectations in the dynamical equations were then approximated using the discrete approximate density of states

g p ̂ ≈ ∑ i = 1 2000 g ( x i ) p ̂ i Δ x .

The exponential family models which did not include the energy were obtained the same way but with β fixed to zero.

Some numerical difficulties arose while inverting the covariance matrix C i j . Thus, C i j − 1 was further approximated using the pseudo-inverse C i j + truncating singular values (of C i j ) falling below 10⁻⁸. Various truncation levels (10⁻² through 10⁻¹⁵) were tested to ensure this choice did not impact numerical accuracy.

Appendix E: Full video of examples

The video provided in the supplementary materials shows a comparison between the true probability density of states in dark blue and the approximate densities using the quasi-equilibrium method (green), the exponential family approximation with dynamics given by the variational method (equivalent to moment closure, yellow), and the exponential family approximation with dynamics given by a gradient flow (salmon). The first part of the video shows, sequentially, the approximation for the model of proteins diffusing on DNA using 4, 8, 16, 32, and 64 internal variables (the legend |ϕ| = N in each part shows the number of internal variables/observables). The second half of the video shows the approximation for the externally driven particle in a periodic landscape for 4, 8, 16, 32, and 64 internal variables.

References

[1] D. Jou, J. Casas-Vázquez, and G. Lebon, “Extended irreversible thermodynamics,” in Extended Irreversible Thermodynamics, Springer, 1996, pp. 41–74.10.1007/978-3-642-97671-1_2Search in Google Scholar

[2] G. Lebon, D. Jou, and J. Casas-Vázquez, Understanding Non-Equilibrium Thermodynamics, vol. 295, Berlin, Heidelberg, Springer, 2008.10.1007/978-3-540-74252-4Search in Google Scholar

[3] B. Coleman, Thermodynamics of Materials with Memory, 1st ed. Vienna, Springer, 1971.10.1007/978-3-7091-2951-7_1Search in Google Scholar

[4] C. Truesdell, “Historical introit the origins of rational thermodynamics,” in Rational Thermodynamics, Springer, 1984, pp. 1–48.10.1007/978-1-4612-5206-1_1Search in Google Scholar

[5] G. A. Maugin and W. Muschik, “Thermodynamics with internal variables. Part i. General concepts,” J. Non-Equilib. Thermodyn., vol. 19, no. 1, pp. 217–249, 1994. https://doi.org/10.1515/jnet.1994.19.3.217.Search in Google Scholar

[6] M. Grmela and H. C. Öttinger, “Dynamics and thermodynamics of complex fluids. i. Development of a general formalism,” Phys. Rev. E, vol. 56, no. 6, pp. 6620–6632, 1997, https://doi.org/10.1103/physreve.56.6620.Search in Google Scholar

[7] A. Mielke, “Formulation of thermoelastic dissipative material behavior using generic,” Continuum Mech. Thermodyn., vol. 23, no. 3, pp. 233–256, 2011, https://doi.org/10.1007/s00161-010-0179-0.Search in Google Scholar

[8] M. Pavelka, V. Klika, and M. Grmela, Multiscale Thermo-Dynamics: Introduction to GENERIC, Berlin/Boston, Walter de Gruyter GmbH & Co KG, 2018.10.1515/9783110350951Search in Google Scholar

[9] F. Garofalo, “An empirical relation defining the stress dependence to minimum creep rate in metals,” Trans. Metall. Soc. AIME, vol. 227, no. 1, p. 351, 1963.Search in Google Scholar

[10] J.-L. Chaboche, “Continuous damage mechanics – a tool to describe phenomena before crack initiation,” Nucl. Eng. Des., vol. 64, no. 2, pp. 233–247, 1981, https://doi.org/10.1016/0029-5493(81)90007-8.Search in Google Scholar

[11] J. Lemaitre, “A continuous damage mechanics model for ductile fracture,” J. Eng. Mater. Technol., vol. 107, no. 1, pp. 83–89, 1985, https://doi.org/10.1115/1.3225775.Search in Google Scholar

[12] Z. P. Bazant, “Mechanics of distributed cracking,” Appl. Mech. Rev., vol. 39, no. 5, pp. 675–705, 1986, https://doi.org/10.1115/1.3143724.Search in Google Scholar

[13] H. C. Öttinger, Beyond Equilibrium Thermodynamics, vol. 5, Hoboken, Wiley-Interscience, 2005.10.1002/0471727903Search in Google Scholar

[14] R. Zwanzig, Nonequilibrium Statistical Mechanics, New York, Oxford University Press, 2001.10.1093/oso/9780195140187.001.0001Search in Google Scholar

[15] H. C. Öttinger, “General projection operator formalism for the dynamics and thermodynamics of complex fluids,” Phys. Rev. E, vol. 57, no. 2, pp. 1416–1420, 1998, https://doi.org/10.1103/physreve.57.1416.Search in Google Scholar

[16] A. Montefusco, M. A. Peletier, and H. C. Öttinger, “A framework of nonequilibrium statistical mechanics. ii. Coarse-graining,” J. Non-Equilib. Thermodyn., vol. 46, no. 1, pp. 15–33, 2021, https://doi.org/10.1515/jnet-2020-0069.Search in Google Scholar

[17] A. Mielke, M. A. Peletier, and J. Zimmer, “Deriving a generic system from a hamiltonian system,” Arch. Ration. Mech. Anal., vol. 249, no. 1, 2025, Art. no. 62. https://doi.org/10.1007/s00205-025-02119-7.Search in Google Scholar

[18] T. Leadbetter, P. K. Purohit, and C. Reina, “A statistical mechanics framework for constructing nonequilibrium thermodynamic models,” PNAS Nexus, vol. 2, no. 12, 2023, Art. no. pgad417, https://doi.org/10.1093/pnasnexus/pgad417.Search in Google Scholar

[19] U. Seifert, “Stochastic thermodynamics, fluctuation theorems and molecular machines,” Rep. Prog. Phys., vol. 75, no. 12, 2012, Art. no. 126001, https://doi.org/10.1088/0034-4885/75/12/126001.Search in Google Scholar PubMed

[20] L. Peliti and S. Pigolotti, Stochastic Thermodynamics: An Introduction, Princeton, Princeton University Press, 2021.Search in Google Scholar

[21] G. L. Eyink, “Action principle in nonequilibrium statistical dynamics,” Phys. Rev. E, vol. 54, no. 4, pp. 3419–3435, 1996, https://doi.org/10.1103/physreve.54.3419.Search in Google Scholar PubMed

[22] T. Leadbetter, P. K. Purohit, and C. Reina, “A statistical mechanics derivation and implementation of non-conservative phase field models for front propagation in elastic media,” J. Mech. Phys. Solids, vol. 204, no. 1, p. 2025, https://doi.org/10.1016/j.jmps.2025.106240.Search in Google Scholar

[23] R. Abeyaratne, S. Dharmaravan, G. Saccomandi, and G. Tomassetti, “Using stochastic thermodynamics with internal variables to capture orientational spreading in cell populations undergoing cyclic stretch,” arXiv preprint arXiv:2507.15694, 2025.Search in Google Scholar

[24] R. Jordan, D. Kinderlehrer, and F. Otto, “The variational formulation of the Fokker–Planck equation,” SIAM J. Math. Anal., vol. 29, no. 1, pp. 1–17, 1998, https://doi.org/10.1137/s0036141096303359.Search in Google Scholar

[25] B. Coleman and W. Noll, “The thermodynamics of elastic materials with heat conduction and viscosity,” Arch. Ration. Mech. Anal., vol. 13, no. 1, pp. 167–178, 1963, https://doi.org/10.1007/bf01262690.Search in Google Scholar

[26] E. Hairer, M. Hochbruck, A. Iserles, and C. Lubich, “Geometric numerical integration,” Oberwolfach Rep., vol. 3, no. 1, pp. 805–882, 2006, https://doi.org/10.4171/owr/2006/14.Search in Google Scholar

[27] B. Gladman, M. Duncan, and J. Candy, “Symplectic integrators for long-term integrations in celestial mechanics,” Celest. Mech. Dyn. Astron., vol. 52, no. 3, pp. 221–240, 1991, https://doi.org/10.1007/bf00048485.Search in Google Scholar

[28] H. Yoshida, “Recent progress in the theory and application of symplectic integrators,” Celest. Mech. Dyn. Astron., vol. 56, nos. 1–2, pp. 27–43, 1993, https://doi.org/10.1007/bf00699717.Search in Google Scholar

[29] B. Leimkuhler and S. Reich, Simulating Hamiltonian Dynamics, Number 14, Cambridge, Cambridge University Press, 2004.10.1017/CBO9780511614118Search in Google Scholar

[30] M. Kobilarov, K. Crane, and M. Desbrun, “Lie group integrators for animation and control of vehicles,” ACM Trans. Graph. (TOG), vol. 28, no. 2, pp. 1–14, 2009, https://doi.org/10.1145/1516522.1516527.Search in Google Scholar

[31] H. C. Öttinger, “GENERIC integrators: structure preserving time integration for thermodynamic systems,” J. Non-Equilib. Thermodyn., vol. 43, no. 2, pp. 89–100, 2018, https://doi.org/10.1515/jnet-2017-0034.Search in Google Scholar

[32] L. Bronstein and H. Koeppl, “A variational approach to moment-closure approximations for the kinetics of biomolecular reaction networks,” J. Chem. Phys., vol. 148, no. 1, 2018, Art. no. 014105, https://doi.org/10.1063/1.5003892.Search in Google Scholar PubMed

[33] B. Turkington, “An optimization principle for deriving nonequilibrium statistical models of hamiltonian dynamics,” J. Stat. Phys., vol. 152, no. 1, pp. 569–597, 2013. https://doi.org/10.1007/s10955-013-0778-9.Search in Google Scholar

[34] H. Risken, The Fokker–Planck Equation, Volume 18 of Springer Series in Synergetics, Berlin Heidelberg, Springer-Verlag, 1996.10.1007/978-3-642-61544-3Search in Google Scholar

[35] E. L. Lehmann and G. Casella, Theory of Point Estimation, New York, Springer Science & Business Media, 2006.Search in Google Scholar

[36] S. Kim et al.., “Probing allostery through dna,” Science, vol. 339, no. 6121, pp. 816–819, 2013, https://doi.org/10.1126/science.1229223.Search in Google Scholar PubMed PubMed Central

[37] J. Singh and P. K. Purohit, “Elasticity as the basis of allostery in dna,” J. Phys. Chem. B, vol. 123, no. 1, pp. 21–28, 2018, https://doi.org/10.1021/acs.jpcb.8b07501.Search in Google Scholar PubMed PubMed Central

[38] M. S. Pinsker, Information and Information Stability of Random Variables and Processes, San Francisco, Holden Day, 1964.Search in Google Scholar

[39] J. Bretagnolle and C. Huber, “Estimation des densités: risque minimax,” Z. Wahrscheinlichkeitstheor. Verw. Geb., vol. 47, no. 1, pp. 119–137, 1979. https://doi.org/10.1007/bf00535278.Search in Google Scholar

[40] C. L. Canonne, “A short note on an inequality between kl and tv,” arXiv preprint arXiv:2202.07198, 2022.Search in Google Scholar

[41] B. Liu, E. Ocegueda, M. Trautner, A. M. Stuart, and K. Bhattacharya, “Learning macroscopic internal variables and history dependence from microscopic models,” J. Mech. Phys. Solids, vol. 178, no. 1, 2023, Art. no. 105329, https://doi.org/10.1016/j.jmps.2023.105329.Search in Google Scholar

[42] W. Qiu, S. Huang, and C. Reina, “Bridging statistical mechanics and thermodynamics away from equilibrium: a data-driven approach for learning internal variables and their dynamics,” J. Mech. Phys. Solids, vol. 203, 2025, Art. no. 106211, https://doi.org/10.1016/j.jmps.2025.106211.Search in Google Scholar

Supplementary Material

This article contains supplementary material (https://doi.org/10.1515/jnet-2025-0071).

Video 1

Received: 2025-06-09

Accepted: 2025-09-26

Published Online: 2025-10-16

This work is licensed under the Creative Commons Attribution 4.0 International License.

https://doi.org/10.1515/jnet-2025-0071

Keywords for this article

stochastic thermodynamics; thermodynamics with internal variables; gradient flow

Creative Commons

BY 4.0