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

We consider thermodynamically consistent autonomous Markov jump processes displaying a macroscopic limit in which the logarithm of the probability distribution is proportional to a scale-independent rate function (i.e., a large deviations principle is satisfied). In order to provide an explicit expression for the probability distribution valid away from equilibrium, we propose a linear response theory performed at the level of the rate function. We show that the first order non-equilibrium contribution to the steady state rate function, $g(x)$, satisfies $u(x)\cdot \nabla g(x) = -\beta \dot W(x)$ where the vector field $u(x)$ defines the macroscopic deterministic dynamics, and the scalar field $\dot W(x)$ equals the rate at which work is performed on the system in a given state $x$. This equation provides a practical way to determine $g(x)$, significantly outperforms standard linear response theory applied at the level of the probability distribution, and approximates the rate function surprisingly well in some far-from-equilibrium conditions. The method applies to a wealth of physical and chemical systems, that we exemplify by two analytically tractable models - an electrical circuit and an autocatalytic chemical reaction network - both undergoing a non-equilibrium transition from a monostable phase to a bistable phase. Our approach can be easily generalized to transient probabilities and non-autonomous dynamics. Moreover, its recursive application generates a virtual flow in the probability space which allows to determine the steady state rate function arbitrarily far from equilibrium.


I. INTRODUCTION
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][2][3][4][5][6][7][8][9][10][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 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 large deviations 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 large deviations 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Ẇ (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 arXiv:2106.05887v2 [cond-mat.stat-mech] 11 Jun 2021 accurate than usual linear response at the level of the probability distribution [16][17][18][19][20][21]. It can also be generalized to transient evolutions and non-autonomous settings. Second, we extend our formalism and derive a general linear response 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.

II. MACROSCOPIC LIMIT AND LARGE DEVIATIONS PRINCIPLE
We consider a Markovian jump process on a discrete space n ∈ 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 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][23][24][25][26][27] in the large volume limit, electronic circuits [28,29] in the limit of large capacitances (see example in section VIII-A), and some coarse-grained models of interacting many-body systems [30]. Then, in the limit Ω → ∞, the probability distribution P (x, t) satisfies a large deviations (LD) principle [31,32]: 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): 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 [33]: Eq.
(3) is a Hamilton-Jacobi equation for the action f (x, t), with the Hamiltonian H(x, p) = ρ ω ρ (x) 1 − e ∆ρ·p and generalized momenta defined as usual as p = ∇f . Therefore, the solutions of Eq. (3) can also be found by integrating the associated Hamiltonian equations [34][35][36]. Eq. (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 Eq. (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 [37] (see Appendix B), and provides information about the lifetime of non-equilibrium metastable states, in full analogy to equilibrium reaction-rate theory [29,36,38,39].

III. THERMODYNAMIC STRUCTURE
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 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 Eq. (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 [40].
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 [28,29] and chemical reaction networks [22][23][24][25][26][27]. According to the previous definitions, for large Ω the LDB conditions become: IV. DETERMINISTIC DYNAMICS From Eq. (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: where the deterministic drift u(x) is given by 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 . Indeed, a second order truncation in the Kramers-Moyal expansion of the master equation in Eq. (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) [41]. An alternative proof of the previous result, based solely on Eq. (3), is given in Appendix B, where we also show that the steady state rate function satisfying Eq. (4) is a Lyapunov function for the deterministic dynamics [37].

V. EQUILIBRIUM STATE AND LINEAR RESPONSE
In this and the following section we will consider the perturbation of Eq. (4) around the equilibrium steady state. For this, it is useful to define ω (0) ρ (x) to be the scaled rates in the absence of non-conservative sources of work (i.e., when W ρ (x) = 0). According to Eq. (6), they satisfy log ω 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 Eq. (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 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 where the vector field u (0) (x) and the scalar fieldẆ (0) (x) are given by in terms of the currents Before proceeding with the proof of Eq. (11), we discuss the interpretation of the previous expressions. In the first place, I ρ (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 (0) (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 (0) (x) are just the deterministic trajectories of the detailed balanced system. The quantity I is the average work rate associated to the transition ρ in the state x (to lower order in W ρ (x)), andẆ (0) (x) is thus the total work rate.
Equation (11) displays an apparent 'gauge invariance', meaning that if g(x) is any given solution, theñ g(x) = g(x) + h(x) will also be a solution provided that u (0) (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 0 , the ensuing trajectory x(t|x 0 ) satisfies lim t→+∞ x(t|x 0 ) = x * . Also, which gives, once integrated over time, Thus, fixing the value g(x * ) at the fixed point, the previous equation allows to determine g(x 0 ) for any other point x 0 . 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 [34]).

VI. PROOF
We now proceed with the proof of Eq. (11). Since we are assuming thermodynamic consistency, we can rewrite Eq. (4) as: Now we use the LDB conditions in Eq. (6) to write: 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 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 Eq. (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, where in the last line we used the detailed balance conditions in Eq. (9) to construct I . The result of Eq. (11) follows immediately.

VII. DYNAMICAL EXTENSIONS
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 Eq. (3), we obtain the following dynamical equation for g(x, t): 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 ( where the time-dependent work rateẆ (0) (x, t) and drift u (0) (x, t) are defined in exactly the same way as before, but in terms of the time-dependent currents I Since the previous equation is a linear partial differential equation for g(x, t), for periodically driven systems one may also apply Floquet theory.

A. CMOS SRAM cell
We consider the model of a low-power static random access memory (SRAM) cell developed in [28,29]. 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 1 and q 2 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 0 /q e , where I 0 is a characteristic current. Under certain physical scalings of the transistors, as discussed in [29], both I 0 and C increase linearly with the size of the device, and therefore we can take C, or equivalently v −1 e , 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 T = (βq 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 1 , v 2 ) corresponding to such steady state can be computed in terms of new variables x = (v 1 − v 2 )/2 and y = (v 1 + v 2 )/2. For example, the rate function of the reduced distribution for x is exactly given by [29]: where L(x, V dd ) = Li 2 (− exp((V dd + x(1 + 2/n))/V T )), and Li 2 (·) 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 * dd = ln(2)V T for n = 1. Expanding the previous expression to first order in V dd , we obtain: −4V dd n n + 2 ln 2 cosh n + 2 2n where we have identified the equilibrium and non-equilibrium contributions as in Eq. (10). We will now recover this last result using the formalism developed here. We begin by writing down the deterministic equations of motion, which read: where the scaled current I(v 1 , v 2 ) is defined in terms of the scaled rates ω ± (v 1 , v 2 ) as In turn, these rates are given by: Also, because of the symmetry of the device, the work contribution of each forward elementary transition is q e V dd [28], and therefore the total work rate iṡ To proceed, we write the equations of motion in the previously defined variables x and y: 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 Eq. (11). Thus, over a trajectory with y(0) = 0, Eq. (11) reduces to: where the superscript (0) means in this case that the currents should be evaluated at V dd = 0. Then, we obtain: = −4V dd n n + 2 ln 2 cosh n + 2 2n in full agreement with Eq. (23).
In Figure 2 we compare the exact rate function of Eq. (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 * dd = ln(2)V T the approximation already developed a bistability). 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.

B. Schlögl model
We consider now the Schlögl model [42], a well known autocatalytic chemical model displaying bistability far from equilibrium. It consists in the following set of chemical reactions, 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: We can set k +1 a = 1 = k +2 b choosing appropriate units for time and concentrations. Then, in this case the LDB conditions impose the relation k −1 = k −2 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 (0) = 1/k −1 = 1/k −2 . The system can display bistability at ∆µ = 0 only for k −2 1.7 (so we will set k −2 = 2 in the following plots). The exact rate function of the non-equilibrium steady state was obtained for example in [42] and reads: We now recover the previous result by employing our formalism. The equilibrium currents are The work per reactions are W 1 = ∆µ and W 2 = 0, the total work rate isẆ (0) (x) = ∆µ x 2 (1 − k −2 x) and u (0) (x) = (1 − k −2 x)(1 + x 2 ). Then, the equation for the first order correction to the rate function is from which we obtain g(x) = −∆µ x 0 y 2 /(1 + y 2 ) dy = ∆µ[arctan(x) − x], in full agreement with Eq. (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 a even qualitatively description for values ∆µ at which the exact rate functions 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.

IX. COMPARISON WITH USUAL LINEAR RESPONSE
The approach that we presented exploits the existence of a macroscopic limit with a LD principle and applies linear response to the LD rate function. A natural question to ask is whether this offers any advantage with respect to usual linear response theory [16][17][18][19][20][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 linear response and valid further away from equilibrium.
Assuming that the first order correction g(x) is already known, the LD approximation for the steady state distri-bution is: In the first line of the previous equation the LD approximation is involved, while in the second line the linear response (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 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 and 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 where c is a constant ensuring dx δ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 VIII A (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 Eq. (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: is given by Eq. (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 Eq. (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 dd = V * dd = V T 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 − i 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 a consequence of 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.

X. LINEAR RESPONSE OF A NON-EQUILIBRIUM STEADY STATE
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 Eq. (17): Let us now consider the perturbation W ρ (x) → W ρ (x) = W ρ (x) + δW ρ (x), and the corresponding response . By expanding the previous equation to first order in δW ρ (x) and g(x), we find The previous equation has the same structure as Eq. (11) for perturbations around equilibrium, but this time the vector field U (x) and the scalar fieldẆ(x) are given by andẆ The previous results of section V 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 U (x) reduces to the equilibrium deterministic drift u(x), anḋ W(x) to the work rateẆ (x). We note that U (x) does not correspond to the field of deterministic trajectories, and also thatẆ(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 deterministic dynamics are shared by the dynamics corresponding to the alternative drift U (x), and also thatẆ(x) vanishes at those fixed points (see below).

XI. VIRTUAL EVOLUTION OF NON-EQUILIBRIUM STEADY STATES
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, Eq. (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 VIII). Then, the steady state rate function f (x|λ) has a parametric dependence on λ, and we can rewrite Eq. (43) as: 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 VIII B. In that case the work parameter is λ = ∆µ and, since the problem is one-dimensional, we can write Eq. (46) as: 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 −2 x) − x (see Eq. (34)). In Figure 5 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 Eq. (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. An important comment regarding the numerical stability of the proposed method is in order. For Eq. (43) to have a smooth and well behaved solution (that is, such that ∇g(x) is finite for all x), the fieldẆ(x) must vanish wherever the field U (x) does. As already noted in the previous section, if x * is a fixed point of the deterministic dynamics, then from Eqs. (44) and (45) we see that U (x * ) = 0 andẆ(x * ) = 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 U (x), the previous minimal consistency condition is indeed satisfied by Eqs. (44) and (45). Importantly, for this to hold, the unperturbed rate function f (x) involved in the definitions of U (x) andẆ(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 U (x) andẆ(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.

XII. CONCLUSION
For thermodynamically consistent stochastic systems displaying a macroscopic limit in which the stationary distribution fulfills a large deviations 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 linear response 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 VIII. Furthermore, we have shown that in the mesoscopic or macroscopic regime the linear response approach at the level of the rate function is more accurate than traditional linear response 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 IX (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 linear response 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 linear response 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 Eq. (41)), and that constrains the strength of non-equilibrium forces for which the correction remains small. 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 P (x + ∆ ρ /Ω, t) = P (x, t)e −∆ρ·∇f (x)+O(Ω −1 ) , (1), dividing both sides by ΩP (x, t), and taking the limit Ω → ∞ in which Ω −1 log(Z(t)) → 0.

Appendix B: Deterministic dynamics
We define x(t) to be a minimum of the rate function f (·, t) at time t. Then, we have that ∂ xi f (x(t), t) = 0 and that the symmetric matrix with elements [∂ 2 xi,xj f (x(t), t)] ij is positive semidefinite. For time t + dt we can write: 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), Eq. (3) in the main text, to obtain: Noticing that the first sum in the last line of the previous equation vanishes when evaluated at x(t), and combining the result with Eq. (B1), up to first order in dt we obtain: Thus, whenever the matrix [∂ 2 xi,xj f (x(t), t)] ij can be assumed to be positive definite, we have u(t)− ρ ω ρ (x(t))∆ ρ = 0, in agreement with Eqs. (7) and (8) in the main text.
If the rate function has a single global minimum, then the instantaneous average value x(t) = Ω −1 n(t) = Ω −1 n nP (n, t) also evolves as d t x(t) = u( x(t) ) in the Ω → ∞ limit. In order to see this, we multiply Eq. (1) in the main text by n and sum, obtaining: 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 Eq. (4), then it follows that [37] d t f ss (x(t)) = ρ ω ρ (x) ∆ ρ · ∇f ss (x(t)) ≤ ρ ω ρ (x(t))[e ∆ρ·fss(x) − 1] = 0, where we used Eq. (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.