Linear response in large deviations theory: a method to compute non-equilibrium distributions

Equilibrium statistical mechanics provides the probability distribution over the states of a system in terms of their energy. Finding an explicit form that generalizes the Gibbs distribution to non-equilibrium stationary conditions has been the object of extensive research in the past and remains a central issue in non-equilibrium physics nowadays [1–11]. The efforts in this direction may be divided into two research lines. The first is based on expanding the non-equilibrium probability around the known equilibrium distribution. Close to equilibrium, the correction to the Gibbs distribution was first related to the dissipation of the system by McLennan [3]. Since the dissipation is an extensive quantity of the trajectory duration and is quadratic in the distance from equilibrium, more recent work [7] has been devoted to clarify which part of the dissipation enter the McLennan formula and the suitable ordering of the long time and near to equilibrium limits. Far from equilibrium, the stationary distribution has been obtained as formal series either based on all the cumulants of the (transient) entropy production or involving the excess dynamical activity induced by the non-equilibrium forces [8]. Even at linear order around equilibrium, the difficulty in obtaining explicit results lies in the calculations of conditional averages over the full set of possible trajectories [6, 12]. The second research line is based on Freidlin–Wentzell large deviations (LDs) theory [13], which states that the logarithm of the stationary distribution of systems with weak noise is proportional to a rate function, called the quasi-potential. Since explicit computations of the rate function are in general very difficult, expansions around a reference quasi-potential have been developed [14, 15]. However, mainly systems affected by Gaussian noise were considered [14, 15], and without leveraging the thermodynamic structure of the underlying dynamics.

In this paper we combine the previous two approaches by considering thermodynamically consistent stochastic dynamics displaying a macroscopic limit in which the probability distribution satisfies a LDs principle. Then, the steady state probability of a state x can be approximated by P_ss( x ) ∝ exp(−Ωf( x )), for a sufficiently large scale parameter Ω, where f( x ) is the rate function or quasi-potential associated to the steady state. For a system in equilibrium at inverse temperature β, the usual Gibbs distribution indicates that f( x ) must be the free energy density ϕ( x ) times β (assuming the free energy is extensive in Ω). Our first result is an explicit expression to compute the first order non-equilibrium correction g( x ) to the steady state rate function. We do so for systems subjected to Poissonian noise, i.e. microscopically described by a Markov jump process, thus complementing the LDs approach and providing a clear thermodynamic interpretation of the results. We show that g( x ) is fully determined by the deterministic dynamics, described by a drift vector field u ( x ), and the scalar field $\dot {W}(\boldsymbol{x})$ that indicates what is the rate at which work is performed on the system for a given state x . This result offers a practical method to determine g( x ), which is applied to exactly solvable problems in electronics and chemistry, and is shown to be more accurate than usual linear response (LR) at the level of the probability distribution [16–21]. It can also be generalized to transient evolutions and non-autonomous settings. Second, we extend our formalism and derive a general LR formula for the quasi-potential around an arbitrary reference state. Finally, we exploit these results to devise a virtual dynamics of the quasi-potential in a generic parameter space that, upon integration, allows one to calculate non-equilibrium rate functions starting from the known system free energy. This method considerably simplifies the complexity of the original problem, which amounts to the solution of an infinite-order algebraic equation in the gradient of the quasi-potential.

We consider a Markovian jump process on a discrete space $\boldsymbol{n}\in {\mathbb{N}}^{k}$ with time-independent transition rates λ_ρ ( n ), associated to jumps n → n + Δ_ρ . Namely, the time evolution of the probability distribution P( n , t) over states is given by the master equation

$$\begin{equation}{\partial }_{t}P(\boldsymbol{n},t)=\sum\limits _{\rho }\left[{\lambda }_{\rho }(\boldsymbol{n}-{\mathbf{\Delta }}_{\rho })P(\boldsymbol{n}-{\mathbf{\Delta }}_{\rho },t)-{\lambda }_{\rho }(\boldsymbol{n})P(\boldsymbol{n},t)\right].\end{equation} \tag{ 1 }$$

We assume that a macroscopic limit exists: there is a scaling parameter Ω such that the transition rates scale as Ω and that the typical values of the density x = n /Ω are finite for Ω → ∞. Important examples satisfying the previous assumptions include chemical reaction networks [22–28] in the large volume limit, electronic circuits [29, 30] in the limit of large capacitances (see example in section 8.1), and some coarse-grained models of interacting many-body systems [31]. Then, in the limit Ω → ∞, the probability distribution P( x , t) satisfies a LDs principle [32, 33]:

$$\begin{equation}P(\boldsymbol{x},t)\asymp {\text{e}}^{-{\Omega}f(\boldsymbol{x},t)}.\end{equation} \tag{ 2 }$$

Indeed, plugging the previous ansatz in the master equation and keeping only the dominant terms in Ω → ∞, we obtain the following dynamical equation for the rate function f( x , t) (see appendix A for details):

$$\begin{equation}{\partial }_{t}f(\boldsymbol{x},t)=\sum\limits _{\rho }{\omega }_{\rho }(\boldsymbol{x})\left[1-{\text{e}}^{{\mathbf{\Delta }}_{\rho }\cdot \nabla f(\boldsymbol{x},t)}\right],\end{equation} \tag{ 3 }$$

where ω_ρ ( x ) = lim_Ω→∞  λ_ρ (Ω x )/Ω are the scaled rates and the gradient operator is with respect to x . Thus, the steady state rate function f_ss( x ) must satisfy [34]:

$$\begin{equation}0=\sum\limits _{\rho }{\omega }_{\rho }(\boldsymbol{x})\left[1-{\text{e}}^{{\mathbf{\Delta }}_{\rho }\cdot \nabla {f}_{\text{ss}}(\boldsymbol{x})}\right].\end{equation} \tag{ 4 }$$

Equation (3) is a Hamilton–Jacobi equation for the action f( x , t), with the Hamiltonian $H(\boldsymbol{x},\boldsymbol{p})={\sum }_{\rho }{\omega }_{\rho }(\boldsymbol{x})\left[1-{\text{e}}^{{\mathbf{\Delta }}_{\rho }\cdot \boldsymbol{p}}\right]$ and generalized momenta defined as usual as p = ∇f. Therefore, the solutions of equation (3) can also be found by integrating the associated Hamiltonian equations [35–37]. Equation (4) represents the stationary version, corresponding to a Hamiltonian dynamics on the manifold of constant null energy H( x , p ) = 0. The steady state rate function solving the latter equation is also known as a 'quasi-potential' or 'non-equilibrium potential'. This is in part because of the analogy between equation (2) and the usual equilibrium Gibbs distribution, but more importantly because f_ss( x ) can be shown to be a Lyapunov function of the deterministic dynamics [38], as the potential energy is in overdamped equilibrium settings (see appendix B). Away from equilibrium, the time derivative of f_ss( x ) along deterministic trajectories was related to the adiabatic/non-adiabatic decomposition of the entropy production rate [39–41] in reference [34]. Also, the quasi-potential provides information about the lifetime of non-equilibrium metastable states, in full analogy to equilibrium reaction-rate theory [30, 37, 42, 43].

Thermodynamic consistency is enforced by assuming that for each transition with rate λ_ρ ( n ) there is a corresponding reverse transition with rate λ_−ρ ( n ) (thus, Δ_−ρ = −Δ_ρ ), and that those rates are related by the local detailed balance (LDB) conditions

$$\begin{equation}\enspace \mathrm{log}\frac{{\lambda }_{\rho }(\boldsymbol{n})}{{\lambda }_{-\rho }(\boldsymbol{n}+{\mathbf{\Delta }}_{\rho })}=-\beta \left[{\Phi}(\boldsymbol{n}+{\mathbf{\Delta }}_{\rho })-{\Phi}(\boldsymbol{n})-{W}_{\rho }(\boldsymbol{n})\right],\end{equation} \tag{ 5 }$$

