Nonequilibrium temperature response for stochastic overdamped systems

The thermal response of nonequilibrium systems requires the knowledge of concepts that go beyond entropy production. This is showed for systems obeying overdamped Langevin dynamics, either in steady states or going through a relaxation process. Namely, we derive the linear response to perturbations of the noise intensity, mapping it onto the quadratic response to a constant small force. The latter, displaying divergent terms, is explicitly regularized with a novel path-integral method. The nonequilibrium equivalents of heat capacity and thermal expansion coefficient are two applications of this approach, as we show with numerical examples.


Introduction
The determination of response functions is arguably one of the most topical issues in statistical physics. Even though its history dates back to the works of Einstein, Nyquist and Onsager [1][2][3][4], it was Kubo [5,6] who subsumed the later developments [7][8][9] under a general theory. For a system slightly driven off equilibrium, the Kubo formula gives the linear response of an observable in terms of the equilibrium time-correlation between the observable itself and the entropy produced by the perturbation. The first systematic application of Kubo's theory-along with kinetic theories based on generalised Boltzmann equations-underscored the endeavor to calculate the transport coefficients of moderately dense gasses [10]. These efforts culminated in the discovery of the algebraic decay in time of the correlation functions entering Kubo formulas [11][12][13], which prevents the existence of transport coefficients in low dimensions.
Later, the possibility to perform progressively more efficient computer simulations and thus to compute response functions numerically, led to the extension of the original theory to thermostatted systems arbitrarily perturbed from an initial equilibrium state [14]. Remarkably, it was established that the (nonlinear) response to an external driving is largely insensitive to the choice of the thermostatting mechanisms [15], represented by the artificial forces required to maintain nonequilibrium steady-state conditions [16].
In contrast to such major achievements, the related theory for the response upon perturbation of nonequilibrium states has progressed far more slowly. Apart from the obvious obstacle represented by the lack of knowledge of nonequilibrium phase-space distributions, further difficulties are met when dealing rigorously with deterministic dynamical systems, owing to the fractal nature of their invariant distribution [17][18][19][20]. Nonequilibrium response theories have rather flourished for stochastic dynamics [21][22][23][24][25][26][27][28][29][30][31][32][33][34][35][36], which is applicable to a wide variety of complex systems in physics as well as in related sciences. However, most of these results are usually restricted to mechanical perturbations and do not consider thermal perturbations. Thus, they do not allow one to compute quantities such as nonequilibrium heat capacities and thermal expansions coefficients, which would arise as the (integrated) linear response to step variations of the temperature, i.e., of the noise intensity in the stochastic dynamical equations. Besides some previous formal results [29,37], only recently there appeared formulas for the thermal response of driven stochastic systems, which are given in terms of correlations between state observables calculated in the unperturbed state. Apparently, the mathematical difficulties entailed by handling noise variations require either to introduce an explicit time-discretisation to avoid divergences in the response [38,39] or to rely on a rescaling of the stochastic dynamics in order to derive regular results [40].
The present work is devoted to show that neither of these expedients is actually necessary. A well-defined thermal response formula can be derived by standard path integral techniques, in close analogy to the case of deterministic perturbations. After introducing the model equations in section 2, we define in section 3 the linear response to a temperature perturbation of a generic observable of the system. In section 4 after a brief explanation of the formal differences from the ordinary response to a deterministic forcing, we tackle the problem first showing that the thermal response is equivalent to a portion of the quadratic (i.e. second-order) response to a constant force. Such expression, which displays divergent terms, is then explicitly regularised in section 5 and is showed to be equivalent to a Kubo formula in equilibrium. In section 6 we illustrate two applications of these results: the energy susceptibility of a driven quenched particle (that is the non-equilibrium specific heat for zero driving) and the thermal expansion coefficient of an anharmonic lattice subjected to large heat flows. Moreover, in the simplest tractable case of a freely diffusing particle we connect our formulas to the Einstein relation. A summary and an outlook are finally given in the conclusions.

