An update on nonequilibrium linear response

The unique fluctuation-dissipation theorem for equilibrium stands in contrast with the wide variety of nonequilibrium linear response formulae. Their most traditional approach is"analytic", which, in the absence of detailed balance, introduces the logarithm of the stationary probability density as observable. The theory of dynamical systems offers an alternative with a formula that continues to work when the stationary distribution is not smooth. We show that this method works equally well for stochastic dynamics, and we illustrate it with a numerical example for the perturbation of circadian cycles. A second"probabilistic"approach starts from dynamical ensembles and expands the probability weights on path space. This line suggests new physical questions, as we meet the frenetic contribution to linear response, and the relevance of the change in dynamical activity in the relaxation to a (new) nonequilibrium condition.


Introduction
An open system in contact with a large environment is in stationary equilibrium for a given reduced scale of description when at each moment its condition realizes minimal free energy. The total entropy production is then zero and the evolution is time-reversal invariant. Small disturbances break the stationary behavior and after some time an new equilibrium is then established mainly via dissipative effects. The comparison with the original equilibrium is the domain of linear response: since about 60 years ago [1,2], for classical, quantum, open and closed systems basically the same formula relates the response to a small perturbation with an equilibrium correlation function. This formula (later called the Kubo formula) exactly picks up the physical interpretation in terms of dissipation, hence the name fluctuation-dissipation theorem [3].
Response theory continues making sense for stimuli to nonequilibrium systems, where entropy is being produced already before the perturbation. Over the years various types of linear response formulae have indeed been obtained for nonequilibria, including rather diverse but specific models such as in climatology and for glassy or coarsening dynamics. For better orientation the present paper identifies some common traits between the various approaches to put them in a more unified framework. Suppose, for example, that one wants a numerical code for predicting the response of an out-of-equilibrium system, i.e. without actually perturbing the system. Without a clear picture of the available theoretical results, the choice of the method, if any, would be limited and possibly suboptimal. To answer requests like these, we aim at classifying the features of several 'extended fluctuation-dissipation theorems' for nonequilibrium, highlight briefly their strengths and weaknesses, and the eventual relations between them.
Reviews and classifications of nonequilibrium response already exist. The paper by Marini Bettolo Marconi et al [4] has reviewed response in statistical physics in the light of recent fluctuation relations. Seifert and Speck [5] have introduced a classification of some fluctuation-dissipation theorems into three classes. Chetrite and Gupta [6] present a more mathematical view. This work includes the case of deterministic dynamics and attempts a concise classification of linear response for nonequilibrium. We will conclude that there are three main classes. The first class of formulae can be derived from the Kubo-Agarwal formula, which itself starts from a Dyson expansion of the perturbed Markov semi-group. A second approach originates from the theory of dynamical systems and can be applied when the stationary distribution is not smooth; it gives rise to a numerical algorithm which in fact also applies for certain stochastic dynamics, as we show. The third class is much more probabilistic and treats noise as an important observable in the linear response. The synthesis is there provided by introducing the excess in dynamical activity in a second 'frenetic' contribution to the traditional Kubo fluctuation-dissipation formula. In that sense, the name fluctuation-dissipation theorem (even 'extended') is not fully suitable for nonequilibrium systems, as e.g. their return to stationary nonequilibrium is not uniquely characterized by dissipation [7]. On the mathematical side, our classification is to work either with generators (for Markov evolutions) or to work with weights on path space (mostly still limited to classical dynamical ensembles). Referring to those, we call it the analytic approach (section 2) versus the probabilistic approach (section 3). Remarks complementing those in previous sections, and conclusions, are found in section 4.