where Φ( n ) is the energy associated to a given state n , and W_ρ ( n ) is a non-conservative work contribution associated to the transition n → n + Δ_ρ . For simplicity we consider isothermal settings at inverse temperature β, but all of the following results hold also in more general situations where the system might interact with reservoirs at different temperatures. In that case, the LDB conditions take a form equivalent to equation (5), with the energy βΦ( n ) being replaced by a generalized Massieu potential, and the work contributions can be expressed in terms of a minimal set of fundamental forces [44].

In the macroscopic limit we will assume that the energy of the system is an extensive quantity and therefore the limit ϕ( x ) = lim_Ω→∞ Φ( n = Ω x )/Ω is well defined. For the work contributions we will abuse notation and write W_ρ ( x ) = lim_Ω→∞  W_ρ ( n = Ω x ). Thus, we consider situations in which the work contributions associated to a transition do not scale with the system size, as is indeed the case in electronic circuits [29, 30] and chemical reaction networks [23–28]. According to the previous definitions, for large Ω the LDB conditions become:

$$\begin{equation}\enspace \mathrm{log}\frac{{\omega }_{\rho }(\boldsymbol{x})}{{\omega }_{-\rho }(\boldsymbol{x})}=-\beta \left[{\mathbf{\Delta }}_{\rho }\cdot \nabla \phi (\boldsymbol{x})-{W}_{\rho }(\boldsymbol{x})\right].\end{equation} \tag{ 6 }$$

From equation (2) we see that as we increase the scale parameter Ω, the distribution P( x , t) becomes increasingly localized around the most probable states, i.e. the mimina of the rate function f( x , t). If we let x (t) be one of the minima of f(⋅, t) for a given time t, then it is possible to show that it evolves according to the following closed deterministic dynamics [22, 30, 37]:

$$\begin{equation}{\mathrm{d}}_{t}\boldsymbol{x}(t)=\boldsymbol{u}(\boldsymbol{x}(t)),\end{equation} \tag{ 7 }$$

where the deterministic drift u ( x ) is given by

$$\begin{equation}\boldsymbol{u}(\boldsymbol{x})=\sum\limits _{\rho }\hspace{2pt}{\enspace \omega }_{\rho }(\boldsymbol{x})\enspace {\mathbf{\Delta }}_{\rho }=\sum\limits _{\rho { >}0}\hspace{2pt}{\enspace I}_{\rho }(\boldsymbol{x})\enspace {\mathbf{\Delta }}_{\rho }.\end{equation} \tag{ 8 }$$

In the last equality we have exploited the fact that, due to the requirement of thermodynamic consistency, every transition has a reversed associated one, and expressed the drift u ( x ) in terms of the average forward currents I_ρ ( x ) = ω_ρ ( x ) − ω_−ρ ( x ). Indeed, a second order truncation in the Kramers–Moyal expansion of the master equation in equation (1) leads to a Fokker–Planck equation with the vector field u ( x ) defined above as the drift. However, note that such expansion is not globally valid, but only locally around the minima of f(⋅, t) [45]. An alternative proof of the previous result, based solely on equation (3), is given in appendix B, where we also show that the steady state rate function satisfying equation (4) is a Lyapunov function for the deterministic dynamics [38].

In this and the following section we will consider the perturbation of equation (4) around the equilibrium steady state. For this, it is useful to define ${\omega }_{\rho }^{(0)}(\boldsymbol{x})$ to be the scaled rates in the absence of non-conservative sources of work (i.e. when W_ρ ( x ) = 0). According to equation (6), they satisfy

$$\begin{equation}\mathrm{log}\frac{{\omega }_{\rho }^{(0)}(\boldsymbol{x})}{{\omega }_{-\rho }^{(0)}(\boldsymbol{x})}=-\beta {\mathbf{\Delta }}_{\rho }\cdot \nabla \phi (\boldsymbol{x}).\end{equation} \tag{ 9 }$$

In that case, the dynamics of the system is said to be (unconditionally) detailed balanced, and for long times a thermal equilibrium state will be reached, which is exactly given by the Gibbs distribution P_eq( n ) ∝ exp(−βΦ( n )). Thus, by equation (2), the stationary rate function corresponding to this situation is f_ss( x ) = βϕ( x ). In a general situation with non-vanishing work contributions W_ρ ( x ), we can always consider the decomposition

$$\begin{equation}{f}_{\text{ss}}(\boldsymbol{x})=\beta \phi (\boldsymbol{x})+g(\boldsymbol{x}),\end{equation} \tag{ 10 }$$

where g( x ) quantifies the deviations with respect to thermal equilibrium. We will show in what follows that, to first order in the quantities βW_ρ ( x ), g( x ) must satisfy the equation

$$\begin{equation}{\boldsymbol{u}}^{(0)}(\boldsymbol{x})\cdot \nabla g(\boldsymbol{x})=-\beta {\dot {W}}^{(0)}(\boldsymbol{x}),\end{equation} \tag{ 11 }$$

where the vector field u ⁽⁰⁾ ( x ) and the scalar field ${\dot {W}}^{(0)}(\boldsymbol{x})$ are given by

$$\begin{equation}{\boldsymbol{u}}^{(0)}(\boldsymbol{x})=\sum\limits _{\rho { >}0}\hspace{2pt}{\enspace I}_{\rho }^{(0)}(\boldsymbol{x}){\mathbf{\Delta }}_{\rho }\qquad {\dot {W}}^{(0)}(\boldsymbol{x})=\sum\limits _{\rho { >}0}\hspace{2pt}{\enspace I}_{\rho }^{(0)}(\boldsymbol{x}){W}_{\rho }(\boldsymbol{x}),\end{equation} \tag{ 12 }$$

in terms of the scaled currents

$$\begin{equation}{I}_{\rho }^{(0)}(\boldsymbol{x})={\omega }_{\rho }^{(0)}(\boldsymbol{x})-{\omega }_{-\rho }^{(0)}(\boldsymbol{x}).\end{equation} \tag{ 13 }$$

Before proceeding with the proof of equation (11), we discuss the interpretation of the previous expressions. In the first place, ${I}_{\rho }^{(0)}(\boldsymbol{x})$ is just the average current along transition ρ according to the detailed balanced dynamics, given that the state of the system is x . Then, the vector field u ⁽⁰⁾ ( x ) is the net drift in state space for state x , as discussed in the previous section. As a consequence, the field lines associated with u ⁽⁰⁾ ( x ) are just the deterministic trajectories of the detailed balanced system. The quantity ${I}_{\rho }^{(0)}(\boldsymbol{x}){W}_{\rho }(\boldsymbol{x})$ is the average scaled work rate associated to the transition ρ in the state x (to lower order in W_ρ ( x )), and ${\dot {W}}^{(0)}(\boldsymbol{x})$ is thus the total scaled work rate (the actual work rate is ${\Omega}{\dot {W}}^{(0)}(\boldsymbol{x})$).

Equation (11) displays an apparent 'gauge invariance', meaning that if g(x) is any given solution, then $~{g}(\boldsymbol{x})=g(\boldsymbol{x})+h(\boldsymbol{x})$ will also be a solution provided that u ⁽⁰⁾ ( x ) ⋅ ∇h( x ) = 0. However, this indeterminacy is removed if one demands that the function g( x ) must be a single-valued function for all x . This can be seen by considering a situation in which the deterministic dynamics has a unique fixed point x * (i.e. u ( x *) = 0). Then, for any initial condition x ₀, the ensuing trajectory x (t| x ₀) satisfies lim_t→+∞  x (t| x ₀) = x *. Also,

$$\begin{equation}{\mathrm{d}}_{t}g(\boldsymbol{x}(t\vert {\boldsymbol{x}}_{0}))={\mathrm{d}}_{t}\boldsymbol{x}(t\vert {\boldsymbol{x}}_{0})\cdot \nabla g(\boldsymbol{x}(t\vert {\boldsymbol{x}}_{0}))={\boldsymbol{u}}^{(0)}(\boldsymbol{x}(t\vert {\boldsymbol{x}}_{0}))\cdot \nabla g(\boldsymbol{x}(t\vert {\boldsymbol{x}}_{0}))=-\beta {\dot {W}}^{(0)}(\boldsymbol{x}(t\vert {\boldsymbol{x}}_{0})),\end{equation} \tag{ 14 }$$