Overdamped Langevin dynamics
The overdamped diffusive system we consider consists of N degrees of freedom, denoted = { } x x x ,... N 1 . For instance, x j may be a component of a particle position vector in d-dimensions, so that = N nd if the system is composed by n particles. The dynamics is given by the overdamped Langevin equation The jth bath temperature T j and mobility m j (which is the inverse of a damping constant) determine the strength of the noise term, while the drift depends on m j and on the mechanical force ( ( )) x F t j . Such structure respects local detailed balance and thus assumes that the baths are noninteracting with each other and always in equilibrium, regardless of the nonequilibrium conditions experienced by the system. Temperatures and mobilities in our formalism do not depend on the coordinates, hence there is no ambiguity in the interpretation of the stochastic equation. Throughout this paper we will always consider the Stratonovich convention, that is the midpoint rule is employed to discretise in time (1) [41], which means that none of the integrals will be of the Ito type and the rules of standard calculus can be applied. See appendix A for more details.
The F i ʼs are generic nonconservative forces that may bring the system arbitrarily far from equilibrium. In the resulting statistical averages, denoted á¼ñ, there is an understood dependence on the initial density of states r ( ) . This may coincide or not with the steady state density. Finally, we introduce the backward generator of the Markovian dynamics (1), written as a sum of 'one-coordinate' operators j , to avoid clutter. It gives the average time derivative of a state observable ( ) Hereafter for any state observable we use the shorthand notation º to indicate the implicit (and possibly explicit) dependence on the time t.

Linear response in path integral formalism
We imagine to perturb the system (1) varying the noise amplitude through a time dependent parameter q ( ) t 1  switched on at time t=0, namely where  i is a constant determining the ith amplitude of the perturbation. This renders (1) for a perturbed degree of freedom into the form Without loss of generality we assume the mobility to be independent of temperature. The extension to the case where m m = Q ( ) i i i does not involve particular difficulties, since the linear response would be just the sum of the temperature response here described plus a standard response to a deterministic perturbation [34,35], which arises linearising the termm F i i . The aim is to calculate the linear response of a generic observable ( )  t to the just introduced temperature change, defined by Here á ñ q ... denotes an average performed in the perturbed dynamics (5) starting from the state r ( ) x 0 0 , which is unaltered by the perturbation. The associated path weight, proportional to the probability of a trajectory with the action functional The last term in (8) appears as the functional Jacobian in deriving the path-weight for [ ] x from the Gaussian path-weight associated to the noise x i , and depends on the convention used to discretise (5) (e.g. it would be absent with the Ito convention). In the following we will also make use of the unperturbed action º q q= |   0 , which amounts to replacing Q j with T j in (8).
Deep physical insights come from separating any action of the form (8) into time-antisymmetric ( ) and time-symmetric (  , 0 ) components: The integrated entropy flux [ ]  x is the antisymmetric part of the action  under the time-reversal transformation . It is defined consistently with thermodynamics as the sum of the individual heat fluxes into the reservoirs, each weighted by the respective bath temperature [41]. The time-symmetric terms have been studied in connection with the notion of dynamical activity, formerly introduced in the context of jump systems [43][44][45] ). The functionals  and  are then the statistical weights of such selected trajectories. Therefore, in the following we will reserve the name dynamical activity for , which was shown to be a good measure of the system activity [46]. Written as it may be seen as a time-integral of a state variable ( ) x V eff that, for systems with interactions deriving from an energy potential ( ) x U and with a global bath temperature T, would read Such quantity was called effective potential [47,48] and is proportional to the escape rate from a configuration x, as the probability to remain in x for a short time Dt is~- . For our nonequilibrium systems we generalise such concept by writing The escape rate of the degree of freedom x j , denoted l j , follows from evaluating the action at fixed x along a very short trajectory of duration Dt 1  , that is, . In the following sections we will sometimes also use the name frenesy for describing correlation functions in the response formulas involving time-symmetrical features. This alternate naming originated in the responsetheory framework [49] and usually refers to quantities akin to -more specifically, to its excess generated by a perturbing force-namely to quantities assessing the system impatience for changing its state (rather than direct measures of the trajectory zigzags). Hopefully the double terminology is guiding the reader through the connections with the recent literature.

