1 Introduction↩︎

Stochastic thermodynamics (ST) is the framework developed to systematically understand the thermodynamics associated with systems where fluctuations play an important physical role in the dynamics, e.g. mesoscopic systems and critical phenomena. The need for this formalism is motivated especially for nonequilibrium systems, where the conventional tools of statistical mechanics do not apply. This framework builds on the extensions of the first and second laws to the stochastic regime, providing a description at the individual trajectory level, as well as universal constraints on the distributions that hold beyond linear response theory. Under the assumption of time-scale separation, the observable (slow) and thermal (fast) degrees of freedom (DoFs) are typically modeled by Markovian dynamics, such as master equations (for discrete-time systems) or Langevin equations (for continuous-time systems) [1].

For continuous-time systems, which are the focus of this paper, the Langevin equation must be complemented by an appropriate fluctuation-dissipation relation (FDR) to ensure thermodynamic consistency. Such an FDR will be shown to provide a well-defined interpretation of entropy production by connecting the first and second laws, as well as guaranteeing the existence of an equilibrium steady state for any closed system. However, there is no systematic approach for determining such relations in a general setting, and the Einstein relation [2] is often assumed by default [1], justified by the requirement that the corresponding Fokker–Planck equation is satisfied by the Boltzmann distribution in the steady state [3], [4]. Moreover, thermodynamically consistent modeling of complex systems remains challenging when accounting for stochastic dynamics, particularly in nonequilibrium regimes, such as active matter systems [5], [6]. We emphasize the meaning of FDRs as a generalization of the Einstein relation, i.e. the relation between dissipative forces and noise strength, and is not to be identified with fluctuation-dissipation or fluctuation-response theorems as in how the system responds to an external perturbation.

The goal of this paper is to present a novel approach to the thermodynamically consistent modeling of stochastic finite-dimensional, continuous-time systems. To this end, we extend the variational formulation of nonequilibrium thermodynamics presented in Ref. [7] to the stochastic regime.

The variational approach in Ref. [7] generalizes classical mechanics to include irreversible processes by extending the Lagrange-d’Alembert principle to incorporate entropy production as a nonlinear nonholonomic constraint, with phenomenological laws embedded directly into the constraint structure. A key insight of this framework is the introduction of thermodynamic displacements, conjugates to the thermal variables, which provide a natural definition of the variational constraint and generalize earlier notions of thermal displacement. The resulting variational structure yields a complete set of evolution equations, formulated as a system of ordinary differential equations (ODEs) for discrete systems, enabling the determination of all system variables at all times. It is important to note that this variational formulation consistently recovers Hamilton’s principle in the absence of irreversible processes, and—as a macroscopic description—assumes local thermodynamic equilibrium.

By extending this framework to the stochastic regime, we generalize it to describe systems driven far from equilibrium. The proposed variational principle for stochastic thermodynamics is constrained by the first and second laws of thermodynamics, under the assumption of time-scale separation. No steady-state condition is imposed, nor is the thermal bath required to be infinite. Instead, results from statistical mechanics are recovered in the appropriate limits, while the formulation reproduces macroscopic nonequilibrium thermodynamics in the vanishing-noise limit. Moreover, as shown in subsequent sections, this approach naturally extends key features of stochastic thermodynamics, including an individual trajectory level description of both configurational and thermal variables, and fluctuation theorems (FTs) in the thermodynamic phase space. Remarkably, once thermodynamic consistency is enforced via the FDRs, the equations for entropy production derived from the variational approach assume the same form as in standard ST, up to a reformulation in an extended phase space.

The paper is organized as follows. Section 2 provides an overview of the framework of ST, which can be skipped by readers already familiar with the topic, except perhaps for subsection 2.4. Section 3 introduces the variational principle for ST and illustrates it with examples of (finite-dimensional) closed and open systems. Section 4 offers a discussion of our findings, and Section 5 summarizes the main conclusions of this study. Additional analytical derivations are organized into five appendices at the end.

2 Stochastic thermodynamics↩︎

2.1 Representations of Langevin processes↩︎

The stochastic dynamics of continuous-time systems have three different but equivalent descriptions [1], [8]: the Langevin equation, the Fokker-Planck equation and the path integral formulation.

In the general case of a \(d\)-dimensional particle system, the Langevin equation for underdamped dynamics can be written as \((\nu=1,...,d)\) [9]: \[\begin{align} &m\dot{x}^\nu = p^\nu, \tag{1}\\ & \dot{p}^\nu = F^\nu(\mathbf{x},\mathbf{p},\lambda) + g^\nu _k(\mathbf{x},\mathbf{p};t)\circ\zeta^k (t),\tag{2} \end{align}\] where \((\mathbf{x},\mathbf{p})\) correspond to position and momentum respectively, \(F^\nu\) is a generalized force, \(\lambda (t)\) denotes a time-dependent driving, and \(\zeta^\nu\) is a gaussian white noise verifying \[\label{eq:gaussian95noise} \begin{eqnarray} &&\langle \zeta^k(t) \rangle=0,\\ &&\langle \zeta^k(t)\zeta^\ell(t') \rangle=2\delta^{k\ell}\delta(t-t'). \end{eqnarray}\tag{3}\] The product \(\circ\) denotes the Stratonovich convention (SC), which will be applied throughout this paper due to its suitability for modeling physical processes [8]. Repeated indices are summed over according to Einstein’s convention. The noise amplitude contribution \(g(\mathbf{x},\mathbf{p};t)\) allows for the possibility of state-dependent noise. The corresponding Fokker-Planck equation (FP) is [10]: \[\begin{align} \label{eq:FP} \partial_t \mathcal{P}&=& -\partial_{x^\nu} \left(\frac{p^\nu}{m}\mathcal{P}\right) -\partial_{p^\nu}\left[\left(F^\nu -g^\nu_k \partial_{p^\mu}g^\mu_k\right)\mathcal{P} -g^\nu_k g^\mu_k \partial_{p^\mu}\mathcal{P}\right] \nonumber\\ &=& -\nabla_{\mathbf{x}}\cdot\mathbf{j}_{x}-\nabla_{\mathbf{p}}\cdot\mathbf{j}_{p} , \end{align}\tag{4}\] where \(\mathcal{P}(\mathbf{x},\mathbf{p},t)\) denotes the probability density function (PDF) of a certain configuration \((\mathbf{x},\mathbf{p})\) at time \(t\), and \(\mathbf{j}\) the probability current. While the Langevin equation 1 2 provides dynamical information of an individual trajectory, the FP equation 4 provides the ensemble level information.

An additional representation which is applied in ST is the path integral (PI) formulation. Similarly as for the Langevin equation, this description allows assigning thermodynamic observables to individual trajectories. Such path integral is the analog of the propagator \(K\) in Quantum Mechanics, where the analog of the Hamiltonian would in this case correspond to the Fokker-Planck operator [11]: \[\begin{eqnarray} &&\mathcal{P}(\mathbf{z},t)=\int _{\mathbb{R}^{2d}}d\mathbf{z}_0 \, K(\mathbf{z},t;\mathbf{z}_0,t_0) \mathcal{P}_0(\mathbf{z}_0),\label{eq:propagator}\\ &&K(\mathbf{z}_t,t;\mathbf{z}_0,t_0) = \int_{\mathbf{z}_0}^{\mathbf{z}_t} \mathcal{D}\mathbf{z}\,P[\mathbf{z}|\mathbf{z}_0] . \end{eqnarray}\tag{5}\] For convenience, we denote the full phase-space variable as \(\mathbf{z} = (\mathbf{x}, \mathbf{p})\). The initial and final times are denoted by \(t_0\) and \(t\), respectively, with corresponding configurations \(\mathbf{z}_0\) and \(\mathbf{z}_t\). The path-probability can be defined as \[P[\mathbf{z}|\mathbf{z}_0]=\mathcal{N}\{\mathbf{z}_n\}\, \exp (-\mathcal{A}[\mathbf{z}]),\] where \(\mathcal{A}[\mathbf{z}]\) is the action functional and \(\mathcal{N}\{\mathbf{z}_n\}\) the normalization factor defined for a discretized temporal grid \(\{t_n = \epsilon n\}_{n=0}^N\). Since to each Langevin equation corresponds a unique FP equation, the PDF \(\mathcal{P}(\mathbf{z},t)\) is uniquely defined for a given interpretation (Stratonovich/Ito). However, this only uniquely defines \(K(\mathbf{z}_t,t;\mathbf{z}_0,t_0)\) in the continuous limit, allowing the freedom to define different normalizations \(\mathcal{N}\{\mathbf{z}_n\}\) or actions \(\mathcal{A}[\mathbf{z}]\) as long as Eq. 5 is fulfilled [12].

As a result, different expressions for the path integral (PI) corresponding to the same interpretation of the noise can be found in the literature [5], [9], [11], [13]–[17]. However, this ambiguity disappears in the context of ST when the contribution of the normalization \(\mathcal{N}\{\mathbf{z}_n\}\) is properly taken into account. In this case, the results at the level of individual trajectories obtained via the PI formalism coincide for its different forms, as well as with those derived from the FP formalism—which is uniquely defined—when ensemble averages are computed, as required for consistency.

A convenient expression for the PI in the context of ST is the Onsager-Machlup (OM) action [5], [13], [14]. For a discretized temporal grid \(\{t_n = \epsilon n\}_{n=0}^N\), the propagator and the action¹ are given, respectively, by: \[\label{eq:PI95multid} \begin{align} & K(\mathbf{z}_t,t;\mathbf{z}_0,t_0) = \int_{\mathbf{z}_0}^{\mathbf{z}_t}\prod_{n=1}^{N-1}\prod_\nu\frac{dz_n^\nu \delta(m\dot{x}^\nu - p^\nu)}{\det(g^{\nu\mu}_{n+1/2})\sqrt{4\pi\epsilon}}\, e^{-\mathcal{A}[\mathbf{z}]},\\ &\mathcal{A}[\mathbf{z}] =\int_{t_0}^{t}d\tau \left[ \frac{1}{2}\partial_{p^\mu} F^\mu +\frac{1}{4}\left(\partial_{p_\nu} g^\mu_i \partial_{p_\mu} g^\nu_i -\partial_{p_\mu} g^\mu_i \partial_{p_\nu} g^\nu_i \right)\right.\\ &\left. \qquad + \frac{1}{4}(\dot{p}^\mu-F^\mu + g_{i}^\mu\partial_{p^\alpha} g_{i}^\alpha)D_{\mu\nu}(\dot{p}^\nu-F^\nu + g_{i}^\nu\partial_{p^\alpha} g_{i}^\alpha)\right], \end{align}\tag{6}\] with \(D^{\nu\mu}= g_{i}^\nu g_{i}^\mu\) denoting the covariance matrix and \(D_{\mu\nu}=(D^{\nu\mu})^{-1}\) its inverse (provided it exists), such that \(D_{\mu\nu}D^{\nu\alpha}=\delta_\mu^\alpha\).

We have already accounted for Eq. 1 via the Dirac delta term, and we denote the constrained integral over paths as \[\mathcal{D}\mathbf{z}=\mathcal{D}\mathbf{x}\mathcal{D}\mathbf{p}\delta^{(d)} (\mathbf{p}-m\mathbf{\dot{x}}).\] This contribution can be understood as the \(D\rightarrow 0\) limit of \(\mathbf{p}=m \mathbf{\dot{x}} + \boldsymbol{\xi}\) with \(\langle \xi^i(t)\xi^j(t') \rangle=2D\delta^{ij}\delta(t-t')\). Due to addition of actions, the PI over a single time step includes the contribution: \[\begin{align} P(\mathbf{z}_{n+1} | \mathbf{z}_n)&\propto& \lim_{D\rightarrow0}\frac{1}{(4\pi D\varepsilon)^{d/2}}\exp \left\{ -\varepsilon\frac{||m \mathbf{\dot{x}} - \mathbf{p}||^2}{4D}\right\} \nonumber\\ &=&\delta^{(d)}(\mathbf{p}-m\mathbf{\dot{x}}), \end{align}\] where the last equality corresponds to the definition of the Dirac delta function as a limit distribution [18].

The advantage of the OM form lies in the fact that its associated measure is symmetric under time reversal. This symmetry arises because the points in the time-reversed trajectory, defined by \(\hat{x}_{n+a} = x_{n+(1-a)}\), coincide with those in the forward trajectory when using the SC, i.e. \(a = 1/2\). As a result, only the terms in the action contribute to the entropy production. In contrast, other expressions of the path integral—such as those in Refs. [11], [16]—involve measures that are not symmetric under time reversal, either due to asymmetric discretization schemes or the presence of terms that are odd under time reversal. These asymmetries must therefore be explicitly accounted for when computing entropy production to avoid inconsistencies [5]. However, we emphasize that these alternative formulations yield covariant actions, in contrast to the OM form [11].

2.2 Stochastic entropy production↩︎

In stochastic thermodynamics (ST), entropy production (EP) is defined as the breaking of time-reversal symmetry (TRSB), and is quantitatively measured by the Kullback-Leibler (KL) divergence—or statistical distance—between the probability distributions of a forward process and its time-reversed counterpart [19]–[21]: \[\label{eq:KL} \langle \Delta s_{tot} \rangle \coloneq D_{KL}[P[\mathbf{z}]||P[\hat{\mathbf{z}}]]=\int_{\mathbf{z}_0}^{\mathbf{z}_t} \mathcal{D}\mathbf{z} P[\mathbf{z}]\text{ln}\frac{P[\mathbf{z}]}{P[\hat{\mathbf{z}}]},\tag{7}\] where—setting \(t_0 = 0\) without loss of generality—\(\mathbf{z}(\tau)=(\mathbf{x}(\tau),\mathbf{p}(\tau))\) denotes the forward trajectory, and \(\hat{\mathbf{z}}(\tau)=(\mathbf{x}(t-\tau),\mathbf{p}(t-\tau))\) the backward or time-reversed trajectory. In the presence of an externally controlled, time-dependent protocol \(\lambda(\tau)\), which reverses as \(\hat{\lambda}(\tau) = \lambda(t-\tau)\), the path probability is then given by \(P[\mathbf{z},\lambda]\) with reverse \(P[\hat{\mathbf{z}},\hat{\lambda}]\). Odd variables are those that change sign under time reversal, such as momentum or velocity, while even variables remain unchanged.

The KL divergence is an invariant (under non-singular transformations), nonnegative flat operator [22], such that the definition 7 already incorporates the second law of ST \(\langle \Delta s_{tot} \rangle \geq 0\). Moreover, as a consequence of this definition, the integral fluctuation theorem (IFT) follows [23]: \[\label{eq:IFT} \langle e^{-\Delta s_{tot}}\rangle = \int \mathcal{D}\mathbf{z} P[\mathbf{z}]\frac{P[\hat{\mathbf{z}}]}{P[\mathbf{z}]}= \int \mathcal{D}\mathbf{z} P[\hat{\mathbf{z}}] = 1\tag{8}\] where the last equality is drawn from the normalization condition and the fact that summing over all paths do not distinguish between forward and backward trajectories [1].