Analytic approach
The framework we consider is that of Markov dynamics (no memory) on regions of R d . That includes deterministic dynamical systems, being essentially first order in time. It also includes jump processes as described with master equations, but we use here the language of diffusion because it reduces naturally to that of deterministic dynamical systems in the limit of zero noise, realized here with a state-independent diffusion constant D → 0. For even greater simplicity we choose overdamped diffusions where the velocity fielḋ for standard white noise ξ t [8]. In section 3.2, we also treat (underdamped or inertial) Langevin processes. In all these cases it makes sense to speak about the so-called backward generator L, working on observables (Heisenberg picture). The expectation of an observable A at time t 0 is then given by when at time zero the states are distributed with probability ρ, possibly singular with respect to the reference volume element dx on state space. We can also abbreviate that as A(t) /dt = L A (t) . When no confusion can arise, we continue to write A(t) = A(x t ) for the (most often random) value of the observable at time t 0, and B(0)A(t) for the time-correlation between observable B at time zero and observable A at time t. For the overdamped diffusion (1) the generator is A perturbation changes the drift F → F h ≡ F + h F 1 where, for simplicity, we avoid inserting a time dependence in the small amplitude h, applied at all times t > 0. This leads to a perturbed backward generator The change in expectations at times t with respect to what we had for the unperturbed dynamics follows from (2): To linear order in h, yielding a linear response including a susceptibility χ (t) as the integration of a response function of the form More generally, when applying a time-dependent perturbation h s at time s > 0 we also have We mostly restrict ourselves to the case where ρ is stationary: the response depends only on the time difference τ = t − s and we can write as the central object of study for the linear response of stationary Markov evolutions within an analytic approach 4 . The direct reading of the right-hand side of (9), further discussed in section 2.2, is that e τ L A evolves the observable A for a time τ and then L 1 acts on the result to evaluate it in state x. However, we start in the next section with the more frequent approach of acting on ρ.

Acting on probabilities
In this section, we focus on manipulations with the stationary probability distribution ρ. The basic step from (9) is partial integration, which means that it is assumed here that ρ has a smooth density with respect to the reference volume element, ρ(dx) = ρ(x) dx. In many cases that appears to be a reasonable physical assumption when the level of description is mesoscopic to macroscopic, independent of whether the system is driven or not (see [9]).