which gives, once integrated over time,

$$\begin{equation}g({\boldsymbol{x}}_{0})=g({\boldsymbol{x}}^{{\ast}})+\beta {\int }_{0}^{+\infty } \mathrm{d}t\enspace {\dot {W}}^{(0)}(\boldsymbol{x}(t\vert {\boldsymbol{x}}_{0})).\end{equation} \tag{ 15 }$$

Thus, fixing the value g( x *) at the fixed point, the previous equation allows to determine g( x ₀) for any other point x ₀. This can be generalized to situations with multistability at equilibrium, where the deterministic dynamics has multiple fixed points. In that case one can construct local solutions g( x ) for each basin of attraction in the same way as before, fixing an arbitrary value for g( x ) at each fixed point. These values must be later adjusted in order for the different local solutions to match continuously with each other at the separatrices between different basins of attraction (the existence of a continuous quasi-potential for diffusive systems with multiple coexisting attractors is dicussed in [35]).

We now proceed with the proof of equation (11). Since we are assuming thermodynamic consistency, we can rewrite equation (4) as:

$$\begin{align}0& =\sum\limits _{\rho { >}0}\left\{{\omega }_{\rho }(\boldsymbol{x})\left[1-{\text{e}}^{{\mathbf{\Delta }}_{\rho }\cdot \nabla {f}_{\text{ss}}(\boldsymbol{x})}\right]+{\omega }_{-\rho }(\boldsymbol{x})\left[1-{\text{e}}^{-{\mathbf{\Delta }}_{\rho }\cdot \nabla {f}_{\text{ss}}(\boldsymbol{x})}\right]\right\}\\ & =\sum\limits _{\rho { >}0}\left[{\omega }_{\rho }(\boldsymbol{x})-{\omega }_{-\rho }(\boldsymbol{x}){\text{e}}^{-{\mathbf{\Delta }}_{\rho }\cdot \nabla {f}_{\text{ss}}(\boldsymbol{x})}\right]\left[1-{\text{e}}^{{\mathbf{\Delta }}_{\rho }\cdot \nabla {f}_{\text{ss}}(\boldsymbol{x})}\right].\end{align} \tag{ 16 }$$

Now we use the LDB conditions in equation (6) to write:

$$\begin{equation}0=\sum\limits _{\rho { >}0}\hspace{2pt}{\enspace \omega }_{\rho }(\boldsymbol{x})\left[1-{\text{e}}^{\beta \left[{\mathbf{\Delta }}_{\rho }\cdot \nabla \phi (\boldsymbol{x})-{W}_{\rho }(\boldsymbol{x})\right]}\enspace {\text{e}}^{-{\mathbf{\Delta }}_{\rho }\cdot \nabla {f}_{\text{ss}}(\boldsymbol{x})}\right]\left[1-{\text{e}}^{{\mathbf{\Delta }}_{\rho }\cdot \nabla {f}_{\text{ss}}(\boldsymbol{x})}\right].\end{equation} \tag{ 17 }$$

We see from the last expression that when W_ρ ( x ) = 0, the choice f_ss( x ) = βϕ( x ) is indeed a solution, and in fact it makes each of the terms in the sum vanish independently. Thus, considering f_ss( x ) = βϕ( x ) + g( x ), we have

$$\begin{equation}0=\sum\limits _{\rho { >}0}\hspace{2pt}{\enspace \omega }_{\rho }(\boldsymbol{x})\left[1-{\text{e}}^{-\beta {W}_{\rho }(\boldsymbol{x})}\enspace {\text{e}}^{-{\mathbf{\Delta }}_{\rho }\cdot \nabla g(\boldsymbol{x})}\right]\left[1-{\text{e}}^{\beta {\mathbf{\Delta }}_{\rho }\cdot \nabla \phi (\boldsymbol{x})}\enspace {\text{e}}^{{\mathbf{\Delta }}_{\rho }\cdot \nabla g(\boldsymbol{x})}\right].\end{equation} \tag{ 18 }$$

Close to equilibrium W_ρ ( x ) is small by definition and g( x ) is assumed to be small and of the same order. It is therefore justified to expand equation (18) to first order in both of them. Doing that, we notice that the first factor in square brackets becomes linear in W_ρ ( x ) and g( x ), and therefore, to first order, all the other factors can be evaluated at zero order,

$$\begin{align}0& =\sum\limits _{\rho { >}0}\hspace{2pt}{\enspace \omega }_{\rho }^{(0)}(\boldsymbol{x})\left[\beta {W}_{\rho }(\boldsymbol{x})+{\mathbf{\Delta }}_{\rho }\cdot \nabla g(\boldsymbol{x})\right]\left[1-{\text{e}}^{\beta {\mathbf{\Delta }}_{\rho }\cdot \nabla \phi (\boldsymbol{x})}\right]\\ & =\sum\limits _{\rho { >}0}\hspace{2pt}{\enspace I}_{\rho }^{(0)}(\boldsymbol{x})\left[\beta {W}_{\rho }(\boldsymbol{x})+{\mathbf{\Delta }}_{\rho }\cdot \nabla g(\boldsymbol{x})\right],\end{align} \tag{ 19 }$$

where in the last line we used the detailed balance conditions in equation (9) to construct ${I}_{\rho }^{(0)}(\boldsymbol{x})={\omega }_{\rho }^{(0)}(\boldsymbol{x})-{\omega }_{-\rho }^{(0)}(\boldsymbol{x})$. The result of equation (11) follows immediately.

By considering a similar splitting of the time-dependent rate function f( x , t) = βϕ( x ) + g( x , t), and repeating exactly the same derivation, this time starting from the dynamical equation of equation (3), we obtain the following dynamical equation for g( x , t):

$$\begin{equation}{\partial }_{t}g(\boldsymbol{x},t)=\beta {\dot {W}}^{(0)}(\boldsymbol{x})+{\boldsymbol{u}}^{(0)}(\boldsymbol{x})\cdot \nabla g(\boldsymbol{x},t).\end{equation} \tag{ 20 }$$

We can also consider a situation with time dependent rates. For a time-dependent scaled energy ϕ( x , t), βϕ( x , t) is the rate function that defines an instantaneous equilibrium state. Considering now the splitting f( x , t) = βϕ( x , t) + g( x , t), we obtain

$$\begin{equation}{\partial }_{t}g(\boldsymbol{x},t)=\beta [{\dot {W}}^{(0)}(\boldsymbol{x},t)-{\partial }_{t}\phi (\boldsymbol{x},t)]+{\boldsymbol{u}}^{(0)}(\boldsymbol{x},t)\cdot \nabla g(\boldsymbol{x},t),\end{equation} \tag{ 21 }$$

where the time-dependent work rate ${\dot {W}}^{(0)}(\boldsymbol{x},t)$ and drift u ⁽⁰⁾ ( x , t) are defined in exactly the same way as before, but in terms of the time-dependent currents ${I}_{\rho }^{(0)}(\boldsymbol{x},t)={\omega }_{\rho }^{(0)}(\boldsymbol{x},t)-{\omega }_{-\rho }^{(0)}(\boldsymbol{x},t)$. Since the previous equation is a linear partial differential equation for g( x , t), for periodically driven systems one may also apply Floquet theory.

8.1. CMOS SRAM cell