The EP for an individual trajectory is obtained by undoing the path average: \[\begin{align} &\Delta s_{tot} = \text{ln}\frac{P[\mathbf{z}]}{P[\hat{\mathbf{z}}]} = \Delta s_m + \Delta s, \tag{9}\\ &\Delta s_m \coloneq \text{ln}\frac{P[\mathbf{z}|\mathbf{z}_0 ]}{P[\hat{\mathbf{z}}|\mathbf{z}_t ]},\quad \Delta s\coloneq \text{ln}\frac{\mathcal{P}(\mathbf{z}_0, t_0)}{\mathcal{P}(\mathbf{z}_t,t)} . \tag{10} \end{align}\] The first term in Eq. 10 represents the stochastic EP in the medium \((\Delta s_m )\), while the second term accounts for the stochastic entropy change of the system \((\Delta s )\). The probability distributions \(\mathcal{P}\) are governed by the corresponding Fokker-Planck equation. Upon taking the ensemble average of the system contribution, one recovers the Shannon entropy \(\mathcal{S}\) [1]: \[\label{eq:stoch95shannon} \begin{eqnarray} s(t)&&= - \ln \mathcal{P}(\mathbf{z}(t),t),\\ \mathcal{S}(t) &&\coloneq \langle s(t) \rangle = -\int_{\mathbb{R}^{2d}} d\mathbf{z}\,\mathcal{P}(\mathbf{z},t)\ln \mathcal{P}(\mathbf{z},t) . \end{eqnarray}\tag{11}\] Taking the time derivative of the total entropy production yields the stochastic entropy production rate (EPR): \[\begin{align} \label{eq:dot95s95tot} \dot{s}_{tot}(t) &=& \dot{s} (t) + \dot{s}_m (t). \end{align}\tag{12}\] We illustrate the EPR by applying the framework to Eq. 2 . First, we define even and odd contributions to the force \(F^\nu_\pm=(F^\nu\pm \hat{F}^\nu)/2\), where the time-reversed force is defined as \(\hat{F}^\nu =F^\nu(\mathbf{x}(u),-\mathbf{p}(u),\lambda(u))\) for \(u=t-\tau\). Then, we can identify the reversible and dissipative contributions to the probability current along the momentum dimension: \[\begin{eqnarray} && j_p^\nu = j^\nu_r + j^\nu_d ,\\ && j^\nu_r = F^\nu_+ \mathcal{P} ,\\ && j^\nu_d = (F^\nu_- -g^\nu_k \partial_{p_\mu}g^\mu_k)\mathcal{P}-D^{\nu\mu} \partial_{p_\mu}\mathcal{P}. \end{eqnarray}\] The system EPR yields: \[\label{eq:stoch95s} \begin{align} \dot{s}&=-\frac{d}{dt}\ln \mathcal{P}= -\frac{\partial_t \mathcal{P}}{\mathcal{P}}-\frac{\partial_{x_\mu} \mathcal{P}}{\mathcal{P}}\frac{p^\mu}{m}-\frac{\partial_{p_\mu} \mathcal{P}}{\mathcal{P}}\dot{p}^\mu \\ &= -\frac{\partial_t \mathcal{P}}{\mathcal{P}}-\frac{\partial_{x_\mu} \mathcal{P}}{\mathcal{P}}\frac{p^\mu}{m}+ \dot{p}^\mu D_{\mu\nu}\left( \frac{j^\nu_d}{\mathcal{P}}-(F^\nu_- -g^\nu_k \partial_{p_\alpha}g^\alpha_k)\right) , \end{align}\tag{13}\] where we applied \[\partial_{p_\mu}\mathcal{P}= D_{\mu\nu} \left((F^\nu_- -g^\nu_k \partial_{p_\alpha}g^\alpha_k)\mathcal{P}-j^\nu_d\right).\] The medium EP is derived by applying the PI 6 in Eq. 10 : \[\begin{align} \Delta s_m &= \text{ln}\frac{P[\mathbf{z}|\mathbf{z}_0 ]}{P[\hat{\mathbf{z}}|\mathbf{z}_t ]}=\text{ln}\frac{\mathcal{N}\{\mathbf{z}_n\}\, \exp\{-\mathcal{A}\}}{\mathcal{N}\{\mathbf{\hat{z}}_n\}\, \exp\{-\hat{\mathcal{A}}\}}=\hat{\mathcal{A}}-\mathcal{A} \\ &= \int_{t_0}^t d\tau \left( (\dot{p}^\mu - F^\mu_+) D_{\mu\nu} (F^\nu_- -g^\nu_k \partial_{p_\alpha}g^\alpha_k)-\partial_{p_\mu}F^\mu_+\right). \end{align}\] Here, the backward action is given by \(\hat{\mathcal{A}}[\mathbf{z},\lambda] = \mathcal{A}[(\mathbf{x}(u), -\mathbf{p}(u), \lambda(u))]\), with \(u=t-\tau\). Since the integration limits coincide under this change of variables, we can relabel the dummy index \(du = d\tau\), allowing direct subtraction of the integrands [24]. Taking the time derivative gives the medium EPR (see [25] for a discussion of time integrals of Markovian variables): \[\label{eq:stoch95s95m} \dot{s}_m = (\dot{p}^\mu - F^\mu_+) D_{\mu\nu} (F^\nu_- -g^\nu_k \partial_{p_\alpha}g^\alpha_k)-\partial_{p_\mu}F^\mu_+ ,\tag{14}\] in agreement with [26] when restricting Eq. 14 to additive noise. Finally, adding both contributions gives the total EPR: \[\label{eq:total95EPR} \begin{align} \dot{s}_{tot} &= \dot{s} + \dot{s}_m = -\frac{\partial_t \mathcal{P}}{\mathcal{P}}-\frac{\partial_{x_\mu} \mathcal{P}}{\mathcal{P}}\frac{p^\mu}{m}+\frac{ D_{\mu\nu}j^\nu_d}{\mathcal{P}}(\dot{p}^\mu-F^\mu_+)\\ &\phantom{=} - \frac{F^\mu_+}{\mathcal{P}} \partial_{p_\mu}\mathcal{P}-\partial_{p_\mu}F^\mu_+ . \end{align}\tag{15}\] It can be seen how the stochastic total EPR is not constrained to be nonnegative, since the second law applies on the ensemble level rather than the individual trajectory level.

Computing the average EPR requires a two-step process, involving an average over all trajectories which are at time \(t\) in a state \(\mathbf{z}\), followed by an ensemble average [1], [23]. The first step involves the integral over paths which gives \(\langle \dot{\mathbf{p}}|\mathbf{x},\mathbf{p},t \rangle=\mathbf{j}_p (\mathbf{x},\mathbf{p},t)/\mathcal{P}(\mathbf{x},\mathbf{p},t)\), required to project the time derivatives defined along a stochastic trajectory into the ensemble space where the PDF \(\mathcal{P}(\mathbf{x},\mathbf{p},t)\) is defined. The second step is taken by averaging over the phase space \(\Omega :=\{(\mathbf{x},\mathbf{p})\}\). As a result, the two-step average yields, for some general function \(g\) [1]: \[\langle\!\langle g(\mathbf{x},\mathbf{p}) \mathbf{\dot{p}} \rangle\!\rangle=\int d\mathbf{x}d\mathbf{p}\, g(\mathbf{x},\mathbf{p})\mathbf{j}_p,\] where \(\mathbf{j}_p\) is the probability current along the momentum dimension defined in 4 . See Appendix 6 for a proof of this relation. In the following derivations, we take into account the natural boundary conditions associated with the PDF [25]: \[\label{eq:natural95BC} \mathcal{P}\big|_{\partial\Omega}= {\boldsymbol{j}}\big|_{\partial\Omega} = 0.\tag{16}\] The average EPR is then: \[\begin{align} \label{eq:dot95S95tot} \dot{S}_{tot} (t)= \dot{\mathcal{S}}(t) + \dot{S}_m (t), \end{align}\tag{17}\] where we define \[\begin{align} &\dot{S}_{tot} \coloneq\langle\!\langle\dot{s}_{tot} \rangle\!\rangle = \int d\mathbf{z}\, \frac{j^\mu_d}{\mathcal{P}} D_{\mu\nu}j^\nu_d, \tag{18} \\ &\dot{S}_m \coloneq \langle\!\langle\dot{s}_m \rangle\!\rangle = \int d\mathbf{z}\, j^\mu_d D_{\mu\nu} (F^\nu_- -g^\nu_k \partial_{p_\alpha}g^\alpha_k)-\langle \partial_{p_\mu}F^\mu_+ \rangle,\tag{19} \end{align}\] again consistent with [26] for additive noise. Notably, only the third term in Eq. 15 survives the averaging process. Provided the covariance matrix \(D^{\nu\mu}\) is positive definite, the average total EPR yields a quadratic form and hence verifies \(\dot{S}_{tot} \geq 0\), with equality being restricted to equilibrium where dissipative currents vanish \(( {\boldsymbol{j}}_d= \underline{0})\).

Moreover, \(\dot{\mathcal{S}}= \langle\!\langle\dot{s} \rangle\!\rangle\), which can be shown by directly differentiating Shannon entropy (see [27] for total time derivatives of integrals): \[\begin{align} \dot{\mathcal{S}}&=& - \int d\mathbf{z}\, \ln \mathcal{P}\, \frac{\partial\mathcal{P}}{\partial t} \nonumber\\ &=& \int d\mathbf{z}\, \ln \mathcal{P}\, [\partial_{x_\mu}\left(\frac{p^\mu}{m}\mathcal{P}\right) + \partial_{p_\mu} (j^\mu_r + j^\mu_d) ]. \end{align}\] Applying integration by parts, the term \(\langle \partial_{\mathbf{x}}\cdot(\mathbf{p}/m) \rangle=0\) vanishes due to independence of phase space variables. The reversible current contribution gives \[\begin{align} \int d\mathbf{z}\, \ln \mathcal{P}\, \partial_{\mathbf{p}} \cdot {\boldsymbol{j}}_r&=\int d\mathbf{z}\, \mathcal{P}\partial_{p_\mu}\frac{j^\mu_r}{\mathcal{P}} = \langle \partial_{p_\mu}F^\mu_+ \rangle. \end{align}\] For the dissipative current contribution: \[\begin{align} &\int d\mathbf{z}\, \ln \mathcal{P}\, \partial_{\mathbf{p}}\cdot {\boldsymbol{j}}_d= -\int d\mathbf{z}\, \frac{j^\mu_d}{\mathcal{P}}\partial_{p_\mu} \mathcal{P}\\ & = -\int d\mathbf{z}\, \frac{j^\mu_d}{\mathcal{P}}D_{\mu\nu}\left( (F^\nu_- -g^\nu_k \partial_{p_\alpha}g^\alpha_k)\mathcal{P}-j^\nu_d \right)\\ &= \int d\mathbf{z}\, \frac{j^\mu_d}{\mathcal{P}} D_{\mu\nu}j^\nu_d -\int d\mathbf{z}\, j^\mu_d D_{\mu\nu} (F^\nu_- -g^\nu_k \partial_{p_\alpha}g^\alpha_k). \end{align}\] Gathering the results, the average system EPR gives: \[\label{eq:dS95multid} \begin{align} \dot{\mathcal{S}}&=\int d\mathbf{z}\, \frac{j^\mu_d}{\mathcal{P}} D_{\mu\nu}j^\nu_d \\ &\phantom{\dot{\mathcal{S}}=}-\int d\mathbf{z}\, j^\mu_d D_{\mu\nu} (F^\nu_- -g^\nu_k \partial_{p_\alpha}g^\alpha_k) +\langle \partial_{p_\mu}F^\mu_+ \rangle\\ &=\langle\!\langle\dot{s}_{tot} \rangle\!\rangle - \langle\!\langle\dot{s}_m \rangle\!\rangle = \langle\!\langle\dot{s} \rangle\!\rangle. \end{align}\tag{20}\] Hence, while the Shannon entropy is a state function of the PDF, the medium entropy—and consequently, the total entropy—is inherently path-dependent. As such, there exists no state function \(S_m\) whose time derivative yields the average medium EP, i.e., \(\nexists S_m | \frac{d}{dt}S_m = \dot{S}_m\). The total and medium entropy are therefore defined explicitly through the trajectory-dependent expressions given by Eqs. 18 and 19 .

2.3 Stochastic energetics↩︎

Stochastic energetics deals with the stochastic extension of the first law and its connections with the second law obtained via ST. We illustrate it with the system Eq. 2 , but we distinguish three different kind of forces, \(\mathbf{F} = -\nabla_{\mathbf{x}}V(\mathbf{x},\lambda)-\underline{\underline{\boldsymbol{\gamma}}} (\mathbf{x},\mathbf{p}) \mathbf{p}/m+\mathbf{f}_+(\mathbf{x},\mathbf{p})\), where \(V(\mathbf{x},\lambda)\) gives the potential energy, \(\underline{\underline{\boldsymbol{\gamma}}}\) the damping coefficient and \(\mathbf{f}_+\) non-conservative even forces, as could be an external manipulation by an agent [26] or feedback cooling [28]. By analogy with the previous section, we identify the reversible and irreversible contributions as \(\mathbf{F}_+ = -\nabla_\mathbf{x} V(\mathbf{x},\lambda) +\mathbf{f}_+\) and \(\mathbf{F}_- =-\underline{\underline{\boldsymbol{\gamma}}}\cdot\mathbf{p}/m\). Following Sekimoto’s approach [8], we can derive the first law from the total energy of the system \(E\): \[\label{eq:1law} \begin{align} dE &= d\left(\frac{|\mathbf{p}|^2}{2m}+V(\mathbf{x},\lambda)\right)=\left(\dot{\mathbf{p}}+\nabla V\right) \cdot\frac{\mathbf{p}}{m}dt+\frac{\partial V}{\partial \lambda}d\lambda \\ & = ( -\underline{\underline{\boldsymbol{\gamma}}} (\mathbf{x},\mathbf{p})\cdot\mathbf{p}/m +\underline{\underline{\mathbf{g}}}(\mathbf{x},\mathbf{p})\circ\boldsymbol{\zeta}) \cdot \frac{\mathbf{p}}{m}dt\\ &\phantom{=} +\mathbf{f}_+ \cdot \frac{\mathbf{p}}{m}dt +\frac{\partial V}{\partial \lambda} d\lambda. \end{align}\tag{21}\] Then, heat and work can be identified as (in units \(k_B= 1\)) [1], [29]: \[\begin{eqnarray} &&\dot{q} \coloneq ( -\underline{\underline{\boldsymbol{\gamma}}} (\mathbf{x},\mathbf{p})\cdot\mathbf{p}/m +\underline{\underline{\mathbf{g}}}(\mathbf{x},\mathbf{p})\circ\boldsymbol{\zeta})\cdot \frac{\mathbf{p}}{m}, \tag{22}\\ &&\dot{w} \coloneq \mathbf{f}_+ \cdot \frac{\mathbf{p}}{m} +\frac{\partial V}{\partial \lambda} \dot{\lambda}. \tag{23} \end{eqnarray}\] This identification stems from the fact that heat involves irreversible forces, that is, forces that contribute to entropy production, meanwhile work involves reversible forces that do not explicitly contribute to the total entropy production [30].