The Kubo-Agarwal formula.
Assuming a smooth density ρ we have that (9) can be rewritten as where the adjoint L 1 is defined by dx( Adjoints are forward generators of the time-evolution on densities, as appears e.g. in master equations. For the diffusion (3) the adjoint of L is the operator of the Fokker-Planck equation (∂ t ρ = Lρ, see [8]) and L 1 ρ = −∇ · (F 1 ρ) so that (10) takes the form which is a specific realization of with observable B(x) = L 1 ρ ρ (x). Note that in general the stationary expectation B = 0 because dx L 1 ρ(x) = 0 from the normalization of ρ. Applications of that Agarwal formula [10] in practice meet the difficulty of needing to know the density ρ (which is usually not available) and the details of the dynamics for L 1 . It is thus a result (from partial integration) on a formal level.
Formula (13) is associated with equilibrium, see formula (2.5) in [11]. Of course we have only used that ρ is smooth. It is easy to verify that in the case of detailed balance with conservative forces F = −∇U at temperature D = 1/β we obtain the Kubo formula [2] for the linear response relation in equilibrium. Indeed, say for (3) with perturbation F 1 = −∇V around equilibrium ρ ∝ exp[−βU ], we obtain so that which is the (classical) Kubo formula. The last equality has used that under detailed balance  [13]. Hänggi and Thomas in [14] find the Agarwal formula in their equation (3.12) for time-dependent processes. In the review [4], formulae (2.70)-(2.73) are the Agarwal formula. We also find formulae such as (12) involving log ρ presented in the book by Risken, formulae (7.10)-(7.13) in [8], and formula (7) in [15]. There also re-started the emphasis on log ρ as 'generalized' potential. [16] derived for tiny displacements of the initial condition along a unit vectorê,

Information potential. The formula that Falcioni et al
(see it also as formula (3.13) in the review [4]) corresponds to the case of diffusion with constant perturbation F 1 = −ê (impulsive constant force in the direction −ê) in (12). Again, as there is no explicit dependence on the noise level the formula can be readily tried for chaotic dynamics with smooth stationary density [17]. The function I ≡ − log ρ is sometimes called the information potential. In fact, that potential gets a prominent place in various works on nonequilibrium linear response that follow the Hatano-Sasa formalism [18], such as in the more recent works [19,20].
Other formulae focusing on I can be found in the works by Prost et al [20] and by Speck and Seifert [5]. Consider the stationary density ρ h for the perturbed dynamics with generator L h and the linear response On the other hand, or This formula with special emphasis on the presence of the information potential appears in [5] (but with a different time-derivative), where I is called stochastic entropy. The result of [20] starts from expanding the Hatano-Sasa identity, which effectively makes a special choice for the observable A. There one imagines the process and hence also the I to depend on a family of parameters λ and one asks how the expected value of A = ∂ λ I at λ = λ * changes under a small change in these parameters around a given λ * . The answer is provided by (20), which is equation (5) in [20]. As already noticed for (13) a disadvantage of incorporating the information potential into the correlation function of (16) and (20) is that it is generally unknown and not directly measurable. Not only do we not know the stationary density, but also no clear physically relevant interpretation in terms of heat, work or thermodynamic potential has been found for I. On the plus side, a clear advantage of (16) and (20) is that no detailed information on the dynamics is needed and that one can use parameterized forms of the stationary density ρ to put in I. Typically, (quasi-)Gaussian approximations are tried for ρ to make formulae (12), (16) and (20) more explicit and practically useful; see e.g. [12,[21][22][23]. Note that there is no reference in formula (12) to the noise strength except through the stationary ρ. If one replaces the ρ in the B of (13) by a Gaussian, its variance will effectively reflect the noise level but needs to be fitted.
Finally, recent works on time-dependent processes with feedback rewrite the positivity of the entropy production in terms of the expected information potential (relative entropies), which appears useful for understanding work relations, see e.g. [24,25]. At any rate, (20) is the climax of the analytic approach for smooth probability densities, preserving the Kubo form (15) most faithfully.

State velocity.
There is a relation between the information potential I and the statespace velocity u. For diffusion processes with generator (3) this state velocity is as can be readily checked from the expression for the stationary probability current j ρ appearing in the Smoluchowski equation ∂ t µ t + ∇ j µ t = 0 expressing conservation of probability µ t . Of course this probability velocity need not be related to physical currents, but an interesting observation writes a nonequilibrium response formula as a co-moving (equilibrium) fluctuation-dissipation theorem.

In a stationary process two quantities A and B have time-translationally invariant correlations B(0)A(t) = B(−t)A(0)
. As we mentioned above, the Kubo formula (or the fluctuation-dissipation theorem) holds under detailed balance, i.e. for an unperturbed evolution which gives rise not only to a stationary but also to a time-symmetric distribution of trajectories. Nonequilibrium means breaking time-reversal invariance. We introduce the operator L * that generates the time-reversed motion to describe for example One can alternatively write , from which we see that L * g = L(ρg)/ρ. For the overdamped diffusion processes that we have considered so far, from (3) and (11), L * = L − 2u · ∇ = −F∇ + D + 2D∇logρ∇. Note now that the Kubo formula (15) would follow from the Agarwal formula when in (13) we take B = −β L * V . Indeed, Detailed balance corresponds to ) and corrections to the Kubo formula will thus arise from choosing in (13) for some function V . In other words the antisymmetric part L−L * 2 V in the formula (24) for B will be responsible for violating the Kubo formula and will show the nonequilibrium aspect. Moreover, L − L * = 2u · ∇ relates to the state velocity (21). Substituting these relations via (24) in (13), we find that for diffusions (3) (with D = 1/β) abbreviated as Equation (26) shows that the equilibrium Kubo form gets 'restored' when describing the system in the Lagrangian frame moving with drift velocity u. The passage to the Lagrangian frame of local velocity d/ds → d/ds − u · ∇ 'removes' the nonconservative forcing from the formula, as explained in [15,26,27]. Still, if we do not know the stationary ρ, the formula (26) implies a statistical average over the unknown vector u.

Acting on the observable
We go back to the original (9) and we move to an algorithm focusing on the evolution of the observable A rather than that of the density ρ. Various contributions to the theory of dynamical systems start exactly from that formula for the formulation of linear response as there is still no assumption needed on the smoothness of the probability ρ.
The main player now is L 1 e t L A in (7), and L need not commute with L 1 . Typically, the generator of the perturbation is L 1 = F 1 · ∇ for states x ∈ R d . We are thus interested in obtaining a useful relation for ∇ e t L A. In other words, we need to start from Ruelle treated a formula with this form for deterministic dynamical systems [11,28], recently applied especially in studies of simple models for climate response [22,29,30]. Indeed, (27) is suited in case the stationary probability law is strange in the sense that its lack of smoothness forbids further partial integration to go back to the Agarwal formula (13). In that same context however, it is useful to further split (27) into two parts [11], one describing fluctuations along the stable directions of motion and one parallel to unstable directions (those with positive Lyapunov exponents). Informally, you would expect that the vector ∇ e τ L A (x) can be given, at least for large τ , in its natural components exp[τ λê(x)]ê · ∇ A(x) for local Lyapunov exponents λê corresponding to the various (stable and unstable) directionsê. Such a decomposition is natural for stationary distributions ρ which belong to the class of Sinai-Ruelle-Bowen (SRB) measures; these have a density in the unstable (expanding) direction (positive Lyapunov exponents) so that there a further partial integration remains possible towards an Agarwal formula (13). That is exactly what is suggested in the hybrid form of section 2.3 of [22]. See [31,32] for an introduction to SRB measures with their physical interpretation.

Algorithm for deterministic and stochastic systems.
The numerical evaluation of linear response for chaotic dynamics has recently been studied by various groups. In particular, Abramov and Majda have developed new computational approaches based on the theory of SRB-measures [29]. We simplify here the presentation to introduce an algorithm that works also for stochastic systems. To exploit (27), no knowledge of ρ is required, but we need to get a hold on ∇ x e L A (x) = ∇ x A(x ) x 0 =x where we have emphasized that the differentiation is on an evolved quantity with respect to the initial condition x. To solve the problem that ∇ (state derivative) and L (generating the time-evolution e t L A(x)) need not commute, we unfold the formula a little further. A practical numerical tool for estimating ∇ x A(x τ ) in deterministic dynamics (where x τ is uniquely determined from x 0 = x) was already presented in [33], but the following numerical method is more efficient for steady states.
Consider first an evolution in discrete time n = 0, 1, . . . (but with a parameter that will allow a continuous time limit → 0), for ξ n anything stationary in time, including possible 'noise' depending on the time n, but that does not depend on the state x. This (ξ n ) is considered frozen so that for its given realization we take (28) as a deterministic dynamics. The main point is that We then need to deal with ∇ x A(x n ) where x n depends on the initial state x through (28). By the chain rule, applied recursively, where ∇ A is a 1 × d row vector (easily computed at all times) and Note now that for each time k, the analytical form of the matrix ∇g k = I + ∇ F does not depend on k or on the 'noise' (ξ n ). Hence, The continuous time version G t (x) can be obtained in a suitable limit n → ∞ for G n (x) with time t = n and time step = t/n. In this limit (28) returns the evolutionẋ t = F(x t ) for the deterministic case, while for stochastic equations an additional suitable rescaling to reproduce, e.g. Gaussian noise with ξ t is needed. At any rate, from where each x s depends deterministically on x for frozen (ξ r , 0 r < s), and T indicates a timeordered integral. The formula (32) or its discretization (31) is ready to be inserted into (27) where the · will first average there over the possible 'noise' and then, with ρ, over the initial condition x. That provides the main algorithm for the analytic approach working on the observable instead of on the probability distribution.
A similar expression can be found in the appendix B of [33] for the case of deterministic dynamics. In that work, however, the estimate of (∇ A)(x n )G t (x) was performed with an adjoint scheme, by integrating numerically backward in time the final value (∇ A)(x n ) with the equation (In [33], they wrote the equation for the column vector rather than for the row vector ∇ A, and a transpose of (∇g) was used.) Although that scheme is equivalent to our (31), it requires a CPU time O(n 2 ) for a perturbation active during all n iterations, as opposed to the O(n) matrix multiplications (2.2.1) for estimating matrices G k for k n. Our scheme is exploiting stationarity: the propagator from time n − k to time n is the same as G k from time 0 to k. In transient regimes we would lose this property and we would also need O(n 2 ) operations. A further study of such a numerical algorithm is contained in [22,29]; in particular appendix A of [29] explains the derivation above (without the generalization to stochastic evolutions). For questions on dealing, in a similar context, with the impact of stochastic perturbations, see [34].

Numerical illustration.
To illustrate the numerical scheme (29) for estimating Ruelle's linear response formula (27) in a simple context, we consider a set of equations introduced in biology to describe circadian cycles, that is, the periodicity of biorhythms, for Drosophila [35,36]. The state space has d = 5 dimensions, with states x = (P 0 , P 1 , P 2 , P N , M). The dynamics couples the concentration M of mRNA with those of four types of proteins, written as P 0 , P 1 , P 2 and P N in [36], where one can find the details of these equations. Denoting we have equations of motionẋ = F(x) of the form with parameters as in previous papers, Integration of (2.2.2) was performed with a simple Verlet scheme with a discrete time step = 0.0025. Starting from random concentrations, the model reaches quickly a cyclic regime, as shown in figure 1.
We check the response in the mRNA concentration M to a change of the rate k S → k S (1 + h), starting from the steady state, i.e. a random phase of the cycle. The dynamical perturbation introduced by that change is h F 1 = (h M k S , 0, 0, 0, 0), since a change in k S affects only the evolution equation of variable P 0 , see (35a). The observable A(x) = M has a gradient (∇ A)(x) = (0, 0, 0, 0, 1), which is coupled to the perturbation via the 'propagator' matrix G n . With the setup of section 2.2.1, we can estimate G n by a sequence of matrix multiplications. Given an evolution in discrete time g(x) = x + F(x) , we are interested in matrix ∇g(x) = I + ∇ F(x) or in coordinates (∇g) i j (x) = δ i j + dF i /dx j . The rows of the matrix (∇ F) are where derivatives with respect to P are denoted as P K P = 1 P+K − P (P+K ) 2 and similarly for coinciding with G 0 . Iteratively, G 1 = (∇g)(x 1 ) · G 0 and so on. By sampling many trajectories, in parallel for perturbed and unperturbed dynamics starting from an initial condition from the steady state, the estimate of δ M(τ )/δh 0 has been obtained by setting a small constant h for τ > 0 and by evaluating the time derivative of ( M(τ ) h − M(τ ) )/ h (we used h = 10 −3 ). The response function instead has been calculated with (29) and (27), and overlaps well with the response, as shown in figure 2. Asymptotic stability in this model [36] probably favors a good performance of formula (27). We computed both the perturbed dynamics and the unperturbed fluctuation formula just to show graphically that the results are the same, but of course there is no additional understanding of the process. Obviously, computing the response only with unperturbed simulations would be convenient as long as the convergence of the method is good. We postpone a more detailed study of the efficiency of this and other numerical schemes to a future work.

Probabilistic approach
Linear response makes sense in general only within a statistical theory. That is to say, sensitive dependence on initial and boundary conditions can create strong and relevant effects beyond linear order on the microscopic scale while macroscopic linearity remains valid. For estimating the mobility we do not investigate the microscopic particle's individual motion at specific times and how it changes under an external field, but we ask for a spatio-temporal averaged current and that includes noise. In fact, for stronger microscopic chaoticity we expect the statistical approach to be more relevant. While various physical variables can show chaotic behavior in their timeevolution, their spatial or temporal averages will typically have much smoother behavior, see also [9] and chapter 6.2 in [37] answering the so-called van Kampen objection.

As noise gets important
In coarsening dynamics of low-temperature spin systems, or in spin glasses, a very long transient regime may exist towards equilibrium, i.e. the nonequilibrium is not imposed by external gradients and the dynamic equations are actually equilibrium ones. It is in this context that an extensive response-literature has been produced in recent years. We mention briefly some results from 2003, whose focus was mostly to develop zero-field algorithms, in other words, exploiting fluctuation-response relations to estimate numerically responses from fluctuations in unperturbed dynamics. Works of Chatelain [38], Ricci-Tersenghi [39] and Crisanti and Ritort [40] introduced new schemes for computing the response of the system from correlations in the spins. They were followed by works of Diezeman [41] and Lippiello, Corberi and Zannetti [42][43][44]. The derivations of these results cannot be reduced in our brief scheme, and various forms of response relations have been obtained. However, the one by Lippiello et al, a specific case of (41) below, emerged as having clear physical significance beyond numerical usefulness, because it matched the form of an earlier study by Cugliandolo et al [45], who proposed a response relation for autocorrelations in Langevin equations. With the present notation (1), perturbation h s x and observable A(x) = x, that would read showing an 'asymmetry' term in addition to the form of the Kubo formula. The derivation of (37) is based on the equality x t ξ t = 2Dβ R(t, s); hence it is centered on the presence of noise ξ t . This was new, certainly with respect to derivations presented in previous sections where eventually noise was just an aspect of the dynamics, not fundamental for the derivations. An equation similar to (37) appears in a more recent study for nonequilibrium steady states [46], based on a path-space formulation by Harada and Sasa [47], see their appendix B. The Harada-Sasa approach has really pioneered the path-space approach of the next subsection. The consistency between the approaches of [42,45,47] indeed shows that these belong to a general framework with new significant physical content. In the following section we discuss the path-space formulation embracing these results and we mention the physical interpretation that goes beyond merely dissipative aspects. Noise becomes visible then as dynamical activity, ruling the time-symmetric fluctuations.

Path-space approach
The origin of dynamical ensembles is the projection of a microscopic Hamiltonian dynamics on the dynamics of reduced variables. That Mori-Zwanzig projection [48,49] originates from making a physical partition of the phase space, depending on the physical situation at hand, in which each microstate X is mapped (many-to-one) to a reduced state x(X ). That induces noise in the reduced dynamics, for example on the mesoscopic level of description. It is then natural to consider dynamical ensembles, i.e. probability distributions on path space where a path refers to a trajectory on the mesolevel. Such was already the approach of Onsager and Machlup starting dynamical fluctuation theory in [50]. These trajectories, under certain limiting conditions (e.g. via a weak coupling limit or via adiabatic elimination), can be described via first-order equations, in which case we meet the Markov processes of the previous section 2. Dynamical ensembles can, however, also describe non-Markovian processes describing important memory effects, see e.g. [51] for the application of a response relation in a visco-elastic medium.
Paths are trajectories ω = (x s , 0 s t) in state space, say looking at the states x s in the time interval [0, t]. For the sake of simplicity we characterize the unperturbed (perturbed) process by the probability weight P(ω) [P h (ω)] with respect to some reference dω. The mathematical idea is to turn perturbed expectations A(t) h into unperturbed ones A(t)P h (ω)/P(ω) . Defining the relative action U(ω) = log P(ω)/P h (ω) and splitting U(ω) = [T (ω) − S(ω)]/2 into a time-symmetric T (ω) and a time-antisymmetric S(ω), we obtain where T 1 and S 1 are the linear contributions in h of T and S, respectively, around h = 0 (where U = 0). This formula is general but the point is that in physical systems with local detailed balance [47] we have a good physical understanding of the path functions T and S [53][54][55][56]. Let us first look under global detailed balance. We already see here that since the observable ] divided by temperature. The first term on the right-hand side of (38c) has thus the same form one finds in equilibrium (39), but in general the unperturbed process is in steady nonequilibrium and S(ω) = βh[V (t) − V (0)] is an excess of entropy production with respect to that already generated by the nonequilibrium dynamics. Hence, the derivative δ A(t)S(ω) /δh s = β d ds V (s)A(t) also mimics the equilibrium Kubo formula.
We now turn to the second term on the right-hand side of (38c). Equation (38c) must be true also for a constant A, for which the response is zero. We thus obtain If δ δh s T 1 (ω) equals a state function B(s), then, since d ds V (s) = L V (s) , one deduces that B(s) = β L V (s), arriving at For diffusion processes (1) with generator (3) and temperature D = 1/β we see that indeed is a state function. The same is true for Markov jump processes (for specific derivations of (41) using stochastic calculus, see [55]). The function L V quantifies the time-symmetric volatility of V under the unperturbed dynamics, also named 'frenesy' [53]. It is related to the 'dynamical activity' in discrete systems, where it quantifies the change in escape rates; T (ω) is the excess in dynamical activity over [0, t]. See [54][55][56] for more details and examples. Thus, there is a new contribution next to the first term in (41) (taking one-half of the Kubo formula (15) and referring as usual to dissipation). In particular, its form now clearly deviates from response formulae such as (20). One then wonders in what sense is the decomposition in formulae such as (41) intrinsically natural. We think it is, as that splitting makes the first term time-antisymmetric and the second term symmetric in time s, which corresponds to the dissipative and the reactive parts of the susceptibility, respectively. It remains to be discovered what is the operational meaning of the notion of dynamical activity, and whether the frenetic contribution is on an equal footing with the entropic part in the response. What is certain, however, is that the dynamical activity and the entropy flux merge for systems close to equilibrium, as it should be indeed for recovering the linear response around equilibrium. So the frenetic part of the response only starts to play its independent role from second order around equilibrium. We then expect that thermodynamic forces in that regime pick up non-gradient contributions that are directly related to the frenetic response.
Note that in (41) one computes averages without needing an explicit knowledge of the stationary probability. On the other hand, some knowledge is required on the dynamics sitting in the generator L. It is still not clear how much kinetic information is truly needed and how practical that gets. We summarize the general physical idea in the formula where Entr [0,t] (ω) is the excess in entropy flux over the time period [0, t] and Esc [0,t] (ω) is the excess in dynamical activity due to the perturbation. The physical challenge is to learn to guess or to find that Esc [0,τ ] (ω) from partial information on the dynamics. For the overdamped dynamics of, e.g., (41) it means to consider the expected rate of change of the perturbing potential V under the original dynamics. The response relation out of equilibrium is no longer a fluctuation-dissipation relation but a fluctuation-dissipation-activation relation.
In fact, the formula can now be turned around and from measuring 'violations' of the fluctuation-dissipation relation one obtains information about the active forces [51]. Harada and Sasa [46,47] already expressed the mean rate of dissipation as the deviation from the equilibrium fluctuation-dissipation relation: such an approach was also used to study active forces for a single F 1 -ATPase molecule [52].