We consider the model of a low-power static random access memory (SRAM) cell developed in [29, 30]. The usual CMOS implementation of this kind of memories involves two inverters or NOT gates connected in a loop, each of which is composed of two MOS transistors, as shown in figure 1(a). The memory is powered by applying a voltage bias 2V_dd, which takes the system out of thermal equilibrium. The stochastic model of this device has two degrees of freedom: the charges q₁ and q₂ at the output of each inverter, which correspond to voltages v_1/2 = q_1/2/C, where C is a value of capacitance characterizing the device. Two Poisson rates are associated to each of the four transistors, corresponding to forward and backward conduction events. In each conduction event, or jump, the voltages v_1/2 can change by the elementary voltage v_e = q_e/C, where q_e is the positive electron charge. Also, all Poisson rates are proportional to I₀/q_e, where I₀ is a characteristic current. Under certain physical scalings of the transistors, as discussed in [30], both I₀ and C increase linearly with the size of the device, and therefore we can take C, or equivalently ${v}_{\text{e}}^{-1}$, as the scaling parameter Ω above. Thus, the macroscopic limit Ω → ∞ must in this case be understood as the limit v_e/V_T → 0 and v_e/V_dd → 0, where ${V}_{\text{T}}={(\beta {q}_{\text{e}})}^{-1}$ is the thermal voltage. A typical steady state distribution for this system (in its bistable phase, see below) is shown in figure 1(b). The rate function f(v₁, v₂) corresponding to such steady state can be computed in terms of new variables x = (v₁ − v₂)/2 and y = (v₁ + v₂)/2. For example, the rate function of the reduced distribution for x is exactly given by [30]:

$$\begin{equation}f(x)=\frac{{x}^{2} + 2{V}_{\text{dd}}\enspace x}{{V}_{\text{T}}}+\frac{2\mathit{\text{n}}{V}_{\text{T}}}{\mathit{\text{n}}+2}\left[L(x,{V}_{\text{dd}}) - L(x,-{V}_{\text{dd}})\right],\end{equation} \tag{ 22 }$$

where $L(x,{V}_{\text{dd}})={\text{Li}}_{2}\left(-\mathrm{exp}(({V}_{\text{dd}}+x(1+2/\mathit{\text{n}}))/{V}_{\text{T}})\right)$, and Li₂(⋅) is the polylogarithm function of second order, with n a parameter characterizing the transistors. In figure 1(c) we compare the distribution obtained from the previous rate function with exact numerical results. The agreement is essentially perfect even if only a few tens of electrons are involved (the scaling parameter is only v_e = 0.1V _T , and V_dd = 1.15V_T). Also, in figure 1(c) we show that there is a transition from a monostable phase into the bistable phase that allows the storage of a bit of information, occurring at the critical powering voltage ${V}_{\text{dd}}^{{\ast}}=\mathrm{ln}(2){V}_{\text{T}}$ for n = 1. Expanding the previous expression to first order in V_dd, we obtain:

$$\begin{equation}f(x)=\underset{\beta \phi (x)}{\underbrace{\frac{{x}^{2}}{{V}_{\text{T}}}}}\underset{g(x)}{\underbrace{-4{V}_{\text{dd}}\frac{\mathit{\text{n}}}{\mathit{\text{n}}+2}\mathrm{ln}\left(2\enspace \mathrm{cosh}\left(\frac{\mathit{\text{n}}+2}{2\mathit{\text{n}}}\frac{x}{{V}_{\text{T}}}\right)\right)+\mathcal{O}({V}_{\text{dd}}^{2})}},\end{equation} \tag{ 23 }$$

where we have identified the equilibrium and non-equilibrium contributions as in equation (10). We will now recover this last result using the formalism developed here.

Zoom In Zoom Out Reset image size

We begin by writing down the deterministic equations of motion, which read:

$$\begin{equation}\begin{aligned}{\mathrm{d}}_{t}{v}_{1}& =I({v}_{1},{v}_{2})-I(-{v}_{1},-{v}_{2})\\ {\mathrm{d}}_{t}{v}_{2}& =I({v}_{2},{v}_{1})-I(-{v}_{2},-{v}_{1}),\end{aligned}\end{equation} \tag{ 24 }$$

where the scaled current I(v₁, v₂) is defined in terms of the scaled rates ω_±(v₁, v₂) as I(v₁, v₂) = ω₊(v₁, v₂) − ω₋(v₁, v₂). In turn, these rates are given by:

$$\begin{equation}\begin{aligned}{\omega }_{+}({v}_{1},{v}_{2})& =({I}_{0}/C)\enspace {\text{e}}^{({V}_{\text{dd}}-{v}_{2}-{V}_{\text{th}})/(\mathit{\text{n}}{V}_{\text{T}})}\\ {\omega }_{-}({v}_{1},{v}_{2})& ={\omega }_{+}({v}_{1},{v}_{2})\enspace {\text{e}}^{-({V}_{\text{dd}}-{v}_{1})/({V}_{\text{T}})}.\end{aligned}\end{equation} \tag{ 25 }$$

Also, because of the symmetry of the device, the work contribution of each forward elementary transition is q_e V_dd [29], and therefore the total work rate is

$$\begin{equation}\dot {W}({v}_{1},{v}_{2})={q}_{\text{e}}{V}_{\text{dd}}\left[I({v}_{1},{v}_{2})+I(-{v}_{1},-{v}_{2})+I({v}_{2},{v}_{1})+I(-{v}_{2},-{v}_{1})\right].\end{equation} \tag{ 26 }$$

To proceed, we write the equations of motion in the previously defined variables x and y:

$$\begin{equation}\begin{aligned}{\mathrm{d}}_{t}x& =[I(x,y)-I(-x,y)-I(-x,-y)+I(x,-y)]/2\\ {\mathrm{d}}_{t}y& =[I(x,y)+I(-x,y)-I(-x,-y)-I(x,-y)]/2,\end{aligned}\end{equation} \tag{ 27 }$$

where the change of variables in the function I(⋅) is implicit. We see from the previous equations that if we consider an initial condition where y(0) = 0, then the ensuing trajectory will satisfy y(t) = 0 at all times, since d_t y|_y=0 = 0 for all x. This is true, in particular, for the equilibrium dynamics (that is, for V_dd = 0) which is the one we must consider in equation (11). Thus, over a trajectory with y(0) = 0, equation (11) reduces to:

$$\begin{equation}{\mathrm{d}}_{t}{x}^{(0)}\enspace {\mathrm{d}}_{x}g(x)=-\beta {\dot {W}}^{(0)}(x),\end{equation} \tag{ 28 }$$

where the superscript (0) means in this case that the currents should be evaluated at V_dd = 0. Then, we obtain:

$$\begin{align}g(x)& =-2\frac{{V}_{\text{dd}}}{{V}_{\text{T}}}{\int }_{0}^{x}\mathrm{d}{x}^{\prime }\enspace {\left.\frac{I({x}^{\prime },0)+I(-{x}^{\prime },0)}{I({x}^{\prime },0)-I(-{x}^{\prime },0)}\right\vert }_{{V}_{\text{dd}}=0}\\ & =-4{V}_{\text{dd}}\frac{\mathit{\text{n}}}{\mathit{\text{n}}+2}\mathrm{ln}\left(2\enspace \mathrm{cosh}\left(\frac{\mathit{\text{n}}+2}{2\mathit{\text{n}}}\frac{x}{{V}_{T}}\right)\right),\end{align} \tag{ 29 }$$

in full agreement with equation (23).

In figure 2 we compare the exact rate function of equation (22) with the first order approximation obtained with our formalism. We see in figure 2(a) that while the approximation is a first order expansion in V_dd/V_T ≪ 1, it is essentially exact for values as high as V_dd/V_T = 0.4. Even more remarkably, the first order approximation is able to predict the transition to bistability, although as we see in figure 2(b) the critical point is underestimated (since at the exact critical point ${V}_{\text{dd}}^{{\ast}}=\mathrm{ln}(2){V}_{\text{T}}$ the approximation already developed a bistability). The accuracy of the first order approximation can be understood by noting that the next order vanishes by symmetry, as well as all even powers of V_dd in the expansion of equation (22). Finally, although the approximation is not expected to work at all for V_dd/V_T ≫ 1, we see in figure 2(c) that it correctly gets the most probable values (which correspond to the fixed points of the deterministic dynamics) as well as the curvature around them, while it moderately fails to describe the height of the barrier separating the two fixed points.

Zoom In Zoom Out Reset image size

8.2. Schlögl model

We consider now the Schlögl model [46], a well known autocatalytic chemical model displaying bistability far from equilibrium. It consists in the following set of chemical reactions,

$$\begin{equation}2\mathit{\text{X}}+\mathit{\text{A}}\underset{{k}_{-1}}{\overset{{k}_{+1}}{{\rightleftharpoons}}}3\mathit{\text{X}};\qquad \mathit{\text{B}}\underset{{k}_{-2}}{\overset{{k}_{+2}}{{\rightleftharpoons}}}\mathit{\text{X}}.\end{equation} \tag{ 30 }$$