In the particular scenario where \(\underline{\underline{\boldsymbol{\gamma}}}(\mathbf{x},\mathbf{p})=\underline{\underline{\boldsymbol{\gamma}}}\) is a constant, the Einstein relation \(T\underline{\underline{\boldsymbol{\gamma}}}=\underline{\underline{\mathbf{D}}}\), where \(\underline{\underline{\mathbf{D}}}\) is the covariance matrix and \(T\) is the absolute temperature, implies for Eq. 14 : \[\begin{align} \dot{s}_m &=& -\frac{1}{T}( -\underline{\underline{\boldsymbol{\gamma}}} \cdot\mathbf{p}/m +\underline{\underline{\mathbf{g}}}\circ\boldsymbol{\zeta})\cdot \frac{\mathbf{p}}{m} -\nabla_{\mathbf{p}}\cdot \mathbf{f}_+ \nonumber \\ &=& -\frac{\dot{q}}{T}- \dot{\mathcal{I}}, \end{align}\] recovering the Clausius relation in agreement with the first law Eq. 21 , as well as an additional contribution which can be identified as the phase space volume dilation rate difference between the forward and backward trajectories: \[\label{eq:entropy95pump} \dot{\mathcal{I}}=\frac{1}{2}\left[\nabla_{\mathbf{p}}\cdot \mathbf{F}(z)-\nabla_{\hat{\mathbf{p}}}\cdot \mathbf{\hat{F}}(z)\right].\tag{24}\] This term, referred to as entropy pumping, is attributed to the velocity-dependent control force due to an external agent [31]. It is then uncovered to be a lower bound on other information measures, such as information flow, accounting for the information flux due to external manipulation [30]. In Ref. [30] the information-theoretic interpretation of entropy pumping was established as a minimal information requirement for the agent to perform the momentum-dependent force \(\mathbf{f}_+\), an interpretation that has been further adopted [26], [28], [32]–[34]. Hence, the EP in the medium includes both dissipative and information terms, as it is well established in the literature [31], [35]–[40]: \[\begin{align} \label{eq:Gen95clausius} \Delta s_m = -\int\frac{dq}{T}-\Delta\mathcal{I}, \end{align}\tag{25}\] and we refer to Eq. 25 as the generalized Clausius relation, see Ref. [26], [30]. As a remark, we define the medium as the collection of all unobserved DoFs, which include those associated with the system (the heat bath) and those external to it, such as the agent performing work. Accordingly, the heat bath is a subset of the medium.

2.4 Thermodynamic consistency↩︎

The average total EP in ST 7 has been shown to always satisfy the second law by definition, due to the properties of the KL divergence. The second law can also be derived by applying Jensen’s inequality to the IFT, which, however, holds for any Markovian stochastic dynamics that do not necessarily represent a physical process [1]. Thus, we can conclude that this second law is not a true thermodynamical law, since it is not influenced by the thermodynamic consistency of the model. This has been previously pointed out in the literature by identifying the total EP as an information measure, termed informatic entropy production rate (IEPR), which, although providing a measure of irreversibility, it cannot be associated with thermodynamic quantities such as heat dissipation or the first law [5], [29].

As a consequence, when applying ST to, for instance, active matter systems, terms without a clear thermodynamic interpretation appear in the EP, identified as unconventional EP [26], [34], [41], [42]. This is because such models, of phenomenological origin, do not explicitly account for the driving mechanisms, therefore lacking thermodynamic consistency [5], [43]. Moreover, the lack of thermodynamic consistency can also be found in models that assume the Einstein relation without taking into consideration multiplicative noise or the associated noise-induced drifts, as it was pointed out in Refs. [5], [14].

Thermodynamic consistency finds different definitions in the literature. In this work, we combine its different meanings and refer to a thermodynamically consistent model as that which verifies:

1. Non-decreasing average total EPR \((\dot{S}_{tot} \geq 0)\) [7], [39], [44].

2. The generalized Clausius relation Eq. 25 such that the thermodynamic entropy is well-defined [5], [6], [26], [30], [41].

3. For isolated systems, the FP satisfies an equilibrium steady state [3], [4], [14], [45].

In the following section, we present a variational formulation of ST which systematically incorporates this notion of thermodynamically consistent modeling.

3 Variational principle of ST↩︎

3.1 Thermodynamic phase space↩︎

Following the approach described in Ref. [7], we derive the dynamical equations (ODEs) for both configurational and thermal variables by extending the phase space \(\Omega\) to include all relevant degrees of freedom. The fundamental thermal variable introduced in this framework is the stochastic thermodynamic entropy \(s\). For example, in a thermomechanical system, the extended phase space becomes \(\Omega:=\{(\mathbf{x},\mathbf{p},s)\}\).

Hence, \(s\) is an independent state variable that accounts for the hidden DoFs of the system and is no longer given as a function of \(\mathcal{P}\) as in 11 , nor is it given as a function of state [24]. This aligns with the fundamental path-dependence of the medium entropy highlighted in Sec. 2.2. Furthermore, this treatment has proven particularly effective in the deterministic (macroscopic) setting for describing nonequilibrium thermodynamics [7], where defining an entropy function is especially challenging [46].

In stochastic descriptions—such as Langevin dynamics—the assumption of time-scale separation renders the heat bath concept intrinsic to the formulation. The thermodynamic entropy \(s\) thereby characterizes the macrostate of the bath, whose dynamics are governed by heat exchange with the system, \(T\,ds = -dq\). This interpretation is consistent with the microscopic dynamical entropy framework introduced in Ref. [47], in which the bath contribution to the total entropy can be inferred from the system’s observable DoFs via the Clausius relation.

In contrast to the standard framework of stochastic thermodynamics, the PDF is now defined over the complete thermodynamic phase space \(\Omega\), acquiring the form \(\mathcal{P}(\mathbf{x},\mathbf{p},s,t)\) following the example of a thermomechanical system. The associated probability measure becomes \(\mathcal{P}(\mathbf{x},\mathbf{p},s,t)d\mathbf{x}d\mathbf{p}ds\). Accordingly, the stochastic system entropy and the Shannon entropy are defined as follows: \[\label{eq:shannon95VP} \begin{align} & s_{sys}(t)= - \ln \mathcal{P}(\mathbf{x}(t),\mathbf{p}(t),s(t),t),\\ &\mathcal{S}(t) = \langle s_{sys}(t) \rangle = -\int_\Omega d\Omega\, \mathcal{P}(\mathbf{x},\mathbf{p},s,t)\ln \mathcal{P}(\mathbf{x},\mathbf{p},s,t), \end{align}\tag{26}\] where \(d\Omega =d\mathbf{x}d\mathbf{p}ds\). The FP equation is denoted by \(\partial_t\mathcal{P}= -\nabla_\Omega\cdot\mathbf{j}\), where \(\nabla_\Omega\) denotes the gradient operator with respect to the thermodynamic phase space variables. In this particular example, \(\nabla_\Omega=(\partial_\mathbf{x},\partial_\mathbf{p},\partial_s)\). Note that, without assuming an infinite heat bath, the distribution \(\mathcal{P}\) encodes not only the statistics of the observable DoFs but also those of the bath macrostate, through its dependence on the thermodynamic entropy \(s\). In subsequent sections it is shown that, under the infinite-bath assumption, the Maxwell–Boltzmann distribution is recovered at equilibrium.

More generally, we denote the thermal variables by \(\varphi_\alpha\), which encompass both observable macroscopic quantities—such as the mass \(N\) of a chemical species—as well as thermodynamic entropy \(s\). These variables span the thermodynamic phase space, which we define as \(\Omega := \{(\mathbf{x}, \mathbf{p}, \varphi_\alpha)\} = \{(\boldsymbol{\omega}, s)\}\), where \(\boldsymbol{\omega}\) collectively denotes the observable DoFs. Associated with each phase space variable is a conjugate variable, referred to as the (generalized) thermodynamic displacement \(\Lambda^\alpha\). These displacements play a central role in our formulation, as they enable the derivation of the general form of the dynamical equations through a variational principle.

An additional central quantity in this formulation is the medium entropy, denoted by \(\Sigma\). This quantity plays a role analogous to the medium entropy \(s_m\) in conventional ST, in that it accounts for the entropy production associated with the medium—as defined in Sec. 2.3. Introducing \(\Sigma\) allows us to distinguish between entropy changes that contribute directly to the system’s thermodynamic entropy \(s\) and those arising from interactions with the external environment. The dynamics associated with \(\Sigma\) are determined by the nonlinear nonholonomic constraints, formulated within the framework of the Lagrange-d’Alembert principle.

3.2 Lagrange-d’Alembert principle↩︎

In Ref. [7] the variational approach to nonequilibrium thermodynamics is formulated by constructing a generalization of the Lagrange-d’Alembert principle of nonholonomic mechanics [48], where the entropy production of the system, written as the sum of the contribution of each of the irreversible processes, is incorporated into a nonlinear nonholonomic constraint. We now extend this formulation to the stochastic regime.

Recall that Hamilton’s Principle defines the dynamics of a physical system by the variational principle: \[\begin{eqnarray} &&\delta \mathcal{A}[\mathbf{x}] = \delta \int_{t_1}^{t_2}dt L (\mathbf{x}, \mathbf{\dot{x}})=0,\\ &&\delta\mathbf{x}(t_1) = \delta\mathbf{x}(t_2)=0, \end{eqnarray}\] where \(\mathbf{x}\in V\) stands for generalized coordinates,² \(\mathbf{\dot{x}} = d\mathbf{x}/dt \in V\) for generalized velocities, and \(L(\mathbf{x},\mathbf{\dot{x}})=K(\mathbf{x},\mathbf{\dot{x}})-U(\mathbf{x})\) is the Lagrangian function, which contains the physical information in terms of the kinetic energy \(K:V\times V\rightarrow \mathbb{R}\) and the potential energy \(U:V\rightarrow \mathbb{R}\). Here, we consider the configuration space \(V\) as a vector space.

In constructing a Lagrange-d’Alembert principle for ST we first extend the Lagrangian to the complete thermodynamic phase space by including the thermal contributions given by the internal energy \(U(\mathbf{x},\varphi_\alpha)\): \[\label{eq:Lagrangian} L(\mathbf{x}, \mathbf{\dot{x}}, \varphi_\alpha)=K(\mathbf{x},\mathbf{\dot{x}})-U(\mathbf{x},\varphi_\alpha).\tag{27}\] This internal energy encodes the properties of the heat bath, such as its temperature \(T(\mathbf{x},\varphi_\alpha)\coloneq \partial U/ \partial s\) or the exerted force on the system \(\mathbf{f}(\mathbf{x},\varphi_\alpha)\coloneq -\partial U/\partial\mathbf{x}\). These relations correspond the total differential of the internal energy, expressed in terms of the entropy and the observable variables \(\boldsymbol{\omega}\) as \(dU=TdS+ \nabla_{\boldsymbol{\omega}}U \cdot d\boldsymbol{\omega}\) [49]. Notably, this definition of temperature arises explicitly from the internal energy and will appear consistently in the fluctuation-dissipation relation. This differs from earlier instances where \(T\) appeared merely as a positive parameter enforcing \(T\underline{\underline{\boldsymbol{\gamma}}}=\underline{\underline{\mathbf{D}}}\), and was only a posteriori interpreted as the bath temperature.

As a technical step, we derive the implicit form of the Euler-Lagrange equations—i.e., the equations for \((\dot{\mathbf{x}}, \dot{\mathbf{p}}, \dot{\varphi}_\alpha)\)—using the Hamilton-Pontryagin principle [50]. This formulation casts the stochastic dynamics as a first-order system on the thermodynamic phase space, ensuring that all time derivatives are well-defined as stochastic differentials and compatible with the chosen integration convention. In the case of the extended Lagrangian 27 , the only second-order kinematic constraint is \(\dot{\mathbf{x}} = \mathbf{v}\), with fixed endpoints for \(\mathbf{x}(t)\). Accordingly, we promote \(\mathbf{p}\) to a Lagrange multiplier enforcing this constraint, by augmenting the action with the term \(\int dt\, \mathbf{p} \cdot (\dot{\mathbf{x}} - \mathbf{v})\), consistently with the variational formulation of implicit Lagrangian systems developed in [50].

However, to consistently recover the reversible dynamics, the variational principle must yield conservation laws for the thermal variables in the absence of irreversible processes, i.e., \(\dot{\varphi}_\alpha = 0\). This is accomplished by coupling each thermal variable \(\varphi_\alpha\) to its conjugate displacement variable \(\Lambda^\alpha\) within the action functional, through a term of the form \(\int dt\, \varphi_\alpha \dot{\Lambda}^\alpha\) [7]. In the case of entropy, where the time evolution of the medium entropy \(\Sigma\) is governed by the nonholonomic constraints, this coupling takes the form \(\int dt\, (s-\Sigma) \dot{\Gamma}\), where \(\Gamma\) is the conjugate variable to the entropy \(s\). The minus sign in this expression ensures that, in the absence of external fluxes, the variational principle gives \(\dot{\Sigma} = \dot{s}\),—as expected for an adiabatically closed system—where the thermodynamic EP accounts for the entire medium entropy.

Thus, the complete set of dynamical equations are obtained by taking variations of the (Hamilton-Pontryagin modified) action given by the extended Lagrangian, e.g. Eq. 27 , \[\begin{align} \delta \int_{t_1}^{t_2} dt \bigg( L (\mathbf{x},\mathbf{v},\varphi_\alpha,s) + \mathbf{p} \cdot ( \mathbf{\dot{x}} - \mathbf{v}) \nonumber\\ + \varphi_\alpha\dot{\Lambda}^\alpha+(s-\Sigma)\dot{\Gamma} \bigg) =0, \end{align}\] subject to the nonholonomic constraints (to be defined). External manipulations can be accounted for by including a potential \(\tilde{U}(\mathbf{x},\lambda)\) for a time-dependent driving, or an external non-conservative force \(\mathbf{f}_+\) via the principle of virtual work.

To motivate the form of the nonholonomic constraints, we draw an analogy with standard ST, in the same spirit as the geometric formulation of nonequilibrium thermodynamics in Ref. [7]. In our framework, the thermodynamic power associated with an irreversible process mirrors the structure of the heat dissipated by a mechanical force, as expressed in Eq. 22 , or equivalently as \(dq = (\mathbf{f}_- + \boldsymbol{\xi}) \cdot d\mathbf{x}\). That is, it takes the form of a generalized friction force—comprising both deterministic and stochastic contributions—contracted with a mechanical displacement.

Translating this structure to the thermodynamic setting, we interpret the thermodynamic fluxes \(J_\alpha\) as generalized friction forces, and the differentials of the thermodynamic displacements \(d\Lambda^\alpha\) as generalized displacements. This yields an analogous expression for dissipated heat: \(dq = J_\alpha\, d\Lambda^\alpha\). By introducing a thermodynamic phase space \(\Omega\) equipped with a set of thermal variables, and its dual \(\Omega^*\) comprising their corresponding conjugate displacements, this ansatz ensures that the structure of dissipation in general irreversible processes inherits the same geometric form as that of mechanical friction.

With this preamble, we have shown how irreversible processes can be incorporated within the same geometric structure, with entropy production providing the underlying organizing principle. We then embed this structure into a variational formulation by employing the Lagrange-d’Alembert principle [7], [48], in which the critical curve condition is subject to two constraints: a kinematic constraint on the solution curve and a variational constraint on the variations to be considered when computing the criticality condition. By construction, these constraints encode the information regarding the entropy production. Thus, we define the kinematic and variational nonlinear nonholonomic constraints, respectively, as: \[\label{eq:nonhol95constraints} \begin{align} & \frac{\partial L}{\partial s}\dot{\Sigma} = (J_\alpha+\xi_\alpha)\cdot \dot{\Lambda}^\alpha+ (J_\beta+\xi_\beta)\cdot (\dot{\Lambda}^\beta-X^\beta_{ext}) + \mathcal{I}_\beta\dot{\Gamma}^\beta, \\ & \frac{\partial L}{\partial s}\delta\Sigma=(J_\alpha+\xi_\alpha)\cdot \delta\Lambda^\alpha+ (J_\beta+\xi_\beta)\cdot\delta\Lambda^\beta+ \mathcal{I}_\beta\delta\Gamma^\beta, \end{align}\tag{28}\] where \(\alpha\) and \(\beta\) denote internal and external processes, respectively, with thermodynamic (irreversible) fluxes \(J_\alpha, J_\beta\) and thermodynamic affinities \(X^\alpha, X^\beta\), together with a thermodynamic affinity \(X^\beta_{ext}\) associated with the exterior. Then, the (generalized) thermodynamic displacements \(\Lambda^\alpha,\Lambda^\beta\) are such that \(\dot{\Lambda}^\alpha= X^\alpha\) and \(\dot{\Lambda}^\beta= X^\beta\) [7]. We recall \(\Sigma\) stands for the medium entropy.