Response to heating as response to a force
We are now in the position to develop the thermal linear response theory, but we immediately find an obstacle.
Since the path weight (7) is normalised to one, ò x P 1, the functional measure q x in (6) contains the noise temperatures Q j (see e.g. [42,50]), and therefore depends itself on the external parameter θ. This is a major difference with respect to an external perturbation of the deterministic forces, which leads to the formal difficulties reported in [38], namely the introduction of an explicit time-mesh to avoid singularities in the results. To overcome this problem we first seek a more manageable expression for the path average. That is obtained through an Hubbard-Stratonovich transformation [51] of the action that, introducing an auxiliary variable y, linearises the quadratic term in (8) and removes the θ dependence from the functional measure of the path weight (see e.g. [50]). By doing so, it is easy to bring (6) is the second-order response function to a constant force perturbation f i of the ith degree of freedom [52], namely , 17 ... f now denotes the average with respect to the perturbed dynamics Formal calculation of response functions to external forces poses no technical difficulty [23,32,34]. After integrating out the auxiliary variable y, it is straightforward to find for (16) Summing up, a standard Hubbard-Stratonovich transformation has allowed us to write the linear response of an observable  to a temperature change as the second-order response to a state-independent force, thus arriving at the intermediate result As anticipated, this result is slightly different from that of a previous approach [38] where the Ito convention was adopted for the path-integrals. Let us add an alternative, intuitive mapping between linear thermal response and quadratic force response, through a less formal derivation of (16). To the purpose, it is sufficient to consider only one degree of freedom. Defining the small parameter m mq and splitting the noise into two independent, zero-mean and white Gaussian noises η and χ, equation (5) reads In view of equation (21), all trajectories can be regarded as generated by the noise η and perturbed by the external random force m c f . The corresponding response is obtained by further averaging over χ. Essentially, we wish to connect the average response to m c f with the response to the deterministic force mf . We thus write the path weight associated to (21) for a single realization of χ, and expand it up to second order in the perturbing force: We recognise the latter as the equal-time ( ) O f 2 term in the path weight associated to (18). So, upon application , the temperature response is also obtained by applying d dq ¢ ( ) t to (23), and we arrive at equality (16).

Regularization of the response
In (20) the divergence caused by the Dirac delta formally compensates the divergence in the squared velocity. This can be heuristically understood recalling that (20), despite being formally expressed in continuous time notation, can be interpreted in terms of discrete, albeit small, time intervals Dt [42,53]. Therefore one has Ḋ x t 1 i 2 , being the dynamics diffusive at short times, and clearly d~D ( ) t 0 1 . However, it would be convenient to recast (20) as an explicit result devoid of singular terms. In the following we perform such operation, first for a single degree of freedom (N = 1), and then extending the result to arbitrary N.

One degree of freedom
With one degree of freedom the parameter  i is superfluous and is thus set to 1. We first focus on the kinetic-like term by starting with the rewriting (valid for > ¢ t t ) 5 This can be achieved recalling that the integral of a total derivative involving the path weight is null. Therefore, we may exploit the identity t t x 0 2 0 corresponding to (8) calculated at q = 0, with N=1. First, we evaluate the second term in (25) making use of the expression for the functional variation of the action derived in appendix C, see(C.3). The entropy variation is shown to vanish, while the variation of [ ]  x expressed in terms of the backward generator  gives Hereafter we restrict to the case in which F does not depend explicitly on time, but only via x. In order to extract from (27) If  is a state observable, i.e., it depends only on the trajectory endpoint, the first term on the right hand side of (28) drops for all ¢ ¹ t t , since it reads Putting all the pieces together we get the compact expression which, plugged in the response formula (20), gives finally This is a regularised version of (20) valid for N=1 and any state observable . We have traded the kinetic-like term and the Dirac delta in (20) with a second-order time derivative and a correlation involving the backward generator. The second-order time derivative, even tough unusual for a linear response formula (but not for a second-order response function [52]), is indeed necessary to obtain the correct result, as it can be easily verified in the analytically solvable case of a particle in free diffusion (see section 6.3).
If one is interested in the response of path-dependent observables (namely,  is a functional of the trajectory up to time t), the first summand in (28) is non-zero and hence (30) has to be supplemented by the term . As an example we may consider the heat exchanged with the thermal bath in a time t, . It turns out that the response formula (30) requires no additional term in this case, since