Here, species A and B will be considered to be chemostated reservoirs (i.e. we fix their concentration), and the only degree of freedom for the problem is the number of molecules of species X. The reactions are assumed to take place in a well-mixed container of volume V, which can be taken to be the macroscopic parameter Ω. Thus, in the macroscopic limit we deal with the concentration x of species X, and the scaled rates (ensuing the mass-action law) read:

$$\begin{equation}{\omega }_{1}(x)={k}_{+1}a{x}^{2},\quad {\omega }_{-1}(x)={k}_{-1}{x}^{3},\quad {\omega }_{2}(x)={k}_{+2}b,\quad {\omega }_{-2}(x)={k}_{-2}x.\end{equation} \tag{ 31 }$$

We can set k₊₁ a = 1 = k₊₂ b choosing appropriate units for time and concentrations. Then, in this case the LDB conditions impose the relation k₋₁ = k₋₂e^−Δμ , where Δμ is the chemical potential difference between species A and B (measured in units such that β = 1 in the following). Thermal equilibrium is realized at Δμ = 0, for which the deterministic concentration of x is x⁽⁰⁾ = 1/k₋₁ = 1/k₋₂. The system can display bistability at Δμ ≠ 0 only for k₋₂ ≳ 1.7 (so we will set k₋₂ = 2 in the following plots). The exact rate function of the non-equilibrium steady state was obtained for example in [46] and reads:

$$\begin{equation}f(x)={\int }_{0}^{x}\mathrm{log}\left(\frac{y{k}_{-2}(1+{\text{e}}^{-{\Delta}\mu }{y}^{2})}{1+{y}^{2}}\right)\mathrm{d}y\end{equation} \tag{ 32 }$$

$$\begin{equation}=x\enspace \mathrm{log}\left(\frac{x{k}_{-2}(1+{\text{e}}^{-{\Delta}\mu }{x}^{2})}{1+{x}^{2}}\right)+2{\text{e}}^{{\Delta}\mu /2}\enspace \mathrm{arctan}\left({\text{e}}^{-\frac{{\Delta}\mu }{2}}x\right)-x-2\enspace \mathrm{arctan}(x)\end{equation} \tag{ 33 }$$

$$\begin{equation}=\underset{\beta \phi (x)}{\underbrace{x\enspace \mathrm{log}({k}_{-2}x)-x}}\underset{g(x)}{\underbrace{\hspace{2pt}+\hspace{2pt}{\Delta}\mu (\mathrm{arctan}(x)-x)+\mathcal{O}({({\Delta}\mu )}^{2})}}.\end{equation} \tag{ 34 }$$

We now recover the previous result by employing our formalism. The equilibrium currents are

$$\begin{equation}{I}_{1}^{(0)}(x)=(1-{k}_{-2}x){x}^{2},\qquad {I}_{2}^{(0)}(x)=1-{k}_{-2}x.\end{equation} \tag{ 35 }$$

The work per reactions are W₁ = Δμ and W₂ = 0, the total work rate is ${\dot {W}}^{(0)}(x)={\Delta}\mu \enspace {x}^{2}(1-{k}_{-2}x)$ and u⁽⁰⁾ (x) = (1 − k₋₂ x) (1 + x²). Then, the equation for the first order correction to the rate function is

$$\begin{equation}(1-{k}_{-2}x)(1+{x}^{2})\enspace {\mathrm{d}}_{x}g(x)=-{x}^{2}(1-{k}_{-2}x){\Delta}\mu ,\end{equation} \tag{ 36 }$$

from which we obtain $g(x)=-{\Delta}\mu \enspace {\int }_{0}^{x}{y}^{2}/(1+{y}^{2})\enspace \mathrm{d}y={\Delta}\mu [\mathrm{arctan}(x)-x]$, in full agreement with equation (34).

In figure 3(a) we compare the equilibrium rate function with the exact and first order approximation for Δμ = 0.3, showing that the first order approximation is indeed accurate. However, in contrast to the previous example, the first order approximation fails to offer an even qualitative description for values of Δμ at which the exact rate function develops a bistability (in fact, the first order approximation never develops a bistability). This model becomes monostable again for larger values of Δμ, and in that case the first order approximation correctly estimates the most probable value and the curvature around it, even if it is not expected to work in that regime, as was also observed in the previous example.

Zoom In Zoom Out Reset image size

The approach that we presented exploits the existence of a macroscopic limit with an LD principle and applies LR to the LD rate function. A natural question to ask is whether this offers any advantage with respect to usual LR theory [16–21], which directly gives the first order correction to the probability distribution (instead of the rate function associated to it), and does not require any macroscopic limit. We now show that our approach is more accurate than usual LR and valid further away from equilibrium.

Assuming that the first order correction g( x ) is already known, the LD approximation for the steady state distribution is:

$$\begin{align}{P}_{\text{ss}}(\boldsymbol{x})& \simeq \frac{1}{Z}\enspace {\text{e}}^{-{\Omega}\beta f(\boldsymbol{x})}\\ & \simeq \frac{1}{Z}\enspace {\text{e}}^{-{\Omega}\left(\beta \phi (\boldsymbol{x})+g(\boldsymbol{x})\right)}\\ & =\frac{1}{{Z}_{\text{eq}}}\enspace {\text{e}}^{-{\Omega}\beta \phi (\boldsymbol{x})}\enspace \frac{{Z}_{\text{eq}}}{Z}\enspace {\text{e}}^{-{\Omega}g(\boldsymbol{x})}\\ & \simeq {P}_{\text{eq}}(\boldsymbol{x})\enspace \frac{{Z}_{\text{eq}}}{Z}\enspace {\text{e}}^{-{\Omega}g(\boldsymbol{x})}.\end{align} \tag{ 37 }$$

In the first line of the previous equation the LD approximation is involved, while in the second line the LR approximation at the level of the rate function is involved. Also, P_eq( x ) = exp(−Ωβϕ( x ))/Z_eq is the LD approximation of the equilibrium state. The partition function Z is

$$\begin{equation}Z=\int \mathrm{d}{\boldsymbol{x}}^{\prime }\enspace {\text{e}}^{-{\Omega}\beta f({\boldsymbol{x}}^{\prime })}=\int \mathrm{d}{\boldsymbol{x}}^{\prime }\enspace {\text{e}}^{-{\Omega}\beta \phi ({\boldsymbol{x}}^{\prime })}\enspace {\text{e}}^{-{\Omega}g({\boldsymbol{x}}^{\prime })},\end{equation} \tag{ 38 }$$

and Z_eq = Z|_{g( x )=0}. For the LR approximation to be sensible we must have g( x ) ≪ βϕ( x ). However, if we consider the more stringent conditions Ωg( x ) ≪ 1, we can write

$$\begin{equation}{P}_{\text{ss}}(\boldsymbol{x})\simeq {P}_{\text{eq}}(\boldsymbol{x})\left(1+{\Omega}\int \mathrm{d}{\boldsymbol{x}}^{\prime }\enspace {P}_{\text{eq}}({\boldsymbol{x}}^{\prime })\left[g({\boldsymbol{x}}^{\prime })-g(\boldsymbol{x})\right]\right)+\mathcal{O}({{\Omega}}^{2}{g}^{2}),\end{equation} \tag{ 39 }$$

and

$$\begin{equation}Z\simeq {Z}_{\text{eq}}\left(1-{\Omega}\int \mathrm{d}{\boldsymbol{x}}^{\prime }\enspace {P}_{\text{eq}}({\boldsymbol{x}}^{\prime })g({\boldsymbol{x}}^{\prime })\right)+\mathcal{O}({{\Omega}}^{2}{g}^{2}).\end{equation} \tag{ 40 }$$

From the previous equation we see that the first order perturbation to the probabilities, δP( x ) = P_ss( x ) − P_eq( x ), is related to the first order perturbation to the LD rate function by