In contrast to the deterministic (macroscopic) formulation, each dissipative flux \(J_\alpha,J_\beta\) must include its stochastic counterpart \(\xi_\alpha,\xi_\beta\) for complete thermodynamic consistency, according to the fluctuation-dissipation theorem (FDT) [3], [4], [45]. This arises as a consequence of coarse-graining: dissipative forces constitute a macroscopic representation of underlying microscopic processes whose DoFs remain unresolved. The lack of information on these DoFs then has to be explicitly accounted for by the noise in order to correctly quantify the associated entropy production.

An additional class of thermodynamic fluxes arising in the stochastic framework corresponds to information fluxes. In the context of Langevin dynamics, Sec. 2.3 showed that external manipulations can induce an entropy flow—termed entropy pumping. As this contribution can be regarded as an external entropy flux \(\mathcal{I}_\beta\)—see Eq. 25 —it can be incorporated into the constraint structure 28 by coupling it to the conjugate displacement of the thermodynamic entropy \((\Gamma)\), yielding additional terms \(\mathcal{I}_p \dot{\Gamma}\) and \(\mathcal{I}_p\, \delta\Gamma\) in the kinematic and variational constraints, respectively, where \[\label{eq:entropy95pump95VM} \mathcal{I}_p \coloneq \frac{1}{2}\left[\nabla_\Omega\cdot \boldsymbol{\mathfrak{F}}-\nabla_{\hat{\Omega}}\cdot \hat{\boldsymbol{\mathfrak{F}}}\right],\tag{29}\] and \(\boldsymbol{\mathfrak{F}}\) collectively denotes all the deterministic forces or fluxes acting on the system (recall \(\hat{\Box}\) denotes time-reversal operation). Equation 29 is just the extension of Eq. 24 to the thermodynamic phase space. Since entropy pumping arises from non-dissipative mechanisms, it is not constrained by the FDT and contributes purely through deterministic terms (i.e., \(\xi_{\beta} = 0\)). This is consistent with the standard framework of ST, as such contributions do not explicitly contribute to the total EP, see Eq. 18 . A concrete example will be provided in Sec. 3.5. See Fig. 1 for an illustration of the different thermodynamic fluxes.

Both type of constraints are systematically related, since the variational constraints can be derived form the kinematic constraints by replacing \(d\Lambda^\alpha\rightsquigarrow\delta\Lambda^\alpha\) and \((d\Lambda^\beta-X^\beta_{ext}dt) \rightsquigarrow \delta\Lambda^\beta\), where the exterior contribution vanishes since \(\delta t = 0\) by definition—variations in Hamilton’s principle are taken at fixed time [48], [51]. As a remark, this principle is not a critical curve condition for the action integral restricted to the space of a curve satisfying the constraints. However, it is emphasized this formulation recovers Hamilton’s Principle in the absence of irreversibility, analogously as in Ref. [7].

Up to this point, we have not specified any FDRs. The FDRs are obtained by treating the second law as an axiom, imposed on the dynamical equations derived from the variational formulation. Within this framework, the second law is stated as: \[\label{eq:2nd95law95VP} \dot{S}_{tot}(t) \coloneq \dot{\mathcal{S}}(t) + \langle\!\langle\dot{\Sigma} (t) \rangle\!\rangle \geq 0,\tag{30}\] where \(\mathcal{S}(t)\) stands for the Shannon entropy as given by Eq. 26 . This definition has the same form as that of ST, i.e., Eq. 17 , except that the inequality is only satisfied under the FDRs and not otherwise, contrary to the standard definition in ST. Moreover, this expression is defined in the complete thermodynamic phase space \(\Omega\).

In subsequent examples we show how, by construction, this is the only condition we need to impose in this formulation to achieve thermodynamic consistency as defined in Sec. 2.4, the rest of conditions arising as a consequence of the second law. Hence, we do not require a specific steady state condition nor the system being either closed or open, such that no thermodynamic equilibrium condition is enforced.

We analyze three different examples: a thermomechanical isolated system (3.3), an interconnected system of multiple reservoirs (3.4), and a thermomechanical open system (3.5).

3.3 Mechanical isolated system↩︎

First, we consider a mechanical isolated system of a particle in a thermal bath, with thermodynamic phase space \(\Omega :=\{(\mathbf{x}, \mathbf{p}, s)\}\) and Lagrangian \[L(\mathbf{x}, \mathbf{v}, s)=\frac{1}{2}m|\mathbf{v}|^2-U(\mathbf{x},s)\] such that \(\partial L/\partial s = -\partial U/\partial s = -T(\mathbf{x},s)\). Note the potential \(V(\mathbf{x})\) can be absorbed in \(U(\mathbf{x},s)\). The variational formulation then involves taking variations of the (Hamilton-Pontryagin modified) action, subject to Lagrange-d’Alembert nonholonomic constraints for the entropy \(\Sigma\): \[\label{eq:VP95mech95closed} \begin{align} &\delta \int_{t_1}^{t_2}dt\left( L (\mathbf{x},\mathbf{v},s) + \mathbf{p} \cdot ( \mathbf{\dot{x}} - \mathbf{v}) + \dot{\Gamma}(s-\Sigma)\right) =0, \\ &\begin{cases} \displaystyle\frac{\partial L}{\partial s}\dot{\Sigma}= (\mathbf{f_-}+ \underline{\underline{\mathbf{g}}}\circ\boldsymbol{\zeta}) \cdot\mathbf{\dot{x}}, \\ \displaystyle\frac{\partial L}{\partial s}\delta \Sigma= (\mathbf{f_-}+ \underline{\underline{\mathbf{g}}}\circ\boldsymbol{\zeta}) \cdot\delta \mathbf{x}, \end{cases} \end{align}\tag{31}\] where \(\mathbf{f_-} = -\underline{\underline{\boldsymbol{\gamma}}} (\mathbf{x},\mathbf{p},s) \mathbf{p}/m\) is a dissipative force with state-dependent damping coefficient, \(\underline{\underline{\mathbf{g}}}(\mathbf{x},\mathbf{p},s)\) is the noise strength and \(\boldsymbol{\zeta}(t)\) is gaussian white noise verifying the relations given by Eq. 3 . \(\Gamma\) denotes the (generalized) thermal displacement, i.e. \(\dot{\Gamma} = T\), an equality that follows from the variational principle.

For simplicity, we assume a 1-dimensional system. Then, the variational principle gives the equations of motion (EoMs): \[\begin{align}\label{eq:EoM95mech95closed} & \displaystyle p=m\dot{x},\\ & \displaystyle\dot{p}=-\partial_x U(x,s) -\gamma(x,p,s)\frac{p}{m}+ g(x,p,s)\circ\zeta,\\ & \displaystyle\dot{s} = \frac{-1}{T(x,s)}\left( -\gamma(x,p,s)\frac{p}{m}+g(x,p,s)\circ\zeta\right)\frac{p}{m}, \end{align}\tag{32}\] with \(\dot{s} = \dot{\Sigma}\), equality that holds in any closed system. See Appendix 7 for the detailed derivation. Taking the total time derivative of \(E = \frac{1}{2}m|\mathbf{v}|^2 + U(\mathbf{x},s)\) we find the first law \[\label{eq:1law95mech95closed} \dot{E} = (\dot{p} + \partial_x U)\frac{p}{m}+ T\dot{s} =0,\tag{33}\] hence energy is conserved and the system is closed. Moreover, the equation for \(\dot{s}\) has the same form as Eq. 22 , following the Clausius relation and confirming that \(s\) consistently represents the thermodynamic entropy, as previously established.

Denoting the dissipative and entropy probability currents, respectively, as \[\begin{eqnarray} j_d =&& \left(-\gamma\frac{p}{m}-g\partial_p g+\frac{p}{m}g\partial_s (g/T)\right)\mathcal{P}\nonumber\\ &&-g^2\partial_p \mathcal{P}+\frac{pg^2}{mT}\partial_s \mathcal{P}, \\ j_s =&& -\frac{p}{mT} { j}_d, \end{eqnarray}\] the system and thermodynamic average EPRs yield: \[\begin{align} &&\dot{\mathcal{S}}= \int d\Omega\,\left\{ \frac{ { j}_d^2}{\mathcal{P}g^2} \right.\nonumber\\ &&\phantom{\dot{\mathcal{S}}=} - \left. \frac{ { j}_d}{\;g^2} \left(-\gamma\frac{p}{m}-g\partial_p g+\frac{p}{m}g\partial_s (g/T)\right)\right\},\\ &&\langle\!\langle\dot{s} \rangle\!\rangle = -\int d\Omega\, \frac{p}{mT} { j}_d. \end{align}\] Enforcing the second law we then derive the FDR: \[\begin{align} &&\dot{S}_{tot}=\dot{\mathcal{S}}+\langle\!\langle\dot{s} \rangle\!\rangle = \int d\Omega\, \frac{ { j}_d^2}{\mathcal{P}D}\geq 0 \iff \\ &&-\gamma\frac{p}{m} = -\frac{D}{T}\frac{p}{m}+g\partial_p g - g\partial_s (g/T)\frac{p}{m}, \label{eq:FDR95mech95closed} \end{align}\tag{34}\] where we recall \(D=g^2\). The total EPR has the same form as that obtained from the standard ST approach in one dimension, see Eq. 18 , but it is instead defined over the thermodynamic phase space \(\Omega\). The corresponding multidimensional version of Eq. 34 , obtained by an analogous computation, is \[-\gamma^{\nu\mu}\frac{p_\mu}{m} = -\frac{D^{\nu\mu}}{T}\frac{p_\mu}{m}+g^\nu_k \partial_{p_\mu}g^\mu_k - g^\nu_k\partial_s (g^\mu_k/T)\frac{p_\mu}{m}.\]

Remarkably, the thermodynamic EP satisfies \[\begin{align} \frac{d}{dt}\langle s \rangle &= \int d\Omega\, s \partial_t \mathcal{P}= -\int d\Omega\, s \nabla_\Omega \cdot\mathbf{j} = \int d\Omega\, \mathbf{j}\cdot\nabla_\Omega s \\ &= \int d\Omega\, j_s=-\int d\Omega\, \frac{p}{mT} { j}_d=\langle\!\langle\dot{s} \rangle\!\rangle, \end{align}\] which implies that \(\dot{S}_{tot}\) is a real time derivative, as we can write \[\label{real95time95rate} \dot{S}_{tot}= \frac{d}{dt}(\mathcal{S}+ \langle s\rangle),\tag{35}\] in contrast with the standard ST setting, i.e. the results of Sec. 2.

Proceeding to analyze the FDR in Eq. 34 , we observe that it connects the dissipative forces with the noise amplitude, and also includes contributions from noise-induced drift terms. In the special case of constant temperature and constant damping coefficient, we recover the standard Einstein relation discussed in Sec. 2.3, namely \(T\gamma = D\).

When \(\gamma,T,g\) are independent of the entropy \(s\), our result reproduces the FDR derived in Ref. [52], where this relation was derived by setting the steady-state distribution to be the Boltzmann distribution. This reference also illustrates how momentum-dependent multiplicative noise can arise, for instance, in systems with Coulomb friction, which requires an extension of the FDR to ensure thermodynamic consistency. Additionally, in this scenario—or more generally, for \(\partial_p {\rm f}_+=0\)—where the observable DoFs decouple from the entropy balance in 32 , \(\dot{s}\) recovers the same expression as \(\dot{s}_m\) in Eq. 14 under the FDR constraint, providing it with a well-defined thermodynamic interpretation (Clausius relation).

Beyond recovering these known results, this framework introduces a novel approach for coupling mechanical and thermal DoFs [24]. This allows for taking into account finite baths and/or coupled system-bath dynamics, thus going beyond the weak coupling limit. The third contribution in Eq. 34 , the noise-induced drift due to \(s\) dependence, is essential to guarantee thermodynamic consistency in the presence of such nonlinear system-bath coupling.