The inertial case
We open a separate subsection on the inertial case of Langevin dynamics because its linear response is much less discussed in the literature despite the obvious interest, for example for models of heat conduction. The main point, however, is that the ideas summarized under (3.2) and (42) remain unchanged. For states (q, p) = (q 1 , q 2 , . . . , q n ; p 1 , p 2 , . . . , p n ) ∈ R 2n of positions and momenta we attach standard white noise ξ i t to each 1 i n, with constant strength D i and a friction coefficient γ i to model heat baths at temperature D i /γ i = T i : The forces F i can contain a nonconservative part but are confining when we want a stationary regime where the particles typically reside in a bounded region. We already inserted the perturbation V (q) with small time-dependent amplitude h s for s 0. The linear response is given by formula (17) in [53]: The first sum again corresponds to the dissipative part from the entropy fluxes in the reservoirs at temperatures T i . The remaining sums give the frenetic contribution. Although not recognized yet, formula (44) can still be rewritten in a similar way as done in (41).
in which the generator L for (43) now reads and is the underdamped version of (3). Note however that the dissipative last term D p f does not contribute in (45) and that these formulae do not work in case D i = 0 for an i with ∂ V /∂q i = 0. In that case the best alternative is probably to apply the algorithm of section 2.2. An important application of (44) is in the modification of the Sutherland-Einstein formula relating transport coefficients such as mobility with fluctuation quantities such as the diffusion constant [59].