$$\begin{equation}\delta P(\boldsymbol{x})/{P}_{\text{eq}}(\boldsymbol{x})=-{\Omega}g(\boldsymbol{x})+c,\end{equation} \tag{ 41 }$$

where c is a constant ensuring ∫d x δP( x ) = 0. Thus, we see that δP( x ) and g( x ) contain the same information, as one can be computed from the other. However, in the macroscopic limit (Ω ≫ 1), the LR approximation at the level of the rate function is expected to be more accurate than the LR approximation at the level of the probabilities, and also valid for stronger perturbations. The reason is that for Ω ≫ 1 the conditions g( x ) ≪ βϕ( x ) might be satisfied even if the conditions Ωg( x ) ≪ 1 are not.

We now illustrate the previous results in the example of section 8.1 (CMOS SRAM cell). For a given value of the powering voltage V_dd, we will compare three probability distributions: (i) the LD approximation with the exact rate function of equation (22): P_ex(x) ∝ exp(−f(x)/v_e) (note that, as shown in figure 1(c), this distribution can be considered exact already for v_e/V_T ≃ 0.1), (ii) the LD approximation combined with the LR approximation of the rate function: P_g (x) ∝ exp(−(x²/V_T + g(x))/v_e), where g(x) is given by equation (29), and (iii) the usual LR approximation at the level of the probability distribution: P_δ (x) = P_eq(x) + δP(x) = P_eq(x) (1 − g(x)/v_e + c) (see equation (41)). In figure 4(a) we see that for V_dd/V_T = 0.2, both P_g (x) and P_δ (x) correctly describe the deviations from the equilibrium distribution. However, further away from equilibrium (in particular, at the critical point ${V}_{\text{dd}}={V}_{\text{dd}}^{{\ast}}={V}_{\text{T}}\enspace \mathrm{ln}(2)$), P_g (x) remains an acceptable approximation while P_δ (x) does not, as shown in figure 4(b). Finally, in figure 4(c) we compare the Hellinger distances H(P_g , P_ex) and H(P_δ , P_ex) between the exact distribution and the LR approximations at the level of the rate function or at the level of the probability distribution, respectively, as a function of the powering voltage V_dd (the Hellinger distance 0 ⩽ H(p, q) ⩽ 1 between distributions p(x) and q(x) is defined by ${H}^{2}(p,q)=1-{\sum }_{i}\sqrt{p({x}_{i})q({x}_{i})}$). We see that P_g (x) is always a better approximation with respect to P_δ (x), and that the gap between the two increases as we go deeper into the macroscopic limit (i.e. if we decrease the value of the elementary voltage v_e), as predicted above. We note that for sufficiently large V_dd the function P_δ (x) stops being positive definite (a common issue of the usual LR approximation at the level of the probability distribution) and therefore the Hellinger distance cannot be computed. Also, in this particular case the distance H(P_g , P_ex) reaches a maximum for large V_dd, after which it decreases again. This is related to the observation made in the previous section that the LR approximation at the level of the rate function correctly describes the most probable values and the curvature around them for large biases, even if a priori it is not expected to work in that regime.

Zoom In Zoom Out Reset image size

In this section we generalize our previous results by considering the problem of perturbing a general non-equilibrium steady state of an autonomous system. The rate function f( x ) of the steady state corresponding to some given work contributions W_ρ ( x ) must satisfy equation (17):

$$\begin{equation}0=\sum\limits _{\rho { >}0}{\omega }_{\rho }(\boldsymbol{x})\left[1-{\text{e}}^{\beta \left[{\mathbf{\Delta }}_{\rho }\cdot \nabla \phi (\boldsymbol{x})-{W}_{\rho }(\boldsymbol{x})\right]}\enspace {\text{e}}^{-{\mathbf{\Delta }}_{\rho }\cdot \nabla f(\boldsymbol{x})}\right]\left[1-{\text{e}}^{{\mathbf{\Delta }}_{\rho }\cdot \nabla f(\boldsymbol{x})}\right].\end{equation} \tag{ 42 }$$

Let us now consider the perturbation ${W}_{\rho }(\boldsymbol{x})\to {W}_{\rho }^{\prime }(\boldsymbol{x})={W}_{\rho }(\boldsymbol{x})+\delta {W}_{\rho }(\boldsymbol{x})$, and the corresponding response f( x ) → f'( x ) = f( x ) + g( x ). By expanding the previous equation to first order in δW_ρ ( x ) and g( x ), we find

$$\begin{equation}\boldsymbol{\mathcal{U}}(\boldsymbol{x})\cdot \nabla g(\boldsymbol{x})=-\beta \dot {\mathcal{W}}(\boldsymbol{x}).\end{equation} \tag{ 43 }$$

The previous equation has the same structure as equation (11) for perturbations around equilibrium, but this time the vector field $\boldsymbol{\mathcal{U}}(\boldsymbol{x})$ and the scalar field $\dot {\mathcal{W}}(\boldsymbol{x})$ are given by

$$\begin{equation}\boldsymbol{\mathcal{U}}(\boldsymbol{x})=\sum\limits _{\rho { >}0}\left({\omega }_{-\rho }(\boldsymbol{x}){\text{e}}^{-{\mathbf{\Delta }}_{\rho }\cdot \nabla f(\boldsymbol{x})}-{\omega }_{\rho }(\boldsymbol{x}){\text{e}}^{{\mathbf{\Delta }}_{\rho }\cdot \nabla f(\boldsymbol{x})}\right){\mathbf{\Delta }}_{\rho },\end{equation} \tag{ 44 }$$

and

$$\begin{equation}\dot {\mathcal{W}}(\boldsymbol{x})=\sum\limits _{\rho { >}0}{\omega }_{-\rho }(\boldsymbol{x})\left({\text{e}}^{-{\mathbf{\Delta }}_{\rho }\cdot \nabla f(\boldsymbol{x})}-1\right)\delta {W}_{\rho }(\boldsymbol{x}).\end{equation} \tag{ 45 }$$

The previous results of section 5 are easily recovered by considering f( x ) → βϕ( x ), W_ρ ( x ) → 0, and using the LDB conditions at equilibrium, since in that way we see that $\boldsymbol{\mathcal{U}}(\boldsymbol{x})$ reduces to the equilibrium deterministic drift u ( x ), and $\dot {\mathcal{W}}(\boldsymbol{x})$ to the work rate $\dot {W}(\boldsymbol{x})$. We note that $\boldsymbol{\mathcal{U}}(\boldsymbol{x})$ does not correspond to the field of deterministic trajectories, and also that $\dot {\mathcal{W}}(\boldsymbol{x})$ cannot be interpreted as the average work rate for a given state, as is the case for perturbations of the equilibrium distribution. However, it is interesting to note that the fixed points of the true unperturbed deterministic dynamics are shared by the dynamics corresponding to the alternative drift $\boldsymbol{\mathcal{U}}(\boldsymbol{x})$, and also that $\dot {\mathcal{W}}(\boldsymbol{x})$ vanishes at those fixed points. These properties are important to analyze the errors of the numerical method introduced in the next section, and are a consequence of the fact that the gradient of the exact rate function f( x ) vanishes at the deterministic fixed points (see appendix B).

In equilibrium statistical mechanics the unnormalized thermal distribution P( x |β) = exp(−βΦ( x )) at a given inverse temperature β can be obtained by solving the evolution equation d_β P( x |β) = −Φ( x )P( x |β), with a known initial condition (typically, the infinite temperature distribution P( x |0) = 1). Thus, the previous differential equation allows to obtain the distribution for an inverse temperature β + dβ given the knowledge of the distribution at inverse temperature β. Interestingly, equation (43) provides the basis for a similar procedure, in which a non-equilibrium steady state distribution is evolved starting, for example, from the equilibrium one. In order to illustrate this idea we will consider a system with a single state-independent work parameter λ (as the examples given in section 8). Then, the steady state rate function f( x |λ) has a parametric dependence on λ, and we can rewrite equation (43) as:

$$\begin{equation}\boldsymbol{\mathcal{U}}(\boldsymbol{x}\vert \lambda )\cdot \nabla [{\partial }_{\lambda }f(\boldsymbol{x}\vert \lambda )]=-\beta \dot {\mathcal{W}}(\boldsymbol{x}\vert \lambda ).\end{equation} \tag{ 46 }$$

Solving this equation we can obtain ∂_λ f( x |λ) up to an irrelevant constant, and if f( x |λ) is known, we can approximate f( x |λ + dλ) ≃ f( x |λ) + dλ ∂_λ f( x |λ). Iterating this procedure it is in principle possible to obtain the rate function f( x |λ) for arbitrary λ starting from the equilibrium distribution f( x |0).

We now apply this method to the Schlögl model discussed in section 8.2. In that case the work parameter is λ = Δμ and, since the problem is one-dimensional, we can write equation (46) as:

$$\begin{equation}{\partial }_{{\Delta}\mu }f(x\vert {\Delta}\mu )=-\beta {\int }_{0}^{x}\enspace \frac{\dot {\mathcal{W}}({x}^{\prime }\vert {\Delta}\mu )}{\mathcal{U}({x}^{\prime }\vert {\Delta}\mu )}\enspace \mathrm{d}{x}^{\prime }.\end{equation} \tag{ 47 }$$

We numerically solve the previous differential equation using a Runge–Kutta method of order 4, starting from the equilibrium rate function f(x|0) = x  log(k₋₂ x) − x (see equation (34)). In figure 5(a) we compare the exact rate function given by equation (34) with the one obtained with our evolution method up to Δμ = 2.7, a bistable point of the system. To obtain the evolved rate function shown in the plot we have sampled the intervals x ∈ [0, 10] and Δμ ∈ [0, 2.7] at 4 × 10³ points each. In the top panel of figure 5(b) we show the results obtained by performing the evolution up to Δμ = 4, this time sampling the interval x ∈ [0, 25] at 10⁴ points and the interval Δμ ∈ [0, 4] at 4 × 10³ points. The color plot indicates the value of the probability distribution P_ev(x|Δμ) ∝ e^{−Ωf(x|Δμ)} (with Ω = 10), where f(x|Δμ) is the rate function obtained with our method at each point in the evolution. For visualization purposes, it was normalized so that the maximum of P_ev(x|Δμ) for a given Δμ is always 1. In the bottom panel of figure 5(b) we show the exact probability distribution P_ex(x|Δμ) computed in the same way as before but in terms of the exact rate function in equation (34). In figure 5(c) we show the Hellinger distance H(P_ev, P_ex) as a function of Δμ. We see that the distance between the distributions computed from the evolved and exact rate functions increases rapidly at the beginning of the evolution, reaches a plateau, and increases again after the bistable region around Δμ ≃ 2.7 is crossed.

Zoom In Zoom Out Reset image size

An important comment regarding the numerical stability of the proposed method is in order. For equation (43) to have a smooth and well behaved solution (that is, such that ∇g( x ) is finite for all x ), the field $\dot {\mathcal{W}}(\boldsymbol{x})$ must vanish wherever the field $\boldsymbol{\mathcal{U}}(\boldsymbol{x})$ does. As already noted in the previous section, if x * is a fixed point of the deterministic dynamics, then from equations (44) and (45) we see that $\boldsymbol{\mathcal{U}}({\boldsymbol{x}}^{{\ast}})=0$ and $\dot {\mathcal{W}}({\boldsymbol{x}}^{{\ast}})=0$. This follows from the LDB conditions and the fact that ∇f( x *) = 0 (see appendix B). Thus, assuming that the deterministic fixed points are the only zeros of $\boldsymbol{\mathcal{U}}(\boldsymbol{x})$, the previous minimal consistency condition is indeed satisfied by equations (44) and (45). Importantly, for this to hold, the unperturbed rate function f( x ) involved in the definitions of $\boldsymbol{\mathcal{U}}(\boldsymbol{x})$ and $\dot {\mathcal{W}}(\boldsymbol{x})$ must correspond exactly to the parameters at which the rates ω_ρ ( x ) are evaluated. However, that correspondence is broken in the intermediate steps of our iterative procedure, since at each step (except at the first one when the equilibrium distribution is perturbed) the fields $\boldsymbol{\mathcal{U}}(\boldsymbol{x})$ and $\dot {\mathcal{W}}(\boldsymbol{x})$ are constructed from an approximate rate function f( x ), which might not exactly satisfy ∇f( x *) = 0. This issue might lead to errors and discontinuities for values of x close to the deterministic fixed points at each step, which will propagate forward in the virtual evolution. That is in fact the main source of errors contributing to the deviations between the exact and evolved rate functions in figure 5(a). A more sophisticated numerical implementation of the proposed method might take as an input the deterministic fixed points { x *} at each step, in order to ensure that they match the minima of the approximated rate function. This point should be considered for applications to more complex systems in many dimensions.

For thermodynamically consistent stochastic systems displaying a macroscopic limit in which the stationary distribution fulfills a LDs principle, we have shown how to compute the correction to the rate function, or quasi-potential, to first order in the forces taking the system out of thermal equilibrium. We have applied this result to two examples: a realistic model of an important kind of electronic memory and the Schlögl model of a bistable chemical reaction. We also generalized our methods to time-dependent settings. In the final part of our work we considered the LR of arbitrary non-equilibrium steady states. Based on that result, we developed a method to incrementally evolve the stationary state rate function in the space of thermodynamic forces, which in principle allows to obtain the rate function arbitrarily far from equilibrium, starting from the known equilibrium one. This might lead to new numerical schemes to study non-equilibrium fluctuations at the mesoscopic level, providing an alternative to the usual stochastic simulations based on the Gillespie algorithm.

Our work makes use of the thermodynamic consistency of the dynamics and of the existence of a macroscopic limit. The results might remain nonetheless an excellent approximation even significantly away from the strict macroscopic limit, as is the case in the examples discussed in section 8. Furthermore, we have shown that in the mesoscopic or macroscopic regime the LR approach at the level of the rate function is more accurate than traditional LR theory at the level of the probability distribution, even if both approaches are first order expansions in the thermodynamic forces. As shown in the example in section 9 (figure 4), what actually happens is that both approaches become more accurate as the scale of the system is increased, but with a larger gap between the accuracy of LR at the level of the rate function and at the level of the probability distribution. This is understood as follows. A fixed change Δ x in the density x = n /Ω requires an increasing number of elementary transitions as the scale of the system is increased. Also, the work contribution during a transition does not scale with the system size, while the associated change in free energy does. Thus, it is natural to expect LR approximations to be more accurate in the macroscopic regime, since they are based on the assumption that work contributions are small with respect to other energy scales. However, while the first order non-equilibrium correction to the rate function is independent of the scale (by definition of the rate function), the first order correction to the probability distribution P( x ) is proportional to the scale parameter (see equation (41)), and that constrains the strength of non-equilibrium forces for which the correction remains small.

We acknowledge funding from the INTER project 'TheCirco' (INTER/FNRS/20/15074473) and CORE project 'NTEC' (C19/MS/13664907), funded by the Fonds National de la Recherche (FNR, Luxembourg), and from the European Research Council, project NanoThermo (ERC-2015-CoG Agreement No. 681456).

All data that support the findings of this study are included within the article (and any supplementary files).

In this section we provide details of the derivation of equation (3) in the main text, starting from the master equation in equation (1). In the first place we represent a given discrete state n of the system by its 'density' x = n /Ω. For example, the components of n are the numbers of molecules in a chemical reaction network or the numbers of charges in an electronic circuit, while the components x are, respectively, the concentrations or the voltages. Thus, a transition n → n + Δ_ρ corresponds to a change Δ_ρ /Ω in the density x . We consider the LDs ansatz P( x , t) = exp(−Ωf( x , t))/Z(t) for the time-dependent probability distribution in terms of x , where Z(t) = ∑ _n  exp(−Ωf( n /Ω, t)) takes care of the normalization. Then, we have

$$\begin{equation}P(\boldsymbol{x}+{\mathbf{\Delta }}_{\rho }/{\Omega},t)=P(\boldsymbol{x},t){\text{e}}^{-{\mathbf{\Delta }}_{\rho }\cdot \nabla f(\boldsymbol{x})+\mathcal{O}({{\Omega}}^{-1})},\end{equation} \tag{ A1 }$$