Addressing the existence of an equilibrium steady state distribution, in order to complete the requirements of thermodynamic consistency provided in Sec 2.4, by imposing the FDR 34 we find that the functional \[\label{eq:SS95functional} \mathcal{P}(x,p,s)=f\left(E(x,p,s)\right)e^s\tag{36}\] satisfies the FP at steady state for any function of the total energy \(f(E)\). Computing the current components gives: \[\begin{align} &\partial_x j_x = \frac{p}{m}\partial_x \mathcal{P}= \frac{p}{m}e^s f'(E)\partial_x U ,\\ & \partial_p j_p = \partial_p \left(-\partial_x U\mathcal{P}-\frac{pg^2}{mT}\mathcal{P}+\frac{pg^2}{mT}\mathcal{P}\right)= -\partial_x j_x,\\ &\partial_s j_s = \partial_s \left(\left[\frac{pg}{mT}\right]^2\mathcal{P}-\left[\frac{pg}{mT}\right]^2\mathcal{P}\right)=0,\\ &\nabla_\Omega \cdot\mathbf{j} = \partial_x j_x - \partial_x j_x = 0, \end{align}\] which implies \(\partial_t \mathcal{P}=-\nabla_\Omega\cdot\mathbf{j}=0\). The distribution 36 is general and holds for coupled system-bath dynamics as well. Moreover, it favors the states that maximize the entropy \(s\), consistent with the maximum entropy principle in statistical mechanics, as well as being of the same form as the equilibrium distribution in fluctuation theory [49], [53]–[55], with an additional prefactor accounting for energy conservation. The correspondence between this equilibrium state and fluctuation theory arises since both assume gaussian fluctuations.

An explicit form of Eq. 36 can be derived by using the fact that the energy distribution \(\mathcal{P}(e)\) is time-independent, given energy is conserved, see Eq. 33 : \[\begin{align} \mathcal{P}_0(e)&=\mathcal{P}(e)=\int_{\Omega}d\Omega\, f(E(x,p,s))e^s\delta(e-E(x,p,s)) \\ &=f(e)\int_\Omega d\Omega\, e^s\delta(e-E(x,p,s)) = f(e)Z(e) \\ & \Rightarrow f(e) = \frac{ \mathcal{P}_0(e)}{Z(e)} . \end{align}\] The phase space integral for \(Z(e)\) can be solved applying the property [18]: \[\label{eq:composite95Dirac95delta} \delta(g(x))=\sum_i \frac{\delta(x-x_i)}{|g'(x_i)|} \quad \text{for}\quad g(x_i)=0,\tag{37}\] for which an explicit form of \(E(x,p,s)\) is required. In the simplest scenario where \(T=\partial_s E\) is a parameter, we assume \(U(x,s)=Ts\), such that the particle and bath energies are separable \(E = H_0 (x,p) + Ts\). Denoting \(s_0=(e-H_0(x,p))/T\): \[\begin{align} Z(e)&=\int_\Omega d\Omega\, e^s\delta(e-E(x,p,s)) =\int_\Omega dxdpds \, e^s\frac{\delta(s-s_0)}{|\partial_s E|}\\ &=\int_{X,P} dxdp \,\frac{1}{|T|}e^{\frac{1}{T}(e-H_0(x,p))}\\ & = \frac{1}{|T|}e^{e/T}\int_{X,P} dxdp \,e^{-H_0(x,p)/T} =\frac{1}{|T|}e^{e/T}\; Z_{MB}. \end{align}\] Hence, \[\begin{align} \label{eq:MB} \mathcal{P}(x,p,s) &=& \mathcal{P}_0(e)\frac{e^s}{Z(e)}\delta(e-E(x,p,s)) \nonumber\\ &\propto& \mathcal{P}_0(E(x,p,s))\frac{e^{s-e/T}}{Z_{MB}} \nonumber \\ &=& \mathcal{P}_0(E(x,p,s))\frac{e^{-H_0(x,p)/T}}{Z_{MB}}. \end{align}\tag{38}\] This corresponds precisely to the Maxwell-Boltzmann (MB) canonical distribution [49], [54], subject to energy conservation. We note that, for a uniform \(\mathcal{P}_0(e)\), the distribution becomes independent of \(s\). This situation corresponds to an infinite bath whose macroscopic state remains unchanged. Remarkably, in this formulation, the MB distribution emerges as a consequence of the FDR derived from the second law of thermodynamics, contrary to the conventional approach, where the FDR is obtained by assuming the MB distribution as the steady-state. Hence, the variational approach recovers the results of equilibrium statistical mechanics without being limited to it. The complete derivations supporting this section are presented in Appendix 7.

3.4 Interconnected system of multiple reservoirs↩︎

In this example, we consider multiple reservoirs interacting via stochastic fluxes. For simplicity, let us assume two reservoirs that exchange mass and heat, such that the whole system is closed. The entropy and mass variables for each reservoir are denoted \(\{(s_i ,N_i)\}\) for \((i=1,2)\), respectively. The Lagrangian corresponds to \(L=-U(s_1,s_2,N_1,N_2)\), which gives the temperature and chemical potential of each reservoir by \[\begin{align} T^k (s_1,s_2,N_1,N_2)= \frac{\partial U}{\partial s_k} , \\ \mu^k (s_1,s_2,N_1,N_2)= \frac{\partial U}{\partial N_k} \end{align}\] accordingly. The mass and entropy fluxes are given by \(\mathcal{J}_{12}=-\mathcal{J}_{21}\) and \(J_{12}=J_{21}\) respectively, with the stochastic counterparts given by the noises \(\zeta(t),\eta(t)\) with amplitudes \(\sigma=\sigma_{12}=-\sigma_{21}\) and \(\kappa=\kappa_{12}=\kappa_{21}\). These symmetries derive from total mass and energy conservation due to the system being closed. For the purpose of the variational formulation, it is convenient to define the flux \(J_{kl}\) for \(k=l\) as [51] \[\label{modeling95assumption} J_{kk}\coloneq -\sum_{l\neq k}J_{kl},\tag{39}\] so that \(\sum_k J_{kl}=0\). The associated constraints are then given by \[\frac{\partial L}{\partial s_k}\delta\Sigma_k=\sum_l(J_{kl}+\kappa_{kl}\circ\eta)\delta\Gamma^l +\sum_l (\mathcal{J}_{l\rightarrow k}+\sigma_{l\rightarrow k}\circ\zeta)\delta W^k .\] Thus, the variational formulation takes the form \[\begin{gather} \delta \int_{t_1}^{t_2} dt \left(L(s_1,s_2,N_1,N_2) + \dot{W}^1 N_1 + \dot{W}^2 N_2 \right.\\ \hfill \left. +\dot{\Gamma}^1(s_1-\Sigma_1)+\dot{\Gamma}^2(s_2-\Sigma_2) \right)=0,\\ \begin{cases} \displaystyle\frac{\partial L}{\partial s_1}\dot{\Sigma}_1=(J_{12}+\kappa_{12}\circ\eta)(\dot{\Gamma}^2-\dot{\Gamma}^1) \\ \hfill -(\mathcal{J}_{12}+\sigma_{12}\circ\zeta)\dot{W}^1 ,\\ \displaystyle \frac{\partial L}{\partial s_2}\dot{\Sigma}_2=(J_{12}+\kappa_{12}\circ\eta)(\dot{\Gamma}^1-\dot{\Gamma}^2) \\ \hfill +(\mathcal{J}_{12}+\sigma_{12}\circ\zeta)\dot{W}^2,\\ \displaystyle\frac{\partial L}{\partial s_1}\delta\Sigma_1=(J_{12}+\kappa_{12}\circ\eta)(\delta\Gamma^2-\delta\Gamma^1)\\ \hfill -(\mathcal{J}_{12}+\sigma_{12}\circ\zeta)\delta W^1, \\ \displaystyle \frac{\partial L}{\partial s_2}\delta\Sigma_2=(J_{12}+\kappa_{12}\circ\eta)(\delta\Gamma^1-\delta\Gamma^2)\\ \hfill +(\mathcal{J}_{12}+\sigma_{12}\circ\zeta)\delta W^2, \end{cases} \end{gather}\] where the variables \(\Gamma^i,W^i\) are conjugate variables with \(s_i\) (or \(\Sigma_i\)) and \(N_i\), respectively. Taking variations then gives the following system: \[\label{eq:EoM95interc} \begin{align} & \dot{s}_k = \dot{\Sigma}_k ,\\ & \dot{N}_1 = -\mathcal{J}_{12}-\sigma_{12}\circ\zeta ,\\ & \dot{N}_2 = \mathcal{J}_{12}+\sigma_{12}\circ\zeta,\\ & T^1\dot{s}_1=(J_{12}+\kappa_{12}\circ\eta)(T^1-T^2)+(\mathcal{J}_{12}+\sigma_{12}\circ\zeta)\mu^1 ,\\ & T^2 \dot{s}_2=(J_{12}+\kappa_{12}\circ\eta)(T^2-T^1)-(\mathcal{J}_{12}+\sigma_{12}\circ\zeta)\mu^2, \end{align}\tag{40}\] where fluxes and noises are allowed to be state-dependent, i.e. \(\mathcal{J}_{12}=\mathcal{J}_{12}(s_1,s_2,N_1,N_2),\,\sigma_{12}=\sigma_{12}(s_1,s_2,N_1,N_2)\), and similarly for \(J_{12},\kappa_{12}\). Moreover, for the sake of generality, we further allow the noises to be cross-correlated. Denoting \(\boldsymbol{\xi}=(\zeta,\eta)\): \[\langle \boldsymbol{\xi}(t):\boldsymbol{\xi}(t') \rangle = 2 \begin{bmatrix} 1 & C\\ C & 1 \end{bmatrix}\delta(t-t'),\] where \(C\) denotes the cross-correlation coefficient. It will be shown how cross-correlation implies a common hidden physical mechanism behind the heat and mass fluxes, giving rise to cross-effects (see Fig. 2).

The first law in this system gives, applying Eqs. 40 : \[\dot{E} = \mu^k \dot{N}_k + T^k \dot{s}_k = 0,\] such that energy is conserved and the system is closed.

To analyze the EPR we define the following currents, according to the FP equation \((\partial_t\mathcal{P}= -\nabla_\Omega\cdot\boldsymbol{j})\) for cross-correlated gaussian white noise [56]: \[\label{eq:currents95inter} \begin{align} &j_{N_1} = -\tilde{\mathcal{J}}_{12}\mathcal{P}- \sigma^2 (\partial_{N_1}\mathcal{P}-\partial_{N_2}\mathcal{P}) \nonumber\\ &\phantom{j_{N_1} =} -\sigma^2 \left(\frac{\mu^2}{T^2}\partial_{s_2}\mathcal{P}-\frac{\mu^1}{T^1}\partial_{s_1}\mathcal{P}\right) \nonumber\\ &\phantom{j_{N_1} =} +C\sigma\kappa\left(\frac{T^1-T^2}{T^1}\partial_{s_1}\mathcal{P}+\frac{T^2-T^1}{T^2}\partial_{s_2}\mathcal{P}\right) ,\\ &j_{N_2} = -j_{N_1},\\ &j_{s_1} = -\frac{\mu^1}{T^1}j_{N_1}+\frac{T^1-T^2}{T^1}j_{th},\\ &j_{s_2} = \frac{\mu^2}{T^2}j_{N_1}+\frac{T^2-T^1}{T^2}j_{th}. \end{align}\tag{41}\] For clarity, we have denoted the thermal contribution as \[\begin{align}\label{eq:j95th95inter} &j_{th} = \tilde{J}_{12}\mathcal{P}-\kappa^2(1-C^2)\left(\frac{T^1-T^2}{T^1}\partial_{s_1} \mathcal{P}+\frac{T^2-T^1}{T^2}\partial_{s_2}\mathcal{P}\right)\\ & \phantom{j_{th} = } -\frac{C\kappa}{\sigma}\left(\tilde{\mathcal{J}}_{12}\mathcal{P}+ j_{N_1}\right), \end{align}\tag{42}\] and absorbed the noise-induced drifts in \[\begin{align} \tilde{\mathcal{J}}_{12} = \mathcal{J}_{12}+ \nu_m, \\ \tilde{J}_{12} = J_{12} + \nu_{th}, \end{align}\] where the explicit expressions for \(\nu_m ,\nu_{th}\) can be found in Appendix 8. The average system and thermodynamic EPRs then yield, respectively,

Note that Eq. 35 is again satisfied since \(\sum_k\langle\!\langle\dot{s}_k \rangle\!\rangle= \frac{d}{dt} \sum_k\langle s_k \rangle\) (see Appendix 8). Adding both contributions and enforcing the second law gives the FDRs: \[\label{eq:FDR95inter} \dot{S}_{tot}\geq 0 \iff \begin{bmatrix} \tilde{\mathcal{J}}_{12} \\ \tilde{J}_{12} \end{bmatrix} = \mathbb{L} \begin{bmatrix} \displaystyle \frac{\mu^1}{T^1} - \frac{\mu^2}{T^2} \\ \displaystyle (T^1-T^2)\left(\frac{1}{T^1}-\frac{1}{T^2}\right) \end{bmatrix},\tag{43}\] with \[\label{eq:L95inter} \mathbb{L} = \begin{bmatrix}\sigma^2 & C\kappa\sigma\\ C\kappa\sigma& \kappa^2 \end{bmatrix},\; \mathbb{L}^{-1} = \frac{1}{(1-C^2)}\begin{bmatrix} \displaystyle\frac{1}{\sigma^2} & \displaystyle-\frac{C}{\kappa\sigma} \\ \displaystyle-\frac{C}{\kappa\sigma} & \displaystyle\frac{1}{\kappa^2} \end{bmatrix},\tag{44}\] being symmetric positive definite (PD) matrices. The average total EPR can then be written explicitly as a quadratic form: \[\dot{S}_{tot} = \int \frac{d\Omega}{\mathcal{P}} \begin{bmatrix} j_{N_2} & j_{th} \end{bmatrix} \mathbb{L}^{-1} \begin{bmatrix} j_{N_2} \\ j_{th} \end{bmatrix}\geq 0,\] with equality only at equilibrium, where dissipative currents vanish.

The FDR 43 shows that, due to cross-correlations, each flux is influenced by both temperature and chemical potential gradients, thereby exhibiting cross-effects. The relationship between the fluxes and the gradients is governed by the symmetric PD matrix \(\mathbb{L}\), which embodies the Onsager reciprocal relations. This structure indicates that, in the macroscopic limit, the model inherently recovers the phenomenological relations of nonequilibrium thermodynamics [51], [57], mediated by the Onsager matrix. Namely, the FDR 43 has the same form as its macroscopic counterpart found in [51] from \[\label{I95macroscopic} I:= \mathcal{J}_{kl}\left(\frac{\mu^k}{T^k} - \frac{\mu^l}{T^l}\right) + J_{kl}(T^l-T^k)\left(\frac{1}{T^l}-\frac{1}{T^k}\right)\geq 0,\tag{45}\] where \(I\) denotes the macroscopic internal EP. In addition, the fluxes include noise-induced drifts, arising from multiplicative noise and/or state-dependent correlations, which are necessary for thermodynamic consistency. These features are particularly relevant in modeling active field theories, as discussed in Refs. [5], [6]. However, the special case \(C=\pm 1\) leads to a degenerate situation in which the two noise terms become identical, introducing a redundancy and rendering the model ill-defined.

Since the system is closed, the functional \[\label{eq:SS95inter} \mathcal{P}(s_1,s_2,N_1,N_2)=f(E(s_1,s_2,N_1,N_2))e^{\sum_k s_k}\tag{46}\] satisfies the FP at steady state, analogously as in Sec. 3.3, see Eq. 36 . This can be shown by enforcing the FDRs 43 in the currents given by Eqs. 41 , together with the relations \[\begin{eqnarray} &&\partial_{N_i}\mathcal{P}= e^{\sum_k s_k}f'(E)\mu^i,\\ &&\partial_{s_i}\mathcal{P}= e^{\sum_k s_k}f'(E)T^i + \mathcal{P}. \end{eqnarray}\] Substituting in the probability currents yields: \[j_{N_i}=0,\;j_{th}=0 \Rightarrow \partial_t \mathcal{P}= -\nabla_\Omega \cdot \boldsymbol{j} = 0.\] Again, the system displays an equilibrium steady state under the FDRs derived through the second law, as we require for the thermodynamic consistency of closed systems. The functional 46 can be further simplified for specific forms of the energy \(E(s_1,s_2,N_1,N_2)\) by applying the composite rule for Dirac delta functions given by Eq. 37 .

For example, assuming that the reservoirs have equilibrated at a temperature \(T=T^1=T^2\) and the grand canonical free energy \(F(N_1,N_2)\) is separable, the energy function takes the form \(E(s_1,s_2,N_1,N_2)=F(N_1,N_2) + T (s_1 + s_2)\), where \(\partial_{N_i}F=\mu^i(N_1,N_2)\). The phase space integral \(Z(e)\) then gives: \[\begin{align} Z(e)&=\int_\Omega d\Omega\, e^{s_1 + s_2}\delta(e-E(s_1,s_2,N_1,N_2)) \\ &\propto e^{e/T}\int_{N_1,N_2} dN_1dN_2 \,e^{-F(N_1,N_2) /T}=e^{e/T}\; Z_{N}. \end{align}\] The explicit form of the equilibrium distribution yields \[\begin{align} \label{ED95exchanger} \mathcal{P}(s_1,s_2,N_1,N_2) &=\mathcal{P}_0(E(s_1,s_2,N_1,N_2))\frac{e^{s_1+s_2}}{Z(E)} \nonumber\\ &\propto \mathcal{P}_0(E(s_1,s_2,N_1,N_2))\frac{e^{-F(N_1,N_2)/T}}{Z_N}. \end{align}\tag{47}\] Analogously to the example in Sec. 3.3, this distribution favors states that minimize the total chemical free energy, subject to energy conservation. Assuming a uniform initial energy distribution, the prefactor can be absorbed into the normalization constant, giving the finite-dimensional analog distribution of Refs. [58], [59]. The complete derivations supporting this section are provided in Appendix 8.

3.5 Mechanical open system↩︎

In this section, we present an example that illustrates the applicability of the framework to systems far from equilibrium. Specifically, we consider a particle driven by an external agent through a force \(\mathbf{f}_+ (\mathbf{x},\mathbf{p},\lambda)\), immersed in a heat bath undergoing external heat transfer. A key distinction from the previous examples is that the system is not constrained to relax to a steady state and may remain far from equilibrium throughout its evolution. The corresponding Lagrangian is given by \[L(\mathbf{x}, \mathbf{v}, s)=\frac{1}{2}m|\mathbf{v}|^2-U(\mathbf{x},s).\] The variational formulation proceeds by taking variations of the (Hamilton–Pontryagin modified) action, subject to the nonholonomic Lagrange–d’Alembert constraints 28 and incorporating the principle of virtual work to account for the external manipulation: \[\begin{align} &\delta \int_{t_1}^{t_2} dt \left( L (\mathbf{x},\mathbf{v},s) + \mathbf{p} \cdot ( \mathbf{\dot{x}} - \mathbf{v}) + \dot{\Gamma}(s-\Sigma) \right) \\ & +\int_{t_1}^{t_2} dt\; \mathbf{f}_+\cdot \delta \mathbf{x}=0, \\ &\begin{cases} \displaystyle\frac{\partial L}{\partial s}\dot{\Sigma}= (\mathbf{f_-}+ \underline{\underline{\mathbf{g}}}\circ\boldsymbol{\zeta}) \cdot\mathbf{\dot{x}}+(\mathcal{J}_{s,h}+\kappa\circ\eta)\cdot(\dot{\Gamma}-T^h)\\ \hfill +\mathcal{I}_p\, \dot{\Gamma}, \\ \displaystyle\frac{\partial L}{\partial s}\delta \Sigma= (\mathbf{f_-}+ \underline{\underline{\mathbf{g}}}\circ\boldsymbol{\zeta}) \cdot\delta \mathbf{x}+(\mathcal{J}_{s,h}+\kappa\circ\eta)\cdot\delta\Gamma \\ \hfill +\mathcal{I}_p\, \delta\Gamma. \end{cases} \end{align}\] The dissipative mechanical forces are interpreted as in Sec. 3.3. We denote the thermal displacement by \(\Gamma\), and the corresponding external entropy flux by \(\mathcal{J}_{s,h}(\mathbf{x}, \mathbf{p}, s)\), with its stochastic counterpart given by \(\kappa(\mathbf{x}, \mathbf{p}, s) \circ \eta(t)\). The external thermodynamic affinity at the boundary is represented by \(T^h\), the temperature of the external reservoir. This should not be confused with the case of an open thermodynamic system with prescribed fluxes. Additionally, we incorporate the entropy pumping due to external manipulation, denoted by \(\mathcal{I}_p\), as defined in Eq. 29 .

We restrict to 1-dimension for convenience, giving the EoMs: \[\begin{align}\label{eq:EoM95open} &p=m\dot{x},\\ & \displaystyle\dot{p}=-\partial_x U-\gamma(x,p,s)\frac{p}{m}+{\rm f}_+ + g(x,p,s)\circ\zeta,\\ & \displaystyle\dot{s} =\frac{-1}{T(x,s)}\left( -\gamma(x,p,s)\frac{p}{m}+g(x,p,s)\circ\zeta\right)\cdot \frac{p}{m} \\ &\phantom{\dot{s} =} \displaystyle +\frac{T^h}{T}\cdot(\mathcal{J}_{s,h}+\kappa\circ\eta), \\ & \displaystyle\dot{\Sigma} =\frac{-1}{T(x,s)}\left( -\gamma(x,p,s)\frac{p}{m}+g(x,p,s)\circ\zeta\right)\cdot \frac{p}{m} \\ & \phantom{\dot{\Sigma} =} \displaystyle +\frac{T^h - T}{T}\cdot(\mathcal{J}_{s,h}+\kappa\circ\eta)-\mathcal{I}_p,\\ &\dot{s}_f =\dot{s}-\dot{\Sigma} = \mathcal{J}_{s,h}+\kappa\circ\eta+\mathcal{I}_p . \end{align}\tag{48}\] Since the system is open, the medium EPR \(\dot{\Sigma}\) includes an additional contribution to the thermodynamic entropy, accounting for entropy exchanges with sources outside the system \(\Omega\). We denote this contribution by \(\dot{s}_f\), representing the entropy flow entering or leaving the system. See Fig. 3 for an illustration of the system described by Eqs. 48 .

Regarding the first law we find \[\begin{align} \label{eq:1law95open} \dot{E} &=\left(\dot{p}+\partial_x U\right)\cdot \frac{p}{m}+\partial_s U\dot{s}={\rm f}_+\cdot \frac{p}{m}+(\mathcal{J}_{s,h}+\kappa\circ\eta)T^h \nonumber\\ &=P^{ext}_W + P^{ext}_{H}. \end{align}\tag{49}\] Here, \(P^{ext}_W,P^{ext}_{H}\) denote the external power input due to mechanical work and heat transfer, respectively. Equation 49 represents the stochastic extension of its deterministic counterpart, as given in Eq. (22) of [51].

Denoting the mechanical and thermal dissipative currents as \[\label{eq:currents95open} \begin{eqnarray} && { j}_d= \left(-\gamma\frac{p}{m}-g\partial_p g+\frac{p}{m}g\partial_s (g/T)\right)\mathcal{P}\nonumber\\ && \phantom{ { j}_d=} -g^2\partial_p \mathcal{P}+\frac{pg^2}{mT}\partial_s \mathcal{P}, \\ && j_{th} = \left[\mathcal{J}_{s,h}-\kappa\partial_s\left(\frac{T^h}{T}\kappa\right)\right]\mathcal{P}-\frac{T^h}{T}\kappa^2\frac{\partial \mathcal{P}}{\partial s}, \end{eqnarray}\tag{50}\] the average system EPR gives \[\begin{align} \dot{\mathcal{S}}=& \int \frac{d\Omega}{\mathcal{P}}\,\left(\frac{ { j}_d^2}{g^2}+\frac{j_{th}^2}{\kappa^2}\right) +\langle \partial_p {\rm f}_+ \rangle \\ & -\int d\Omega\, \frac{ { j}_d}{\;g^2} \left(-\gamma\frac{p}{m}-g\partial_p g+\frac{p}{m}g\partial_s (g/T)\right)\\ & -\int d\Omega\,\frac{j_{th}}{\kappa^2}\left[\mathcal{J}_{s,h}-\kappa\partial_s\left(\frac{T^h}{T}\kappa\right)\right]. \end{align}\] Given that \(\mathcal{I}_p = \partial_p {\rm f}_+\), the average medium EPR yields \[\begin{align} \langle\!\langle\dot{\Sigma} \rangle\!\rangle &=\langle\!\langle\dot{s} \rangle\!\rangle-\langle\!\langle\dot{s}_f \rangle\!\rangle \\ &= \int d\Omega\, \left(-\frac{p}{mT} { j}_d+ \frac{T^h-T}{T}j_{th}\right)-\langle \partial_p {\rm f}_+ \rangle. \end{align}\] Gathering both contributions and enforcing the second law the FDRs are obtained \[\begin{align} &\dot{S}_{tot} = \dot{\mathcal{S}}+ \langle\!\langle\dot{\Sigma} \rangle\!\rangle= \int \frac{d\Omega}{\mathcal{P}}\,\left(\frac{ { j}_d^2}{D}+\frac{j_{th}^2}{\kappa^2}\right) \geq 0 \iff \\ &\begin{cases} \displaystyle -\gamma\frac{p}{m}=-\frac{D}{T}\frac{p}{m}+g\partial_p g-\frac{p}{m}g\partial_s (g/T),\\ \displaystyle \mathcal{J}_{s,h}=\kappa^2\left(\frac{T^h -T}{T}\right)+\kappa\partial_s\left(\frac{T^h}{T}\kappa\right). \end{cases} \label{eq:FDR95open} \end{align}\tag{51}\] The FDR for the mechanical DoFs coincides with that of a closed mechanical system discussed in Sec. 3.3, see Eq. 34 . This is because, in the present example, we have not included cross-effects between the dissipative force and the external heat flux, although such couplings could, in principle, be incorporated. Additionally, the FDR associated with the external heat flux resembles that of interconnected systems introduced in Sec. 3.4, see Eq. 43 . In the macroscopic limit, it again converges to the results of nonequilibrium thermodynamics, i.e., the FDR 51 is consistent with its macroscopic counterpart given by [51] \[I:= -\frac{1}{T}\left(-\gamma\frac{p}{m}\right)\cdot \frac{p}{m}+\mathcal{J}_{s,h}\left(\frac{T^h -T}{T}\right)\geq 0.\]

However, the implications of these FDRs extend further. By applying the path integral (PI) formulation to the system described by Eqs. 48 , one can compute the ratio of probabilities for forward and backward trajectories in the thermodynamic phase space \(\Omega\). The inverse of the covariance matrix is of the form \[\begin{align} D_{\mu\nu}=\begin{bmatrix} \displaystyle\frac{1}{g^2}+\frac{(p/m)^2}{\kappa^2(T^h)^2} & \displaystyle\frac{T(p/m)}{\kappa^2(T^h)^2}\\ \displaystyle\frac{T(p/m)}{\kappa^2(T^h)^2} & \displaystyle\frac{T^2}{\kappa^2(T^h)^2} \end{bmatrix}. \end{align}\] Let us define the even and odd vectors \[\begin{align} \boldsymbol{u}_+ &= \begin{bmatrix} \dot{p} +\partial_x U - {\rm f}_+\\ \displaystyle\frac{p}{mT}(-\gamma\frac{p}{m}-\nu_p)- \frac{T^h}{T}( \mathcal{J}_{s,h}-\kappa\partial_s(\kappa T^h /T)) \end{bmatrix}, \\ \boldsymbol{u}_- &= \begin{bmatrix} \displaystyle\gamma p/m+\nu_p \\ \dot{s} \end{bmatrix}, \end{align}\] with the noise-induced drift \(\nu_p = g\partial_p g-\frac{p}{m}g\partial_s (g/T)\). We denote \(F^p ,F^s\) as the deterministic forces in the momentum and entropy Eqs. 48 . Decomposing \(D_{\mu\nu}\) into the symmetric and anti-symmetric contributions \(D_{\mu\nu}^+ ,D_{\mu\nu}^-\), respectively, the OM action and its time reversal are given by \[\begin{align} \mathcal{A}[\mathbf{u}] =& \int_{t_1}^{t_2}d\tau\left( \frac{1}{4}(\mathbf{u_+}+\mathbf{u_-})^T D_{\mu\nu}(\mathbf{u_+}+\mathbf{u_-}) \right.\nonumber\\ & \left. +\frac{1}{2}(\partial_p F^p + \partial_s F^s)\right),\\ \hat{\mathcal{A}}[\mathbf{u}] =& \int_{t_1}^{t_2}d\tau\left( \frac{1}{4}(\mathbf{u_+}-\mathbf{u_-})^T (D_{\mu\nu}^+ - D_{\mu\nu}^-)(\mathbf{u_+}-\mathbf{u_-})\right. \nonumber\\ & \left. +\frac{1}{2}(-\partial_p \hat{F}^p + \partial_s \hat{F}^s)\right), \end{align}\] where \(\mathbf{u}=(x,p,s)\). The ratio of logarithms 10 then gives \[\begin{align} \hat{\mathcal{A}}[\mathbf{u}]- \mathcal{A}[\mathbf{u}] =& \int_{t_1}^{t_2}d\tau \left(\frac{1}{4}\left( -2\mathbf{u_+}D_{\mu\nu}^-\mathbf{u_+}-2\mathbf{u_-}D_{\mu\nu}^-\mathbf{u_-}\right.\right. \nonumber\\ &-\left.4\mathbf{u_+}D_{\mu\nu}^+ \mathbf{u_-}\right) -\partial_pF_+^p - \partial_s F_-^s \bigg ). \end{align}\] Thus, after some algebra, \[\begin{align} \Delta\mathfrak{s}_m &= \int_{t_1}^{t_2}d\tau\left( \frac{1 }{g^2}\left( -\gamma\frac{p}{m}+g\circ\zeta\right)\left(-\gamma\frac{p}{m}-\nu_p \right)\right. \nonumber\\ & \left. + \frac{1}{\kappa^2}(\mathcal{J}_{s,h}+\kappa\circ\eta)( \mathcal{J}_{s,h}-\kappa\partial_s(\kappa T^h /T)) - \partial_pF_+^p \right), \end{align}\] where we denoted the medium entropy defined through the PI in the extended phase space \(\Omega\) as \(\mathfrak{s}_m\), in order to distinguish it from \(s_m\) in the restricted phase space of Sec. 2.2. Differentiating and enforcing the FDRs 51 , we find \[\begin{align} \dot{\mathfrak{s}}_m =& \frac{d}{d\tau}\text{ln}\frac{P[\mathbf{u}|\mathbf{u}_0 ]}{P[\hat{\mathbf{u}}|\mathbf{u}_t ]} \nonumber \\ =& -\frac{1}{T}\left( -\gamma\frac{p}{m}+g\circ\zeta\right)\cdot \frac{p}{m} + \frac{T^h -T}{T}(\mathcal{J}_{s,h}+\kappa\circ\eta) \nonumber\\ &- \partial_p {\rm f}_+ , \end{align}\] which comparing with Eqs. 48 gives \(\dot{\mathfrak{s}}_m = \dot{\Sigma}\). This result is remarkable for a number of reasons. Firstly, it verifies the consistency of our interpretation of \(\dot{\Sigma}\) as the medium EPR. Additionally, it shows how the variational approach, specifically the FDRs, provides a well-defined physical meaning to \(\dot{\mathfrak{s}}_m\), by satisfying the generalized Clausius relation 25 . Moreover, this result implies that \(\Delta s_{tot} = \Delta s_{sys}+\Delta\Sigma\), with \[\Delta s_{sys}=\text{ln}\frac{\mathcal{P}(\mathbf{u}_0, t_0)}{\mathcal{P}(\mathbf{u}_t,t)} ,\quad \Delta \Sigma =\text{ln}\frac{P[\mathbf{u}|\mathbf{u}_0 ]}{P[\hat{\mathbf{u}}|\mathbf{u}_t ]},\] is compatible, or can be derived, through the KL divergence 7 , when the system is constrained by the FDRs. As a consequence, we can prove the following relation, for arbitrary (normalized) initial and final distributions \(\mathcal{P}_0 ,\mathcal{P}_t\) [1]: \[\begin{align} \label{eq:masterFT} \langle e^{-\ln \mathcal{P}_0/\mathcal{P}_t - \Delta\Sigma} \rangle &= \langle e^{- \Delta\Sigma} \mathcal{P}_t /\mathcal{P}_0 \rangle \nonumber \\ &= \int \mathcal{D}\mathbf{u}\,\mathcal{P}_0 P[\mathbf{u}|\mathbf{u}_0]e^{- \Delta\Sigma} \mathcal{P}_t /\mathcal{P}_0 \nonumber \\ & = \int \mathcal{D}\mathbf{u}\, P[\hat{\mathbf{u}}|\mathbf{u}_t ]\mathcal{P}_t = 1, \end{align}\tag{52}\] where the last equality follows due to normalization, analogously as for the IFT, see 8 in Sec. 2.2. The relation 52 , referred to as the master fluctuation theorem (master FT), serves as a unifying result from which many known FTs in ST can be derived as special cases; see Ref. [1] for a detailed discussion. This demonstrates that a significant portion of the ST literature is encompassed within the variational framework presented here, rather than constituting an independent framework. Notably, the resulting FTs remain valid even in the presence of nonlinear couplings between thermal and mechanical DoFs, thereby broadening the domain of applicability of current formulations. Full derivations supporting this section are provided in Appendix 9.

A technical remark is in order: when applying the same computation to systems lacking external irreversible processes, such as in the example of Sec. 3.3, one often encounters a singular covariance matrix due to degeneracies arising from energy conservation. In these cases, the master FT can still be recovered by employing the Moore–Penrose pseudo-inverse [60], as illustrated in Appendix 10. Consequently, the above results remain valid for closed systems as well.

4 Discussion↩︎

Through the examples presented, we have illustrated the application of the variational formulation of stochastic thermodynamics and shown how, by construction, enforcing the second law yields a modeling framework that is fully thermodynamically consistent in the sense defined in Sec. 2.4. Specifically, requiring the non-negativity of the average total entropy production rate ensures that two key conditions are satisfied: the thermodynamic entropy is well-defined through a generalized Clausius relation; and for isolated systems, the Fokker–Planck equation admits an equilibrium steady-state distribution consistent with detailed balance. These properties are not imposed ad hoc but instead emerge naturally from the variational construction.

Furthermore, although derived from a geometrically motivated ansatz, the dynamical equations governing entropy production—including those at the level of individual trajectories—retain the same form as in conventional stochastic thermodynamics, up to a reformulation in an extended phase space. Remarkably, this formulation has been shown to coincide with the Kullback–Leibler divergence approach, thereby ensuring consistency with fluctuation theorems and compatibility with the broader theoretical landscape.

In Secs. 3.3 and 3.4, we showed that the equilibrium distribution arises as a special case of a general equilibrium functional, which recovers classical results—such as the Maxwell-Boltzmann distribution—in appropriate limits. For mechanical systems, the corresponding FDRs generalize the standard formulations by enabling a new class of couplings between thermal and mechanical DoFs, in which the system parameters explicitly depend on the heat bath macrostate \(s\). In the free energy variational framework of Ref. [61], it is shown how a Legendre transform of the action shifts this entropy dependence into an explicit temperature dependence of the parameters, allowing systematic inclusion of cases where phenomenological coefficients—such as viscosity—vary with temperature. Such couplings become particularly relevant in systems where this temperature dependence is non-negligible, as demonstrated in e.g. Ref. [62].

We emphasize that, due to the intrinsic dependence of the probability distribution on thermodynamic entropy, steady states of the Fokker-Planck equation necessarily correspond to equilibrium states, as they imply a stationary (time-independent) entropy. Nonetheless, this does not contradict previous findings on nonequilibrium steady states; the key difference is that those steady states refer to marginal distributions where the entropy state variable has been integrated out.

A further advantage of this approach is its capacity to model complex systems as interconnections of primitive subsystems. As demonstrated in Sec. 3.4, this modular construction naturally yields FDRs for coupled mass and heat transfer that define symmetric, positive-definite matrices, analogous to Onsager’s matrix [57]. In the presence of cross-correlated noise, these matrices consistently capture cross-effects in accordance with nonequilibrium thermodynamic theory.

The variational formulation also extends to open systems operating far from equilibrium. In the macroscopic limit, the obtained FDRs are consistent with phenomenological results from classical nonequilibrium thermodynamics [7], [51]. Importantly, this formalism allows model parameters—such as temperature, transport coefficients or friction—to be state-dependent in the complete phase space \(\Omega\), thus offering significant flexibility. Nevertheless, it is important to remark the limitation of this construction which relies on a Langevin description, and therefore assumes a separation of time scales between the slow dynamics of interest and the fast DoFs that have been integrated out.

A number of technical subtleties arise in specific limits. For example, while we have focused on the underdamped regime to retain the full dynamical structure, taking the overdamped limit within the variational framework must be done with care. In particular, the operations of enforcing the second law and applying the overdamped approximation do not, in general, commute: the overdamped approximation should be performed prior to the second law constraint. However, even then, the overdamped approximation does not provide a generally reliable account of entropy production and may fail to capture key contributions [63].

Furthermore, when multiplicative noise is present, the noise-induced drift terms lead to FDRs that take the form of partial differential equations rather than algebraic constraints. While closed-form solutions may be obtainable in specific cases—as in Ref. [52] for Coulomb friction—this is not guaranteed in general. In such cases, the model can be systematically extended by incorporating these drifts into an effective dissipative force (e.g., \({\rm f}_- \rightarrow {\rm f}_- + \nu_p\)), resulting in an algebraic FDR. This strategy has been employed in various contexts [5], [6], [14], and fits naturally within the variational formulation developed here.

Finally, when applying this framework to models that omit relevant energy sources, inconsistencies manifest in the FDR. For example, in systems exhibiting negative friction, the FDR yields \(g^2/T = \gamma < 0\), which is incompatible with \(g^2 \geq 0\) unless one assumes a negative temperature \(T < 0\). This apparent violation can be resolved by recognizing that such effective energy injection arises from the coupling to an additional degree of freedom.

5 Conclusion↩︎

In summary, the results presented here show that the variational formulation of stochastic thermodynamics provides a unified and systematic geometric framework for constructing models that are thermodynamically consistent by design. Crucially, this approach remains valid far from equilibrium, applies to both closed and open systems, and naturally accommodates nonlinear system-bath interactions beyond weak coupling. It systematically integrates microscopic fluctuations with macroscopic thermodynamic structure, offering a flexible and modular formalism grounded in the second law. These features make it broadly applicable, from isolated mechanical systems to open, interacting, and actively driven environments.

Remarkably, fluctuation-dissipation relations have emerged as fundamental constraints for thermodynamic consistency beyond detailed balance. In particular, they have been shown to ensure compliance with the second law of thermodynamics and to provide a well-defined interpretation of entropy production. Within this context, the variational approach provides a systematic framework for deriving such relations from first principles.

Ongoing work extends this framework to infinite-dimensional systems (fields), where the underlying geometric structure of the variational principle enables thermodynamically consistent modeling of complex fluids, active matter, and other continuum systems. Promising research directions include the development of variational integrators for structure-preserving numerical schemes [64]–[66], the application of Lagrangian reduction by symmetry to capture emergent phenomena [67], [68], and the incorporation of non-Markovian dynamics within this formalism.

6 Two-step average proof↩︎

In order to prove \[\langle\!\langle g(\mathbf{x},\mathbf{p}) \mathbf{\dot{p}} \rangle\!\rangle=\int d\mathbf{x}d\mathbf{p}\, g(\mathbf{x},\mathbf{p})\mathbf{j}_p,\] we only need to prove \(\langle \dot{\mathbf{p}}|\mathbf{x},\mathbf{p},t \rangle=\mathbf{j}_p (\mathbf{x},\mathbf{p},t)/\mathcal{P}(\mathbf{x},\mathbf{p},t)\), since \[\begin{align} \langle\!\langle g(\mathbf{x},\mathbf{p}) \mathbf{\dot{p}} \rangle\!\rangle&=\langle g(\mathbf{x},\mathbf{p}) \langle \dot{\mathbf{p}}|\mathbf{x},\mathbf{p},t \rangle \rangle\\ &= \int d\mathbf{x}d\mathbf{p}\,\mathcal{P}(\mathbf{x},\mathbf{p},t)g(\mathbf{x},\mathbf{p})\frac{\mathbf{j}_p (\mathbf{x},\mathbf{p},t)}{\mathcal{P}(\mathbf{x},\mathbf{p},t)}\\ &= \int d\mathbf{x}d\mathbf{p}\,g(\mathbf{x},\mathbf{p})\mathbf{j}_p . \end{align}\] We also restrict to the 1-dimensional case for clarity. Following the approach in [1], we take the 2nd order central difference consistent with Stratonovich convention: \[\begin{align} \label{eq:center95difference} \langle \dot{p}|x,p,t \rangle=&&\lim_{\epsilon\rightarrow 0}\frac{1}{2\epsilon}\left[\langle p(t+\epsilon)-p(t)|x,p,t \rangle\right. \nonumber\\ && + \left.\langle p(t)-p(t-\epsilon)|x,p,t \rangle\right], \end{align}\tag{53}\] where the average over paths is performed over one single time step \(\Delta t=\epsilon\) through the discretized propagator 6 , conditioned on \((x(t)=x,p(t)=p)\). We assume \(\dot{p}\) is taken at some fixed position, such that we can treat \(x\) as a constant. The change of variables to \(\Delta p =p(t+\epsilon)-p(t)\) simplifies the conditional probability through the PI integral to \[\begin{align} & P(p(t+\epsilon)=p+\Delta p\,|\,p(t)=p) \\ & =\frac{1}{\sqrt{4 \pi g^2 \epsilon}}\exp\left[ -\left\{\frac{1}{4 g^2\epsilon}\left(\Delta p-(f-gg')\epsilon\right)^2+\frac{\epsilon}{2}f' \right\} \right]\\ &= \frac{1}{\sqrt{4 \pi g^2 \epsilon}}\exp\left[ -\left\{\left(\alpha-\beta\right)^2+\frac{\epsilon}{2}f' \right\} \right], \end{align}\] with \(\alpha=\Delta p/\sqrt{4 g^2 \epsilon}\) and \(\beta=(f-gg')\epsilon/\sqrt{4 g^2 \epsilon}\). Hence, the first term of Eq. 53 yields: \[\begin{align} &\langle p(t+\epsilon)-p(t)\,|\,x,p,t \rangle= \int_\mathbb{R} d(\Delta p)\Delta p\, P(p+\Delta p|p)\\ & = \exp\left(-\frac{\epsilon}{2}f'\right)\sqrt{\frac{4g^2\epsilon}{\pi }}\int_\mathbb{R}d\alpha\,\alpha e^{ -\left(\alpha-\beta\right)^2}\\ &=\exp\left(-\frac{\epsilon}{2}f'\right)\sqrt{4 g^2 \epsilon}\beta= \exp\left(-\frac{\epsilon}{2}f'\right)(f-gg')\epsilon. \end{align}\] where the properties of the gaussian integral have been applied [25].

For the second term of Eq. 53 , the endpoint conditioning \(p(t)=p\) is applied using Bayes theorem [1], [18]. Since the associated Langevin equation is a Markov process, we apply Bayes theorem in the context of Markov chains as [69]: \[\begin{align} &P(p(t-\epsilon)=p-\Delta p\,|\,p(t)=p)\\ &=P(p(t)=p\,|\,p(t-\epsilon)=p-\Delta p) \frac{\mathcal{P}(x,p-\Delta p,t-\epsilon)}{\mathcal{P}(x,p,t)}\\ &\approx P(p(t)=p|p(t-\epsilon)=p-\Delta p)\\ &\qquad \times\left(1-\Delta p \partial_p \ln\mathcal{P}(x,p,t)-\mathcal{O}(\epsilon) \right) , \end{align}\] where we Taylor expanded \(\mathcal{P}(x,p-\Delta p,t-\epsilon)\). Thus, \[\begin{align} &\langle p(t)-p(t-\epsilon)|x,p,t \rangle \\ & = \int_\mathbb{R} d(\Delta p)\Delta p\, P(p|p-\Delta p)\left(1-\Delta p \partial_p \ln\mathcal{P}(x,p,t)-\mathcal{O}(\epsilon) \right)\\ & = (1-\mathcal{O}(\epsilon) )\int_\mathbb{R} d(\Delta p)\Delta p\, P(p|p-\Delta p) \\ &\quad - \partial_p \ln\mathcal{P}(x,p,t)\int_\mathbb{R} d(\Delta p)\Delta p^2\, P(p|p-\Delta p)\\ & = \exp\left(-\frac{\epsilon}{2}f'\right)\left[(f-gg')\epsilon- \partial_p \ln\mathcal{P}(x,p,t)2g^2\epsilon - \mathcal{O}(\epsilon^2) \right], \end{align}\] where the properties of the gaussian integral have been again applied. Adding all terms and taking the limit: \[\begin{align} \langle \dot{p}|x,p,t \rangle&=\lim_{\epsilon\rightarrow 0}\frac{1}{2\epsilon}\exp\left(-\frac{\epsilon}{2}f'\right)\left[2(f-gg')\epsilon\right. \\ &\quad \left.- \partial_p \ln\mathcal{P}(x,p,t)2g^2\epsilon-\mathcal{O}(\epsilon^2)\right]\\ &=(f-gg')- g^2\partial_p\ln\mathcal{P}(x,p,t)\\ &= \frac{1}{\mathcal{P}}\left[(f-g\partial_pg)\mathcal{P}- g^2 \partial_p \mathcal{P}\right]=\frac{j_p(x,p,t)}{\mathcal{P}(x,p,t)}, \end{align}\] where the last equality is given by definition 4 , and the proof is completed.

7 Mechanical isolated system derivations↩︎

In this section the detailed derivations of Sec. 3.3 are provided. First, we derive the EoMs 32 . The variational principle in this case is of the form: \[\label{equation95q95s} \begin{align} &\delta \int_{t_1}^{t_2} dt\left( L (x,v,s) + p \cdot ( \dot{x} - v) + \dot{\Gamma}(s-\Sigma)\right) =0 , \\ &\begin{cases} \displaystyle\frac{\partial L}{\partial s}\dot{\Sigma}= ({\rm f}_-+ g\circ\zeta) \cdot \dot{x}, \\ \displaystyle\frac{\partial L}{\partial s}\delta \Sigma= ({\rm f}_-+ g\circ\zeta) \cdot\delta x. \end{cases} \end{align}\tag{54}\] Taking variations of the Lagrangian and taking into account the fixed endpoints condition gives \[\begin{align} \int_{t_1}^{t_2}dt\left( \frac{\partial L}{\partial x}\delta x +\frac{\partial L}{\partial v}\delta v +\frac{\partial L}{\partial s}\delta s +\delta p(\dot{x} -v) - \dot{p} \delta x - p\delta v \nonumber \right.\\ - \delta\Gamma(\dot{s} - \dot{\Sigma}) + \dot{\Gamma}(\delta s - \delta\Sigma)\bigg)=0, \end{align}\] such that, applying the nonholonomic constraints, \[\begin{eqnarray} && \delta p:v=\dot{x},\\ && \delta v: p = \frac{\partial L}{\partial v} = mv,\\ && \delta\Gamma:\dot{s} = \dot{\Sigma},\\ && \delta s:\dot{\Gamma} = -\frac{\partial L}{\partial s} = T,\\ && \delta x: \dot{p} = \frac{\partial L}{\partial x} + {\rm f}_-+ g\circ\zeta, \end{eqnarray}\] which finally gives the system 32 .

To this the system corresponds the FP equation \(\partial_t \mathcal{P}= -\nabla_\Omega \cdot \mathbf{j}\) with currents \[\begin{align} j_x =& \mathcal{P}\frac{p}{m},\\ j_p =& \left(-\partial_x U-\gamma\frac{p}{m}-g\partial_p g+\frac{p}{m}g\partial_s (g/T)\right)\mathcal{P}\nonumber\\ & -g^2\partial_p \mathcal{P}+\frac{pg^2}{mT}\partial_s \mathcal{P}, \\ j_s =& -\frac{p}{mT} { j}_d. \end{align}\] We denote the reversible and irreversible currents as \({ j}_r= -\partial_x U \mathcal{P}\) and \({ j}_d= j_p - { j}_r\), respectively. The average system EPR gives, applying natural boundary conditions 16 : \[\begin{align} \dot{\mathcal{S}}&=-\frac{d}{dt}\int d\Omega\, \mathcal{P}\ln \mathcal{P}=-\int d\Omega\, \ln \mathcal{P}\partial_t \mathcal{P}\\ & = \int d\Omega\, \ln \mathcal{P}\nabla_\Omega\cdot \mathbf{j} = -\int d\Omega\, \frac{\mathbf{j}}{\mathcal{P}}\cdot\nabla_\Omega\mathcal{P}\\ &=-\int d\Omega \left[\frac{p}{m}\partial_x \mathcal{P}-\partial_x U\partial_p \mathcal{P}+\frac{ { j}_d}{\mathcal{P}}(\partial_p\mathcal{P}-\frac{p}{mT}\partial_s\mathcal{P})\right]\\ & =\int d\Omega\left[ \mathcal{P}\left(\partial_x\frac{p}{m} -\partial_p\partial_x U \right) -\frac{ { j}_d}{\mathcal{P}}(\partial_p\mathcal{P}-\frac{p}{mT}\partial_s\mathcal{P})\right]\\ &= \int d\Omega\, \frac{ { j}_d}{\mathcal{P}g^2} \left[ { j}_d-\left(-\gamma\frac{p}{m}-g\partial_p g+\frac{p}{m}g\partial_s (g/T)\right)\mathcal{P}\right]. \end{align}\] The average thermodynamic EPR can be derived by either applying the two-step average \(\langle\!\langle\dot{s} \rangle\!\rangle = \int d\Omega \,j_s\) or from the ensemble description: \[\begin{align} \frac{d\langle s \rangle}{dt} &= \int d\Omega\, s \partial_t \mathcal{P}= -\int d\Omega\, s \nabla_\Omega\cdot\mathbf{j} = \int d\Omega\, \mathbf{j}\cdot \nabla_\Omega s \nonumber\\ &= \int d\Omega\, j_s=-\int d\Omega\, \frac{p}{mT} { j}_d=\langle\!\langle\dot{s} \rangle\!\rangle. \end{align}\] The total EPR then yields: \[\begin{align} \dot{S}_{tot} &= \dot{\mathcal{S}}+\langle\!\langle\dot{s} \rangle\!\rangle= \int d\Omega\, \frac{ { j}_d^2}{\mathcal{P}g^2} \\ &-\int d\Omega\, { j}_d\left\{\frac{p}{mT}+\frac{1}{g^2}\left(-\gamma\frac{p}{m}-g\partial_p g+\frac{p}{m}g\partial_s (g/T)\right)\right\}, \end{align}\] and the FDR 34 is found by setting to zero the expression enclosed by the brackets \(\{...\}\).

8 Interconnected systems derivations↩︎

In this section the detailed derivations of Sec. 3.4 are provided. Firstly, taking variations of the Lagrangian and taking into account the fixed endpoints condition gives \[\begin{align} \int_{t_1}^{t_2}dt\sum_i\left( \frac{\partial L}{\partial N_i}\delta N_i +\frac{\partial L}{\partial s_i}\delta s_i+\dot{W}^i\delta N_i - \dot{N}_i\delta W^i \nonumber \right.\\ - \delta\Gamma^i(\dot{s}_i - \dot{\Sigma}_i) + \dot{\Gamma}^i(\delta s_i - \delta\Sigma_i)\bigg)=0. \end{align}\] Applying the nonholonomic constraints, \[\begin{align} & \delta N_i: -\frac{\partial L}{\partial N_i}=\dot{W}^i =\mu^i,\\ & \delta s_i:-\frac{\partial L}{\partial s_i} =\dot{\Gamma} ^i= T^i,\\ & \delta\Gamma^i:\dot{s}_i = \dot{\Sigma}_i,\\ & \delta W^i: \dot{N}_i=(\mathcal{J}_{j\rightarrow i}+\sigma_{j\rightarrow i}\circ\zeta), \end{align}\] which finally gives the system 40 . To this system corresponds the FP with probability currents given by Eqs. 41 . The full expression of the mass and heat fluxes including the noise-induced drifts \(\nu_m ,\nu_{th}\) are

We now derive the total EPR. Firstly, the average thermodynamic EPR can be derived by either applying the two-step average or from the ensemble description, recalling the natural boundary conditions 16 : \[\begin{align} & \frac{d}{dt}\sum_k\langle s_k \rangle = \int d\Omega\, (s_1+s_2) \partial_t \mathcal{P}= -\int d\Omega\, (s_1+s_2) \nabla_\Omega\cdot\mathbf{j} \nonumber\\ &= \int d\Omega\, \mathbf{j}\cdot \nabla_\Omega (s_1+s_2) = \int d\Omega\, (j_{s_1}+j_{s_2}) =\sum_k\langle\!\langle\dot{s}_k \rangle\!\rangle \nonumber\\ &= \int d\Omega\left\{j_{N_1}\left(\frac{\mu^2}{T^2}-\frac{\mu^1}{T^1}\right)+j_{th}\left(\frac{T^1-T^2}{T^1}+\frac{T^2-T^1}{T^2}\right) \right\}. \end{align}\] Then, the system EPR is computed as follows,

where we used Eqs. (41 , 42 ). Adding both contributions and recalling the matrices \(\mathbb{L},\mathbb{L}^{-1}\) as defined in 44 , the average total EPR gives the form, with the shorthand notation \(\boldsymbol{\mathfrak{J}}=(j_{N_2},j_{th})\): \[\begin{align} \dot{S}_{tot} &= \int \frac{d\Omega}{\mathcal{P}}\boldsymbol{\mathfrak{J}}^T\mathbb{L}^{-1}\boldsymbol{\mathfrak{J}} \nonumber\\ & - \int d\Omega\,\boldsymbol{\mathfrak{J}}\cdot\left\{ \mathbb{L}^{-1} \begin{bmatrix} \tilde{\mathcal{J}}_{12} \\ \tilde{J}_{12} \end{bmatrix}- \begin{bmatrix} \displaystyle\frac{\mu^1}{T^1}-\frac{\mu^2}{T^2} \\ \displaystyle\frac{T^1-T^2}{T^1}+\frac{T^2-T^1}{T^2} \end{bmatrix}\right\}. \end{align}\] The FDR 43 is then found by setting to zero the expression enclosed by the brackets \(\{...\}\). The matrix \(\mathbb{L}\) has eigenvalues \[\lambda_{\pm}=\frac{1}{2}\left(\kappa^2+\sigma^2\pm\sqrt{(\sigma^2-\kappa^2)^2+4C^2\kappa^2\sigma^2}\right),\] which satisfy \(\lambda_{\pm}>0\) for \(|C|<1\), ensuring that \(\mathbb{L}\) is positive definite. In this case, its inverse \(\mathbb{L}^{-1}\) is also positive definite, since its eigenvalues are \(1/\lambda_{\pm}>0\), which guarantees that \(\dot{S}_{tot} \geq 0\). Recall that the extreme case \(C=\pm 1\) render the model ill-defined due to redundancies.

9 Mechanical open system derivations↩︎

In this section the detailed derivations of Sec. 3.5 are provided. First, we derive the EoMs 48 . The variational principle in this case is of the form: \[\begin{align} &\delta \int_{t_1}^{t_2}dt\left( L (x,v,s) + p \cdot ( \dot{x} - v) + \dot{\Gamma}(s-\Sigma)\right) \\ & +\int_{t_1}^{t_2}dt\; {\rm f}_+\cdot \delta x =0, \\ &\begin{cases} \displaystyle\frac{\partial L}{\partial s}\dot{\Sigma}= (f_-+ g\circ \zeta) \cdot\dot{x}+(\mathcal{J}_{s,h}+\kappa\circ\eta)\cdot(\dot{\Gamma}-T^h)\\ \hfill +\mathcal{I}_p\, \dot{\Gamma}, \\ \displaystyle\frac{\partial L}{\partial s}\delta \Sigma= (f_-+ g\circ \zeta) \cdot\delta x+(\mathcal{J}_{s,h}+\kappa\circ\eta)\cdot\delta\Gamma \\ \hfill +\mathcal{I}_p\, \delta\Gamma. \end{cases} \end{align}\] Taking variations of the Lagrangian and taking into account the fixed endpoints condition gives \[\begin{align} \int_{t_1}^{t_2}dt\left( \frac{\partial L}{\partial x}\delta x +\frac{\partial L}{\partial v}\delta v +\frac{\partial L}{\partial s}\delta s +\delta p(\dot{x} -v) - \dot{p} \delta x - p\delta v \nonumber \right.\\ - \delta\Gamma(\dot{s} - \dot{\Sigma}) + \dot{\Gamma}(\delta s - \delta\Sigma) + {\rm f}_+\cdot \delta x\bigg)=0, \end{align}\] such that, applying the nonholonomic constraints, \[\begin{eqnarray} && \delta p:v=\dot{x},\\ && \delta v: p = \frac{\partial L}{\partial v} = mv,\\ && \delta s:\dot{\Gamma} = -\frac{\partial L}{\partial s} = T,\\ && \delta\Gamma:\dot{s} - \dot{\Sigma} = \mathcal{J}_{s,h}+\kappa\circ\eta +\mathcal{I}_p,\\ && \delta x: \dot{p} = \frac{\partial L}{\partial x} + {\rm f}_+ + {\rm f}_-+ g\circ\zeta, \end{eqnarray}\] which finally gives the system 48 .

To this the system corresponds the FP equation \(\partial_t \mathcal{P}= -\nabla_\Omega\cdot \mathbf{j}\) with currents \[\begin{align} j_x &= \mathcal{P}\frac{p}{m},\\ j_p &= \left(-\partial_x U+ {\rm f}_+ \right)\mathcal{P}+ j_d, \\ j_s &= -\frac{p}{mT} { j}_d+ \frac{T^h}{T}j_{th}, \end{align}\] with \(j_d ,j_{th}\) as defined in Eqs. 50 . Since in this case \(\dot{s} \neq \dot{\Sigma}\), the second law relies on the average medium EPR, see Eq. 30 , for which we apply \(\langle \dot{p}|x,p,s,t \rangle=j_p /\mathcal{P}\) and \(\langle \dot{s}|x,p,s,t \rangle=j_s /\mathcal{P}\): \[\begin{align} \langle\!\langle\dot{\Sigma} \rangle\!\rangle&=\langle\!\langle\dot{s} \rangle\!\rangle-\langle\!\langle\dot{s}_f \rangle\!\rangle \\ &= \int d\Omega\, \left(-\frac{p}{mT} { j}_d+ \frac{T^h-T}{T}j_{th}\right)-\langle \partial_p {\rm f}_+ \rangle. \end{align}\] The average system EPR gives, applying natural boundary conditions 16 :

Adding both contributions the total average EPR gives \[\begin{align} \dot{S}_{tot} &= \dot{\mathcal{S}}+ \langle\!\langle\dot{\Sigma} \rangle\!\rangle\\ &= \int \frac{d\Omega}{\mathcal{P}}\,\left(\frac{ { j}_d^2}{g^2}+\frac{j_{th}^2}{\kappa^2}\right)\\ & -\int d\Omega\, { j}_d\left\{\frac{p}{mT}+\frac{1}{g^2}\left(-\gamma\frac{p}{m}-g\partial_p g+\frac{p}{m}g\partial_s (g/T)\right)\right\}\\ & +\int d\Omega\,j_{th}\left\{\left(\frac{T^h-T}{T}\right)-\frac{1}{\kappa^2}\left[\mathcal{J}_{s,h}-\kappa\partial_s\left(\frac{T^h}{T}\kappa\right)\right]\right\}. \end{align}\] The FDRs 51 are found by enforcing the second law, setting to zero the expressions enclosed by the brackets \(\{...\}\).

10 Master FT for closed systems↩︎

The covariance matrix of Eqs. 32 corresponds to \[D^{\nu\mu}=\begin{bmatrix} g^2 &\displaystyle -\frac{pg^2}{mT}\\ \displaystyle -\frac{pg^2}{mT} & \displaystyle \left(\frac{pg}{mT}\right)^2 \end{bmatrix},\] such that \(\det D^{\nu\mu}=0\) and hence the inverse does not exist. However, we can derive the Moore–Penrose pseudo-inverse [60] \[D_{\mu\nu}^\dagger =\frac{1}{g^2(T^2+(p/m)^2)^2}\begin{bmatrix}T^4 &\displaystyle -T^3 \frac{p}{m}\\ \displaystyle-T^3 \frac{p}{m} & \displaystyle \left(T\frac{p}{m}\right)^2 \end{bmatrix},\] and apply it in the OM action to compute the medium EPR—see Ref. [70] for a rigorous treatment of OM functionals involving pseudo-inverses in Gaussian measures. We denote the even and odd contributions \[\begin{align} \boldsymbol{u}_+ &= \begin{bmatrix} \dot{p} +\partial_x U \\ \displaystyle\frac{p}{mT}(-\gamma\frac{p}{m}-\nu_p) \end{bmatrix}, \\ \boldsymbol{u}_- &= \begin{bmatrix} \displaystyle\gamma p/m+\nu_p \\ \dot{s} \end{bmatrix}, \end{align}\] with \(\nu_p = g\partial_p g-\frac{p}{m}g\partial_s (g/T)\). Decomposing \(D_{\mu\nu}^\dagger\) into the symmetric and anti-symmetric contributions \(D_{\mu\nu}^{\dagger,+} ,D_{\mu\nu}^{\dagger,-}\), respectively, the OM action and its time reversal are given by \[\begin{align} \mathcal{A}[\mathbf{u}] =& \int_{t_1}^{t_2} d\tau \bigg(\frac{1}{4}(\mathbf{u_+}+\mathbf{u_-})^T D_{\mu\nu}^\dagger (\mathbf{u_+}+\mathbf{u_-}) \nonumber\\ & +\frac{1}{2}(\partial_p F^p + \partial_s F^s)\bigg),\\ \hat{\mathcal{A}}[\mathbf{u}] =& \int_{t_1}^{t_2} d\tau \bigg(\frac{1}{4}(\mathbf{u_+}-\mathbf{u_-})^T (D_{\mu\nu}^{\dagger,+} - D_{\mu\nu}^{\dagger,-} )(\mathbf{u_+}-\mathbf{u_-}) \nonumber\\ & +\frac{1}{2}(-\partial_p \hat{F}^p + \partial_s \hat{F}^s)\bigg). \end{align}\] The ratio of logarithms 10 then gives (recall the measure contribution cancels): \[\begin{align} \hat{\mathcal{A}}[\mathbf{u}]- \mathcal{A}[\mathbf{u}] =& \int_{t_1}^{t_2} d\tau\bigg(\frac{1}{4}\left( -2\mathbf{u_+}D_{\mu\nu}^{\dagger,-} \mathbf{u_+}-2\mathbf{u_-}D_{\mu\nu}^{\dagger,-} \mathbf{u_-}\right. \nonumber\\ &-\left.4\mathbf{u_+}D_{\mu\nu}^{\dagger,+} \mathbf{u_-}\right) -\partial_pF_+^p - \partial_s F_-^s \bigg), \end{align}\] which after some algebra simplifies to \[\begin{align} \Delta\mathfrak{s}_m = \int_{t_1}^{t_2} d\tau \,\frac{1 }{g^2}\left( -\gamma\frac{p}{m}+g\circ\zeta\right)\left(-\gamma\frac{p}{m}-\nu_p \right) . \end{align}\] Differentiating and enforcing the FDR 34 , we find \[\begin{align} \dot{\mathfrak{s}}_m &= \frac{d}{d\tau}\text{ln}\frac{P[\mathbf{u}|\mathbf{u}_0 ]}{P[\hat{\mathbf{u}}|\mathbf{u}_t ]} = -\frac{1}{T}\left( -\gamma\frac{p}{m}+g\circ\zeta\right)\cdot \frac{p}{m} \nonumber\\ &= \dot{\Sigma}. \end{align}\] The last equality holds by definition, as \(\dot{s} = \dot{\Sigma}\) for closed systems; see Eqs. 32 . Consequently, \(\Delta\Sigma\) is consistent with the logarithmic ratio of path probabilities under the constraint of the FDR, and therefore satisfies the master FT 52 , as established by the argument presented in Sec. 3.5. The same result holds upon inclusion of an external manipulation \({\rm f}_+\), demonstrating that the master FT remains valid for both closed and open systems, including those influenced by external irreversible processes or manipulations.

References↩︎