Nonlinear responses
There is also a growing number of works on higher-order terms around equilibrium. In fact part of the book by Evans and Morriss is devoted to that [60]. Other references include [61], section 10 in [58] or the more recent [4,[62][63][64]. One typical start is the fluctuation symmetry in the distribution of the entropy flux, transient from the reference equilibrium system as also explained in [65] for thermostated dynamics.
As was emphasized in [66], the main point is probably not to be able to write formal expansions and formulae, but to find useful structures and unifying interpretations. We will not deal with that here, except for mentioning one particular general relation between second-order and first-order terms that has been largely unnoticed [67], and an instance of which has appeared as identity (25) in [66].
A equals A up to flipping the sign of the momenta. All terms are explicitly expressed as correlation functions in the equilibrium reference process (indicated with the superscript · 0 ). The derivation of (47) uses linear response around nonequilibrium for non-state functions (i.e. the correlation function in the right-hand side).

Non-state functions
In the case of observables that depend explicitly on time or are functions of the trajectory over several times (such as products A 1 (t 1 )A 2 (t 2 ) · · · A n (t n )), some of the formulae must change. The basic techniques remain however in place. A typical application is how the heat depends on a change in parameters, e.g. for estimating nonequilibrium heat capacities [68]. Heat is not a state function but varies with the trajectory and the applied protocol.

Note on effective temperature
A traditional approach to violation of the equilibrium fluctuation-dissipation theorem is to imagine an effective temperature in the otherwise unchanged Kubo formula [69]. That idea has had most success with mean field type systems, but it has remained more unclear how the effective temperature can provide a consistent and general tool for realistic systems. Still today it serves as a paradigm for interpreting experimental results, see e.g. [70]. One possible approach for the future would be to associate an effective kinetic temperature with the ratio between the frenetic and the entropic contribution in (42), see appendix A in [55].

The transient case
A lot of attention has been devoted to the linear response behavior for relaxational processes.
We mentioned some of that in section 3.1, but there is no way to be complete. For example, the interest in ageing and glassy dynamics has greatly stimulated the search for modified fluctuation-dissipation relations [71]. Here we emphasize that the analytic approach and in particular the methods of section 2.1 lose their simple structure when the unperturbed reference is time-dependent (and not stationary as was assumed). When one is not starting at time zero from the stationary distribution, one needs to take the average in the unperturbed transient regime. Also causality is automatically verified there as the second term in (41) equals the first term when s > t. This unification in expressions for the transient and the steady regime is only natural, as even the stationary regime is physically and ultimately but a very long transient.

General properties of linear response
Some formal properties of linear response go basically unchanged from equilibrium to nonequilibrium contexts. For example, sum rules and Kramers-Kronig relations only depend on the presence of an underlying Hamiltonian dynamics for the total system plus environment or on causality. Of course, specific expressions can differ but there is physically nothing new, see also [72,73]. Nevertheless, there do exist essential differences. We have already alluded to the fact that the name 'fluctuation-dissipation theorem' is no longer so appropriate because of the importance and complementary character of changes in the dynamical activity (and not only in the dissipation). In fact, the word 'fluctuation' also becomes less correct as the time-correlations in the expressions for R(t, s) no longer express a symmetrized time-correlation. To make that point, let us evaluate the equilibrium Kubo formula (15) when the perturbing potential V equals the observable A: which is the variance (or the 'diffusion') of the displacement V (x t ) − V (x 0 ). That is sometimes called a generalized Einstein relation. In fact all linear transport coefficients for equilibrium can be expressed like that, as was understood already in 1960 by Helfand [74]. For nonequilibrium linear response, that relation including its positivity gets violated, see also [7,59]. In the words of the previous subsection, negative effective temperatures become possible in nonequilibrium. The origin lies in the frenetic contribution, e.g. the second term on the right-hand side of (41), which can overrun the first dissipative term. Examples of negative effective temperature in active matter are already known [75][76][77].

Outlook
Response theory is primarily about predicting the reaction of a system in terms of its unperturbed behavior. As we mostly have in mind response in time, that involves temporal correlations. Therefore, dynamical fluctuation theory can be expected to be most prominent.
(See [78] for the distinction between static and dynamic fluctuation theory: statics looks at deviations around the law of large numbers at single times, for example for the average over many copies of the system, while dynamical fluctuations are around the law of large times, for deviations around time-averages.) Around equilibrium, dynamical fluctuations are governed, just like static fluctuations, by the entropy and dissipation functions. That further enables connecting different response coefficients, such as in Onsager reciprocity or via Maxwell relations. Nonequilibrium makes a more drastic difference between static and dynamical fluctuations. In particular, no useful connections between different types of responses have been discovered for nonequilibrium processes. We think of the analogue of relations between compressibilities, heat capacities and conductivities. In our opinion, that challenge will require finding experimental access to quantities like the dynamical activity which is complementing entropic characterizations in the description of nonequilibrium processes.

Conclusions
We have presented a concise guide (with plenty of references) to the multi-faceted world of linear response for systems out of equilibrium. From there we have discovered similarities but also some previously unnoticed missing pieces. One such piece was the extension to stochastic systems of Ruelle's formulation, which we have introduced, together with an efficient algorithm. That also allows for numerical calculations where explicit knowledge of the density of states is not required, as opposed to other (analytic) formulations that we have described. However, the resulting linear response formula contains a correlation function whose physical meaning is not very clear.
In contrast, a rich physical picture emerges from a probabilistic approach based on pathspace weights, where the stationary distribution is also not needed. Besides the quite different mathematical apparatus compared to the more standard analytical approach, the probabilistic way indeed emerges as the one that currently offers more relations with dynamical fluctuation theory: one has to study also how the system correlates with the activity of the perturbing potential, a time-symmetric quantity complementary to the time-antisymmetric fluxes of entropy. The combination fluctuation-dissipation therefore does not suffice to characterize the linear response of nonequilibrium systems.