and

$$\begin{equation}{\mathrm{d}}_{t}P(\boldsymbol{x},t)=-P(\boldsymbol{x},t)[{\Omega}\enspace {\mathrm{d}}_{t}f(\boldsymbol{x},t)+{\mathrm{d}}_{t}\enspace \mathrm{log}(Z(t))].\end{equation} \tag{ A2 }$$

Equation (3) is then obtained by inserting the previous two expressions in the master equation of equation (1), dividing both sides by ΩP( x , t), and taking the limit Ω → ∞ in which Ω⁻¹ log(Z(t)) → 0.

We define x (t) to be a minimum of the rate function f(⋅, t) at time t. Then, we have that ${\partial }_{{x}_{i}}f(\boldsymbol{x}(t),t)=0$ and that the symmetric matrix with elements ${[{\partial }_{{x}_{i},{x}_{j}}^{2}f(\boldsymbol{x}(t),t)]}_{ij}$ is positive semidefinite. For time t + dt we can write:

$$\begin{equation}\begin{aligned}0& ={\partial }_{{x}_{i}}f(\boldsymbol{x}(t+\mathrm{d}t),t+\mathrm{d}t)\\ & ={\partial }_{{x}_{i}}f(\boldsymbol{x}(t)+\boldsymbol{u}(t)\enspace \mathrm{d}t+\mathcal{O}(\mathrm{d}{t}^{2}),t+\mathrm{d}t)\\ & ={u}_{j}(t)\enspace {\partial }_{{x}_{j},{x}_{i}}^{2}f(\boldsymbol{x}(t),t)\enspace \mathrm{d}t+{\partial }_{t}{\partial }_{{x}_{i}}f(\boldsymbol{x}(t),t)\enspace \mathrm{d}t+\mathcal{O}(\mathrm{d}{t}^{2}),\end{aligned}\end{equation} \tag{ B1 }$$

where we have introduced the velocity u (t) = d_t x (t), and repeated indices are summed. Now we employ the differential equation satisfied by f( x , t), equation (3) in the main text, to obtain:

$$\begin{align}{\partial }_{t}{\partial }_{{x}_{i}}f(\boldsymbol{x},t)& ={\partial }_{{x}_{i}}{\partial }_{t}f(\boldsymbol{x},t)\\ & ={\partial }_{{x}_{i}}\left[\sum\limits _{\rho }\hspace{2pt}{\omega }_{\rho }(\boldsymbol{x})\left[1-{\text{e}}^{{({\mathbf{\Delta }}_{\rho })}_{j}{\partial }_{{x}_{j}}f(\boldsymbol{x},t)}\right]\right]\\ & =\sum\limits _{\rho }\hspace{2pt}{\partial }_{{x}_{i}}{\omega }_{\rho }(\boldsymbol{x})\left[1-{\text{e}}^{{({\mathbf{\Delta }}_{\rho })}_{j}{\partial }_{{x}_{j}}f(\boldsymbol{x},t)}\right]+\sum\limits _{\rho }\hspace{2pt}{\omega }_{\rho }(\boldsymbol{x})\left[-{\text{e}}^{{({\mathbf{\Delta }}_{\rho })}_{j}{\partial }_{{x}_{j}}f(\boldsymbol{x},t)}{({\mathbf{\Delta }}_{\rho })}_{j}{\partial }_{{x}_{i},{x}_{j}}^{2}f(\boldsymbol{x},t)\right].\end{align} \tag{ B2 }$$

Noticing that the first sum in the last line of the previous equation vanishes when evaluated at x (t), and combining the result with equation (B1), up to first order in dt we obtain:

$$\begin{equation}{\left[\boldsymbol{u}(t)-\sum\limits _{\rho }\hspace{2pt}{\omega }_{\rho }(\boldsymbol{x}(t)){\mathbf{\Delta }}_{\rho }\right]}_{j}{\partial }_{{x}_{j},{x}_{i}}^{2}f(\boldsymbol{x}(t),t)=0.\end{equation} \tag{ B3 }$$

Thus, whenever the matrix ${[{\partial }_{{x}_{i},{x}_{j}}^{2}f(\boldsymbol{x}(t),t)]}_{ij}$ can be assumed to be positive definite, we have u (t) − ∑_ρ ω_ρ ( x (t))Δ_ρ = 0, in agreement with equations (7) and (8) in the main text.

If the rate function has a single global minimum, then the instantaneous average value $\left\langle \boldsymbol{x}(t)\right\rangle ={{\Omega}}^{-1}\left\langle \boldsymbol{n}(t)\right\rangle ={{\Omega}}^{-1}{\sum }_{\boldsymbol{n}}\hspace{2pt}\boldsymbol{n}P(\boldsymbol{n},t)$ also evolves as ${\mathrm{d}}_{t}\left\langle \boldsymbol{x}(t)\right\rangle =\boldsymbol{u}(\left\langle \boldsymbol{x}(t)\right\rangle )$ in the Ω → ∞ limit. In order to see this, we multiply equation (1) in the main text by n and sum, obtaining:

$$\begin{equation}{\partial }_{t}\sum\limits _{\boldsymbol{n}}\hspace{2pt}\boldsymbol{n}P(\boldsymbol{n},t)=\sum\limits _{\rho ,\boldsymbol{n}}\hspace{2pt}\boldsymbol{n}\left[{\lambda }_{\rho }(\boldsymbol{n}-{\mathbf{\Delta }}_{\rho })P(\boldsymbol{n}-{\mathbf{\Delta }}_{\rho },t)-{\lambda }_{\rho }(\boldsymbol{n})P(\boldsymbol{n},t)\right]\end{equation} \tag{ B4 }$$

$$\begin{equation}=\sum\limits _{\rho ,\boldsymbol{n}}\left[(\boldsymbol{n}+{\mathbf{\Delta }}_{\rho }){\lambda }_{\rho }(\boldsymbol{n})P(\boldsymbol{n},t)-\boldsymbol{n}{\lambda }_{\rho }(\boldsymbol{n})P(\boldsymbol{n},t)\right]\end{equation} \tag{ B5 }$$

$$\begin{equation}=\sum\limits _{\rho ,\boldsymbol{n}}\hspace{2pt}{\mathbf{\Delta }}_{\rho }{\lambda }_{\rho }(\boldsymbol{n})P(\boldsymbol{n},t),\end{equation} \tag{ B6 }$$

where we shifted n → n + Δ_ρ in the first sum on the rhs of (B4). Dividing both sides by Ω, passing to integral over x and recognizing that P peaks around the single minimum of f as Ω → ∞ we obtain the desired result.

Finally, we mention that if f_ss( x ) is the steady state rate function, i.e. the solution to equation (4) in the main text, then it follows that [34, 38]

$$\begin{equation}{\mathrm{d}}_{t}{f}_{\text{ss}}(\boldsymbol{x}(t))=\sum\limits _{\rho }\hspace{2pt}{\omega }_{\rho }(\boldsymbol{x})\enspace {\mathbf{\Delta }}_{\rho }\cdot \nabla {f}_{\text{ss}}(\boldsymbol{x}(t)){\leqslant}\sum\limits _{\rho }\hspace{2pt}{\omega }_{\rho }(\boldsymbol{x}(t))[{\text{e}}^{{\mathbf{\Delta }}_{\rho }\cdot \nabla {f}_{\text{ss}}(\boldsymbol{x})}-1]=0,\end{equation} \tag{ B7 }$$

where we used equation (7) and x ⩽ e^x − 1. Thus, the function f_ss( x ) decreases monotonically along deterministic trajectories and, since the rate function can always be considered to be bounded by below, it is a Lyapunov function of the deterministic dynamics. As a consequence, if x * is a fixed point of the deterministic dynamics, then it is also a minimum of f_ss( x ), and therefore we have ∇f_ss( x *) = 0 (assuming f_ss( x ) is continuous). A thermodynamic interpretation of d_t ϕ( x (t)) is given in [34].