Many degrees of freedom
The procedure is easily extended to a system composed of > N 1 degrees of freedom. Equations (24),(25) and (28) are still valid replacing x with x i , and taking the action (corresponding to (8) where we reverted to the notation accommodating the particle labels. Equation (27) is then generalised to (see In the following we focus on systems with two-body potential interactions, deferring the more general result (valid for arbitrary d, generic driving and interactions) to appendix C. Yet, the results reported here are general enough to describe the thermal response of a broad class of non-equilibrium systems, such as heat conducting lattices in contact with different heat baths (equation (38)), and aging systems (equation (40)). Under the above assumption, the variation of [ ]  x in(33) is given by where a state observable  was considered. Finally, using the explicit form of the entropy variation (C.10), we find for the response function ( as the total generator of the dynamics in the complete state space. Consequently, for isothermal systems the response formula takes the simpler form ( which is a straightforward generalisation of (30) to a many-body system. As noted above, if  is a path-dependent observable one needs to include in the response formula the additional term coming from the first summand of (25). For the example of the total heat flux into the reservoirs, and thus vanishes when the interactions derive from a two-body potential.

Susceptibility
Upon integration of (38) we get an equation for the susceptibility of the system where we recall that integrals are in the Stratonovich sense and ( )  T i was introduced in (35). The term S 1 is the standard correlation between observable and entropy production, appearing with a 1/2 prefactor with respect to the equilibrium version (see next section), in which it would be the only correlation relevant for determining the linear response. The term S 2 is a novel correlation between observable and a time-antisymmetric quantity, proportional to the functional variation of the bath entropy x , which may be non-zero only if ¹ T T j i for some j. The remaining correlations, the frenetic terms [49] K 1 and K 2 , collect correlations between the observable and time-symmetric dynamical features. As in previous studies of force perturbations, both Sʼs and Kʼs contain, respectively, the entropy and frenesy [49] in excess due to the perturbation.
In order to correctly evaluate the time derivative of the correlation in K 2 , when dealing with data it is important to avoid taking discrete-time derivatives with ¢ > t t because cusps are not unusual in correlation functions for ¢  t t. To 2 is associated with K 2 s .  Kubo 2 and the corresponding susceptibility is

A steady state formula and its reduction to the Kubo formula at equilibrium
is the heat transferred to the system in the time interval [ ] t 0, . This formula shows that the temperature response in equilibrium is totally determined by the correlation between observable and the entropy ( )  t T paid by the reservoir to change the system energy. When a global perturbation is applied to an isothermal steady state regime, say with = "  i 1 i , equation (40) may be recast in an alternative form, that correctly reduces to the Kubo formula (46) in equilibrium, as we show in the following. In the derivation we stay in a generic steady state condition until the very end, so that in turn we obtain another quite general formula for the response function, equation (51) below, in which the genuine nonequilibrium contribution is well distinguished from the Kubo correlation. A possible practical issue of such elegant separation is that it can be computed explicitly only if one knows the microscopic probability density of states.
We start noticing that the last term in(40) is in equilibrium half of the expected result: The remaining frenetic terms yield an analogous contribution at equilibrium. To show that, we first use that the system is in a stationary state. This implies that correlations are functions of the time difference only, hence ¢ s s is the state velocity, that is the probability current J s j associated to x j , over the steady state density of the system r s [49,54]. We will ultimately exploit the time-reversal invariance of equilibrium states, which formally manifests in the equality = *   , as the probability currents v j are by definition absent at equilibrium.
Together with stationarity, we used that *  is the adjoint of , and the equality m in the last passage. We then turn to the second and third summand in (40), starting with the rewriting m Here we introduced the commutator acting as, e.g., Finally, at equilibrium the Kubo formula (46) is retrieved by setting = " v j 0 j and using the rewriting (48) for potential forces. Equation (51) is a thermal response counterpart of previous results for the steady-state force response based on the notion of state velocity [26,27].
6. Examples 6.1. Specific heat for a quenched toy system In this first example we want to highlight that this framework is valid not only for steady states but also for transient regimes. There is to recall an understood dependence of the statistical averages á¼ñ on the initial density of states r 0 .
Let us consider a paradigmatic model of nonequilibrium overdamped systems, namely a single particle in a periodic potential = ( ) U x x cos and subject to an additional constant force f, for simplicity with mobility m = 1.
, in the evolution equation (1) of the unperturbed system. The backward operator acts on the force as sin . To generate a transient condition we choose to thermalise the particle at ¹ T T 0 and to switch to T only at t=0, when the perturbation is also applied. In this way, even for f=0 one cannot apply the Kubo formula for equilibrium systems, as the initial state in not in equilibrium at temperature T. Due to the periodic potential, as an arbitrary procedure for obtaining a well defined r ( ) x 0 , we shift to the interval p [ ] 0, 2 any x obtained from a long simulation run. However, averages such as á ¢ ñ ( ) ( )  x t t 2 need to be computed with x interpreted as a nonperiodic coordinate. We adopted a Heun scheme [41] to integrate the stochastic equation, because it yields trajectories that are consistent with the Stratonovich path-weights used in our theory.
In figure 1 we show examples of susceptibilities of the internal energy ( =  U ) to a change of T for = T 5 0 and T = 0.3, both for f=0 and f = 0.7. We compare the susceptibility c q ( ) t U, from (43) with that computed directly as h T 100 active from t=0 on. We note that, for f=0, the force F is potential and thus the heat exchanged with the bath reduces to an energy difference, Therefore, the susceptibility of the energy gives in the long-time limit the specific heat C of the system: (46) were valid, twice the entropic term (44a) would yield the response. One can note that this is not the case, rather all terms in the response formula are relevant for determining the correct form of the susceptibility.

If a Kubo formula
In these examples, in particular, the term (44d) is especially important. Being the derivative of a correlation function, it is however the noisiest one. One could resort to some high-frequency filtering for better results. In the example of the following subsection we will show that (44e) is a good alternative to (44d) in case one is dealing with steady states.

Thermal expansion in a temperature gradient
In equilibrium at a given temperature T, the correlation function between the heat absorbed by a system and its length may be used to predict the thermal expansion response. In this example we show how this picture breaks down out of equilibrium, where, as exposed in the previous sections, one needs to know also correlations between length and time-symmetric observables, given by (44c) and (44d) or (44e), as well as the new entropic form (44b) due to temperature unbalances. This example specialises to steady state conditions but, with respect to the previous examples, it includes the more general setup of multiple heat baths, in which one can exploit the general formulation with perturbation amplitudes  i . Let us consider the N degrees of freedom arranged in a one-dimensional chain. The system has an energy term is a pinning potential on the first site, and x i ʼs represent the displacements from the average positions. The length of the system in excess with respect to the length at zero temperature, º -X x x N 1 , increases on average for increasing T i ʼs due to the asymmetric two-body potential u(r) (see the inset of figure 2(b)). As a paradigm of nonequilibrium conditions, the system is driven by a set of temperatures varying linearly from T 1 to > T T N 1 . We study the response of the length X to temperature variations, in the form of (a) a global constant increase of the temperatures given by a constant =  obtained with a constant h=0.005 turned on at t=0. From figure 2 one also sees that the entropic and frenetic terms have opposite trends, between each other and with switched roles in the two cases, complementing each other to sum up to the correct response level. In figure 2(b) we also show the response c q X, obtained by an

Free diffusion of one degree of freedom
Let us consider the equations of motion (1) for free diffusion of a single degree of freedom, . The noise prefactor mT 2 comes from assuming the bath to be in equilibrium. In this way the mean square displacement of a free particle in a time t is simply , the response of the mean velocity to a small force is the free-particle mobility μ, and the Einstein relation m = D T between diffusion constant D and mobility is found. One can note that the susceptibility of the observable = ( ) ( )  t x t 2 to a change of T is expected to be mt 2 , hence the corresponding response function is m 2 . We show how our formalism reduces to this result.
For free diffusion all terms in (30) drop but the one involving the second derivative. In this case, the response function can be calculated directly from its definition (6) and one can thus prove analytically that both sides of (30) are equal to the same quantity. As we argued above, the response of the mean square displacement to the perturbation Here we used that the initial condition is independent of the perturbation and noise, thus only the noise autocorrelation contributes. On the other hand, the response formula (30) becomes  give an incorrect result as the system is not in a steady state. It is also trivial to verify (16), namely that this result coincides with the second order response to a state-independent force, giving rise to the dynamics

Conclusions
For overdamped stochastic systems far from equilibrium we have obtained the linear response function of generic state observables to a change in the temperature of the Langevin heat baths. Improving a previous result [38], we need not express the response in terms of a finite time mesh, being all the divergencies appearing in the continuous limit removed, and being all terms in the susceptibility standard integrals or derivatives. This was achieved by deriving a sort of Dyson-Schwinger equation [42], i.e., a relation between unperturbed correlation functions involving an arbitrary observable. This method complements and expands our recent results [40] obtained via a different approach, in which the additional noise stemming from the perturbation was turned into mechanical forces by means of a space rescaling. As in many previous examples, in order to describe a nonequilibrium systems, one needs to know more than just the entropy production. The additional information concerns the knowledge of dynamical quantities that are even under the reversal of the arrow of time (squares of forces, etc). Among them we have recognised the change of the time-integral of the effective potential (i.e., the total escape rate integrated along trajectories) upon variation of the perturbed degree of freedom, d , which complements, perhaps surprisingly, the usual entropy production entering Kubo formula.
For the common scenario of isothermal systems in a steady state, we have also shown how to convert the results in a formula that separates the Kubo term from a nonequilibrium additional correlation that includes the state velocity, see (51). Such version is complementary to the others in the sense that it requires the knowledge of the density of states rather than that of dynamical details.
Future developments of this framework should include multiplicative noise, i.e. those cases where the temperature experienced by the particle depends on their positions.

Appendix A. Stochastic convention for path weights
In this context of temperature response, even if equations have a noise prefactor that does not depend on the system's state x, it turns out that the choice of using Stratonovich path-weights rather than Ito ones is not trivial. As discussed previously [39], by differentiating with respect to temperature one proves a response formula that depends on the choice of the path-weight. One can check that the formulas in this paper are indeed different from those found adopting the Ito convention [38]. The adoption of the Stratonovich convention in the path weight (8) is reflected in the Stratonovich productẋF in (20). If we used an Ito convention in (8), thenẋF in (20) would also be of the Ito type and the same equation would not match the result obtained in the Stratonovich convention.
Ultimately, the path-weight, and thus the corresponding discretisation of (1), have to be chosen consistently with the process that generates the sampled trajectory via (1). By sampled trajectory we mean for example a sequence , 2 , , , of configurations sampled stroboscopically every time step Dt . The Ito convention is by construction suitable for numerical data generated by integration of (1) with The variation of the bath entropy is identically zero, unless F is an explicit function of time Since the dynamical activity is independent ofẋ, its variation is simply the derivative of the escape rate from ¢ ( ) i.e., the mean trajectory satisfies Newton's equation with an effective force mF. In the weak-noise limit T 1  , such trajectory becomes the most probable one, being the minimiser of the action. This expression could be obtained directly by applying the backward generator  to the Langevin equation (1), and using that ξ does not depend on x.

C.2. Many degrees of freedom
For > N 1, thanks to the independency of the different thermal noises, the action (32) is simply the sum of 'single-coordinate' actions: with ( )  j following the structure (26). Nevertheless, its variation is not just equal to (C.3) but in general it will contain additional terms owing to the interactions between different degrees of freedom. One indeed finds modified expressions for the variation of the total entropy flux into the (unperturbed) reservoirs We remark that for systems in d=1 (C.9) does not impose any limitation on the driving, that is, one-body nonconservative forces can be present as well, they simply do not enter in(C.7), which concerns only the interactions between different particles. Instead, in > d 1, different indexes i and j in (C.7) may refer to the coordinates of the same particle, thus (C.7) cannot be simplified to (C.10) in the presence of generic non-conservative forces.
It is worth noting that when the equality ¶ = ¶