An update on the nonequilibrium linear response

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]).

2.1.1. The Kubo–Agarwal formula

Assuming a smooth density ρ we have that (9) can be rewritten as

$$\begin{equation} R(\tau) = \int\, {\rm d} x\, \rho(x)\,\frac{{\mathcal L}_1\rho}{\rho}(x)\, \mathrm{e}^{\tau L}A\,(x), \end{equation} \tag{ 10 }$$

where the adjoint ${\mathcal L}_1$ is defined by $\int {\rm d} x ({\mathcal L}_1\rho)(x)\,A(x) = \int {\rm d} x \rho (x)\,(L_1 A)(x)$. 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 ($\partial _t\rho = {\mathcal L} \rho$, see [8])

$$\begin{equation} {\cal L} \rho = -\nabla\cdot(F\rho)+\Delta(D \rho) \end{equation} \tag{ 11 }$$

and ${\cal L}_1 \rho = -\nabla \cdot (F_1\rho)$ so that (10) takes the form

$$\begin{equation} R(\tau) = -\left\langle [\nabla \cdot F_1(0) + F_1(0)\cdot\nabla \log \rho(0)]\,A(\tau) \right\rangle, \end{equation} \tag{ 12 }$$

which is a specific realization of

$$\begin{equation} \fbox {$ \displaystyle R(\tau) = \left\langle \frac{{\mathcal L}_1\rho}{\rho}(0)\,A(\tau) \right\rangle $} = \langle B(0)\,A(\tau)\rangle \end{equation} \tag{ 13 }$$

with observable $B(x) = \frac {{\mathcal L}_1\rho }{\rho }(x)$. Note that in general the stationary expectation 〈B〉 = 0 because $\int {\rm d} x\, {\mathcal L}_1 \rho (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 ${\mathcal 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₁ = −∇V around equilibrium $\rho \propto \exp [-\beta U]$, we obtain

$$\begin{equation} \frac{{\cal L}_1 \rho}{\rho} = -\beta \nabla V\cdot \nabla U + \Delta V = \beta\, LV \end{equation} \tag{ 14 }$$

so that

$$\begin{equation} \left\langle \frac{{\mathcal L}_1\rho}{\rho}(0)\,A(\tau) \right\rangle = \beta\,\langle LV(0)\,A(\tau)\rangle = \beta\,\frac{{\rm d}}{{\rm d} \tau}\langle V(0)\, A(\tau) \rangle, \end{equation}\noindent \tag{ 15 }$$

which is the (classical) Kubo formula. The last equality has used that under detailed balance 〈LV (0) A(τ)〉 = 〈V (0) LA(τ)〉 (see further details in section 2.1.3). In a way, the Agarwal formula (13) repeats Kubo's original derivation while stopping short before specifying ρ.

Also others have re-found the Agarwal formula, such as in theorem 2, formulae (2.22) and (2.23), of chapter 2 in [12]. Weidlich gives a quantum version: his equation (2.17) replaces ${\mathcal L}_1 \rho \rightarrow \frac {\mathrm {i}}{\hbar }[\rho ,H_1]$ with the commutator of the perturbing Hamiltonian H₁ [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.

2.1.2. Information potential

The formula that Falcioni et al [16] derived for tiny displacements of the initial condition along a unit vector $\hat {e}$,

$$\begin{equation} \fbox {$ \displaystyle R(\tau) = \left\langle \hat{e}\cdot\nabla\log \rho (0)\,A(\tau) \right\rangle $} \end{equation} \tag{ 16 }$$

(see it also as formula (3.13) in the review [4]) corresponds to the case of diffusion with constant perturbation $F_1=-\hat {e}$ (impulsive constant force in the direction $-\hat {e}$) 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 ${\mathcal I} \equiv -\log \rho$ 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 ${\mathcal 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 $\rho ^h(x) - \rho (x) = h\,\rho _1(x) + {\mathcal O}(h^2)$. From stationarity ${\mathcal L}^h\rho ^h =0$ we obtain ${\mathcal L}_1\rho + {\mathcal L}\rho _1=0$. For the Agarwal formula (13) we need

$$\begin{equation} B = \frac{{\mathcal L}_1\rho}{\rho} = -\frac{{\mathcal L}\rho_1}{\rho}. \end{equation} \tag{ 17 }$$

On the other hand,

$$\begin{equation} -\partial_h {\mathcal I}^h = \partial_h \log \rho^h = \lim_{h\to 0}\frac{\rho^h - \rho}{h\,\rho^h} = \frac{\rho_1}{\rho} \end{equation} \tag{ 18 }$$

(∂_h is understood in h = 0) so that for inserting in (13), $B= {\mathcal L}(\rho \,\partial _h {\mathcal I}^h) / \rho$. We conclude that

$$\begin{eqnarray} R(\tau) &=& \int {\rm d} x {\mathcal L}(\rho\,\partial_h {\mathcal I}^h)(x)\mathrm{ e}^{\tau L}A (x)\nonumber\\ &=& \int {\rm d} x \,\rho(x)\,\partial_h {\mathcal I}^h(x)\,L\,\mathrm{e}^{\tau L}A (x) \nonumber\\ &=& \frac{{\rm d}}{{\rm d} \tau}\int {\rm d} x \,\rho(x)\,\partial_h {\mathcal I}^h(x)\,\mathrm{e}^{\tau L}A (x) \end{eqnarray} \tag{ 19 }$$

or

$$\begin{equation} \fbox {$ \displaystyle R(\tau) = \frac{{\rm d}}{{\rm d} \tau}\langle \partial_h {{\mathcal I}}^h(0)\,A(\tau)\rangle $}. \end{equation} \tag{ 20 }$$

This formula with special emphasis on the presence of the information potential appears in [5] (but with a different time-derivative), where ${\mathcal 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 $\cal I$ to depend on a family of parameters λ and one asks how the expected value of $A= \partial _\lambda \cal 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 ${\mathcal 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 $\cal {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–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.

2.1.3. State velocity

There is a relation between the information potential $\mathcal {I}$ and the state-space velocity u. For diffusion processes with generator (3) this state velocity is

$$\begin{equation} u \equiv \frac{j_\rho}{\rho} = F - D\nabla \log \rho = F + D\nabla {\mathcal I} \end{equation} \tag{ 21 }$$

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

$$\begin{equation} \left\langle B(-t)A(0) \right\rangle = \int \rho({\rm d} x)\, A(x) \left(\mathrm{e}^{t L^{*}}B\right)(x),\quad t\geqslant 0. \end{equation}\noindent \tag{ 22 }$$

One can alternatively write $\int \rho ({\rm d} x)\,(Lf)(x)\,g(x) = \int \rho ({\rm d} x)\, f(x)\,(L^*g)(x)$, from which we see that $L^{*}g = {\mathcal L}(\rho g)/\rho$. 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,

$$\begin{equation} \langle L^*V(0) A(\tau)\rangle = \langle V(0) L A(\tau)\rangle = \frac{{\rm d}}{{\rm d} \tau}\langle V(0) A(\tau)\rangle. \end{equation}\noindent \tag{ 23 }$$

Detailed balance corresponds to L = L* (i.e. 〈B(−t)A(0)〉 = 〈B(t)A(0)〉) and corrections to the Kubo formula will thus arise from choosing in (13)

$$\begin{equation} B(x) = -\frac{L+L^*}{2}V\kern-1pt(x) = -L^*V(x) - \frac{L-L^*}{2}V\kern-1pt(x), \end{equation}\noindent \tag{ 24 }$$

for some function V . In other words the antisymmetric part $\frac {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/β)

$$\begin{eqnarray} &&\fl R(t,s) = \beta \frac{{\rm d}}{{\rm d} s}\int {\rm d} x\, \rho(x)\,V(x)\, \mathrm{e}^{(t-s)L}A(x) -\beta\int {\rm d} x\, \rho(x)\, u(x) \cdot \nabla V(x)\, \mathrm{e}^{(t-s)L}A(x),\end{eqnarray}\noindent \tag{ 25 }$$

abbreviated as

$$\begin{equation} \fbox {$ \displaystyle R(t,s) = \beta\frac{{\rm d}}{{\rm d} s}\left\langle V(s)A(t) \right\rangle - \beta \left\langle [u(s) \cdot\nabla V(s)] A(t) \right\rangle $}. \end{equation}\noindent \tag{ 26 }$$

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.

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₁ e^tLA in (7), and L need not commute with L₁. Typically, the generator of the perturbation is L₁ = F₁·∇ for states $x\in {\mathbb R}^d$. We are thus interested in obtaining a useful relation for ∇ e^tLA. In other words, we need to start from

$$\begin{equation} \fbox {$ \displaystyle R(\tau) = \int\,\rho({\rm d} x) F_1(x)\cdot \nabla \,\mathrm{e}^{\tau L}A\,(x) $}. \end{equation} \tag{ 27 }$$

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^τLA (x) can be given, at least for large τ, in its natural components $\exp [\tau \lambda _{\hat {e}}(x)]\,\hat {e}\cdot \nabla A(x)$ for local Lyapunov exponents $\lambda _{\hat {e}}$ corresponding to the various (stable and unstable) directions $\hat {e}$. 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.

2.2.1. 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^τLA (x) = ∇_x〈A(x_τ)〉_x0 = 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^tLA(x)) need not commute, we unfold the formula a little further. A practical numerical tool for estimating ∇_xA(x_τ) in deterministic dynamics (where x_τ is uniquely determined from x₀ = 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),

$$\begin{eqnarray} &&x_{n+1}= g_n(x_n) = x_n + v_n(x_n) ,\quad x_0=x\nonumber\\ &&\mbox{with} \ v_n(x)= \epsilon[F(x) + \xi_n], \end{eqnarray} \tag{ 28 }$$

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 ∇_x〈A(x_n)〉_x = 〈∇_xA(x_n)〉_x where the 〈·〉 averages over the 'noise' (ξ_n). We then need to deal with ∇_xA(x_n) where x_n depends on the initial state x through (28). By the chain rule, applied recursively,

$$\begin{eqnarray} \nabla_x A(x_n) &=& (\nabla A)(x_n)\cdot\nabla_x \big(g_{n-1}\circ\cdots\circ g_0\big)\,(x)\nonumber\\ &=& (\nabla A)(x_n)\cdot {\mathbb G}_{n-1}\,(x)\nonumber\\ &=& (\nabla A)(x_n)\,(\nabla g_{n-1})(x_{n-1}) \cdot \nabla_x\big(g_{n-2}\circ\cdots\circ g_0\big)\,(x)\nonumber\\ &=&(\nabla A)(x_n)\,(\nabla g_{n-1})(x_{n-1})\cdot {\mathbb G}_{n-2}(x)=\cdots, \end{eqnarray} \tag{ 29 }$$

where ∇A is a 1 × d row vector (easily computed at all times) and ${\mathbb G}_k(x) \equiv \nabla (g_k\circ g_{k-1}\circ \cdots \circ g_0)\,(x)$ is a d × d matrix obeying the recursive relations

$$\begin{equation} {\mathbb G}_n(x) = \nabla g_{n}(x_n) \cdot {\mathbb G}_{n-1}(x), \end{equation} \tag{ 30a }$$

$$\begin{equation} {\mathbb G}_0(x)= \nabla g_{0}(x). \end{equation} \tag{ 30b }$$

Note now that for each time k, the analytical form of the matrix $\nabla g_k = {\mathbb I} + \epsilon \nabla F$ does not depend on k or on the 'noise' (ξ_n). Hence,

$$\begin{eqnarray} &&\fl \nabla_x A(x_n) = (\nabla A)(x_n)\cdot\left[{\mathbb I} + \epsilon(\nabla F)(x_{n-1})\right]\cdot\left[{\mathbb I} + \epsilon(\nabla F)(x_{n-2})\right]\cdots\left[\vphantom{x_{n-2}}{\mathbb I} + \epsilon(\nabla F)(x)\vphantom{x_{n-2}}\right].\end{eqnarray}\noindent \tag{ 31 }$$

The continuous time version ${\mathbb G}_t(x)$ can be obtained in a suitable limit n → ∞ for ${\mathbb G}_n(x)$ with time t = n  and time step  = t/n. In this limit (28) returns the evolution $\dot {x}_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 ${\mathbb I} + \epsilon (\nabla F)(x) \simeq \mathrm {e}^{\epsilon (\nabla F)(x)}$ we get formally

$$\begin{equation} \nabla_x A(x_\tau) = (\nabla A)(x_\tau)\cdot{\bf{T}}\exp\left[\int_0^\tau \mathrm{d}s (\nabla F)(x_s) \right], \end{equation} \tag{ 32 }$$

where each x_s depends deterministically on x for frozen (ξ_r,0 ⩽ r < s), and T indicates a time-ordered 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 $(\nabla A)(x_n){\mathbb G}_t(x)$ was performed with an adjoint scheme, by integrating numerically backward in time the final value (∇A)(x_n) with the equation

$$\begin{equation} \partial_t (\nabla A)(x_t) + (\nabla A)(x_t)\cdot (\nabla g)(x_t) = 0. \end{equation} \tag{ 33 }$$

(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 ${\mathcal O}(n^2)$ for a perturbation active during all n iterations, as opposed to the ${\mathcal O}(n)$ matrix multiplications (2.2.1) for estimating matrices ${\mathbb G}_k$ for k ⩽ n. Our scheme is exploiting stationarity: the propagator from time n − k to time n is the same as ${\mathbb G}_k$ from time 0 to k. In transient regimes we would lose this property and we would also need ${\mathcal 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].

2.2.2. 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₀, P₁, P₂, P_N, M). The dynamics couples the concentration M of mRNA with those of four types of proteins, written as P₀, P₁, P₂ and P_N in [36], where one can find the details of these equations. Denoting

$$\begin{equation} \left[\frac{P}{K}\right] \equiv \frac{P}{P+K}, \end{equation} \tag{ 34 }$$

we have equations of motion $\dot x=F(x)$ of the form

$$\begin{equation} \frac{{\rm d} P_0}{{\rm d} t} = k_S M - \nu_1 \left[\frac{P_0}{K_1}\right]+ \nu_2 \left[\frac{P_1}{K_2}\right], \end{equation} \tag{ 35a }$$

$$\begin{equation} \frac{{\rm d} P_1}{{\rm d} t} = \nu_1 \left[\frac{P_0}{K_1}\right]- \nu_2 \left[\frac{P_1}{K_2}\right]- \nu_3 \left[\frac{P_1}{K_3}\right] + \nu_4 \left[\frac{P_2}{K_4}\right], \end{equation} \tag{ 35b }$$

$$\begin{equation} \frac{{\rm d} P_2}{{\rm d} t} = \nu_3 \left[\frac{P_1}{K_3}\right] - \nu_4 \left[\frac{P_2}{K_4}\right] + \nu_d \left[\frac{P_2}{K_d}\right] -k_1P_2+k_2P_N, \end{equation} \tag{ 35c }$$

$$\begin{equation} \frac{{\rm d} P_N}{{\rm d} t} = k_1P_2-k_2P_N, \end{equation} \tag{ 35d }$$

$$\begin{equation} \frac{{\rm d} M}{{\rm d} t} = \nu_S \left[\frac{K_I^n}{P_N^n}\right] - \nu_m \left[\frac{M}{K_m}\right], \end{equation} \tag{ 35e }$$

with parameters as in previous papers,

$$\begin{eqnarray*} & &\left(\nu_S = 0.5,\, \nu_m = 0.3,\, \nu_1 = 6,\, \nu_2=3,\, \nu_3=6,\, \nu_4=3\right)\, {\rm nM\,h}^{-1}\\ & &\left(K_m=0.2,\, K_I=2,\, K_1=1.5,\, K_2=2,\, K_3=1.5,\, K_4=2\right)\,{\rm nM}\\ & &\left(k_S=2,\, k_1=2,\, k_2=1\right)\,{\rm h}^{-1},\quad n=4. \end{eqnarray*}$$

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.

Zoom In Zoom Out Reset image size

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₁ = (h M k_S,0,0,0,0), since a change in k_S affects only the evolution equation of variable P₀, 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 ${\mathbb G}_n$. With the setup of section 2.2.1, we can estimate ${\mathbb G}_n$ by a sequence of matrix multiplications. Given an evolution in discrete time g(x) = x + F(x), we are interested in matrix $\nabla g(x) = {\mathbb I} + \nabla F(x) \epsilon$ or in coordinates (∇g)_ij(x) = δ_ij + dF_i/dx_j . The rows of the matrix (∇F) are

$$\begin{equation} \left(-\nu_1\left[\frac{P_0}{K_1}\right]_{P_0} , \nu_2 \left[\frac{P_1}{K_2}\right]_{P_1}, 0, 0, k_S\right), \end{equation} \tag{ 36a }$$

$$\begin{equation} \left(\nu_1\left[\frac{P_0}{K_1}\right]_{P_0} , -\nu_2 \left[\frac{P_1}{K_2}\right]_{P_1}\!\!\!-\nu_3 \left[\frac{P_1}{K_3}\right]_{P_1}, \nu_4 \left[\frac{P_2}{K_4}\right]_{P_2} , 0,0\right), \end{equation} \tag{ 36b }$$

$$\begin{equation} \left(0, \nu_3 \left[\frac{P_1}{K_3}\right]_{P_1}, -\nu_4 \left[\frac{P_2}{K_4}\right]_{P_2}\!\!\!-\nu_d \left[\frac{P_2}{K_d}\right]_{P_2}\!\!\!-k_1, k_2, 0\right), \end{equation} \tag{ 36c }$$

$$\begin{equation} \left(0,0,k_1,-k_2, 0\right),\end{equation} \tag{ 36d }$$

$$\begin{equation} \left(0,0,0,\nu_S\left[\frac{K_I^n}{P_N^n}\right]_{P_N}, -\nu_m\left[\frac{M}{K_m}\right]_{M} \right), \end{equation} \tag{ 36e }$$

where derivatives with respect to P are denoted as $\left[\frac{P}{K}\right]_P = \frac {1}{P+K} - \frac {P}{(P+K)^2}$ and similarly for $\left[\frac {K_I^n}{P_N^n}\right]_{P_N}$. This matrix with x₀ = (P₀,P₁,P₂,P_N,M)₀ yields $(\nabla g)(x_0) = {\mathbb I} + (\nabla F)(x_0)\epsilon$, coinciding with ${\mathbb G}_0$. Iteratively, ${\mathbb G}_1 = (\nabla g)(x_1) \cdot {\mathbb 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₀ 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⁻³). 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.

Zoom In Zoom Out Reset image size

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–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_sx and observable A(x) = x, that would read

$$\begin{equation} \fbox {$ \displaystyle R(t,s) = \frac{1}{2 D} \left\{ \left(\frac{{\rm d}}{{\rm d} s} - \frac{{\rm d}}{{\rm d} t}\right)\langle x_{s} x_{t}\rangle -(\langle F(x_{t}) x_{s}\rangle - \langle F(x_{s}) x_{t}\rangle)\right\} $}, \end{equation} \tag{ 37 }$$

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.

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 ${\mathcal U}(\omega) = \log P(\omega)/P^h(\omega)$ and splitting ${\mathcal U}(\omega) = [{\mathcal T}(\omega)-{\mathcal S}(\omega)]/2$ into a time-symmetric ${\mathcal T}(\omega)$ and a time-antisymmetric ${\mathcal S}(\omega)$, we obtain

$$\begin{equation} \langle A(t)\rangle^h - \langle A(t) \rangle = \int {\rm d} \omega \,P(\omega) A(x_t) \big(\mathrm{ e}^{-{\mathcal U}(\omega)} - 1\big) \end{equation} \tag{ 38a }$$

$$\begin{equation} \hspace{84pt}= \left\langle A(t) \,\big(\mathrm{e}^{-{\mathcal U}(\omega)} - 1\big) \right\rangle \end{equation} \tag{ 38b }$$

$$\begin{equation} \hspace{84pt}= \textstyle\frac{1}{2}\left\langle A(t){\mathcal S}_1(\omega) \right\rangle - \frac{1}{2} \left\langle A(t){\mathcal T}_1(\omega) \right\rangle+{\mathcal O}(h^2), \end{equation} \tag{ 38c }$$

${\mathcal T}_1$

${\mathcal S}_1$

${\mathcal T}$

${\mathcal S}$

${\mathcal U}=0$

${\mathcal T}$

${\mathcal S}$

$$\begin{eqnarray} \langle A(t)\rangle^h - \langle A(t) \rangle &=& \langle A(t) - A(0)\rangle^h \nonumber\\ &=& \frac{1}{2} \left\langle [A(t)-A(0)]{\mathcal S}_1(\omega) \right\rangle \nonumber\\ &=& \left\langle A(t)\,{\mathcal S}_1(\omega) \right\rangle, \end{eqnarray} \tag{ 39 }$$

because $\langle [A(t)-A(0)]){\mathcal T}_1(\omega)\rangle = 0$ and $\langle A(t)\,{\mathcal S}_1(\omega)\rangle = -\langle A(0)\,{\mathcal S}_1(\omega)\rangle$ by time-reversal symmetry of the reference equilibrium process. The result (39) should equal the Kubo formula (15), and indeed it does as will become clear. When the perturbation is from a potential −h_sV (x_s), it is well established that the time-antisymmetric ${\mathcal S}(\omega)={\mathcal S}_1(\omega)$ must be the path-dependent entropy flux into the environment as caused by the perturbation potential V [57, 58]. If the environment is at uniform temperature, for small constant h for times t > 0 the change in entropy is ${\mathcal S}(\omega)=\beta h [V(t)-V(0)]$, namely dissipated energy h[V (t) − V (0)] 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 ${\mathcal S}(\omega)=\beta h [V(t)-V(0)]$ is an excess of entropy production with respect to that already generated by the nonequilibrium dynamics. Hence, the derivative $\delta \langle A(t){\mathcal S}(\omega)\rangle /\delta h_s=\beta \frac{{\rm d}}{{\rm d} s}\langle V(s) A(t)\rangle$ 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

$$\begin{equation} \beta \frac{{\rm d}}{{\rm d} s}\left\langle V(s) \right\rangle = \frac{\delta}{\delta{h_s}} \left\langle {\mathcal T}_1(\omega) \right\rangle = \left\langle \frac{\delta}{\delta{h_s}}{\mathcal T}_1(\omega) \right\rangle. \end{equation} \tag{ 40 }$$

If $\frac {\delta }{\delta {h_s}}{\mathcal T}_1(\omega)$ equals a state function B(s), then, since $\frac {{\rm d}}{{\rm d} s}\langle V(s)\rangle = \langle LV(s)\rangle$, one deduces that B(s) = βLV (s), arriving at

$$\begin{equation} \fbox {$ \displaystyle R(t,s) = \frac{\beta}{2}\frac{{\rm d}}{{\rm d} s} \left\langle V(s)A(t) \right\rangle - \frac{\beta}{2} \left\langle L V(s)A(t) \right\rangle $}. \end{equation} \tag{ 41 }$$

For diffusion processes (1) with generator (3) and temperature D = 1/β we see that indeed LV (x) = (F·∇V)(x) + DΔV (x) is a state function. The same is true for Markov jump processes (for specific derivations of (41) using stochastic calculus, see [55]). The function LV 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; $\cal T(\omega)$ is the excess in dynamical activity over [0,t]. See [54–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

$$\begin{equation} \langle A(t)\rangle^h - \langle A(t)\rangle = \textstyle\frac{1}{2}\langle \mbox{Entr}^{[0,t]}(\omega)\,A(t)\rangle - \frac{1}{2}\langle \mbox{Esc}^{[0,t]}(\omega)\,A(t)\rangle, \end{equation}\noindent \tag{ 42 }$$

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₁-ATPase molecule [52].

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¹,q²,...,qⁿ; $p^1,p^2,\ldots ,p^n)\in {\mathbb R}^{2n}$ of positions and momenta we attach standard white noise ξⁱ_t to each 1 ⩽ i ⩽ n, with constant strength Dⁱ and a friction coefficient γⁱ to model heat baths at temperature Dⁱ/γⁱ = Tⁱ:

$$\begin{eqnarray} \dot{q}^i &=& p^i, \nonumber\\ \dot{p}^i &=& F^i(q) -\gamma^i p^i + h_t \frac{\partial V}{\partial q^i} + \sqrt{2D^i}\,\xi^i_t. \end{eqnarray} \tag{ 43 }$$

The forces Fⁱ 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]:

$$\begin{eqnarray} \fl R(t,s) &=& \sum_{i}\frac{1}{2T^i}\left\langle \frac{\partial V}{\partial q^i}(q_s) \,p^i_s\,A(t) \right\rangle -\sum_i\frac{1}{2D^i} \Bigg\{ \left\langle \frac{\partial V}{\partial q^i}(q_s)\,F^i(q_s)\,A(t) \right\rangle\nonumber\\ \fl& & - \frac{{\rm d}}{{\rm d} s} \left\langle \frac{\partial V}{\partial q^i}(q_s)\, p^i_s\,A(t) \right\rangle + \sum_j \left\langle \frac{\partial^2\,V}{\partial q^j\partial q^i}(q_s)\,p^j_s\,p^i_s\,A(t) \right\rangle \Bigg\}. \end{eqnarray} \tag{ 44 }$$

The first sum again corresponds to the dissipative part from the entropy fluxes in the reservoirs at temperatures Tⁱ. 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). Supposing Dⁱ = γ/β we must replace there $V\rightarrow \dot {V}/\gamma =p\cdot \nabla _q V/\gamma$,

$$\begin{equation} \fbox {$ \displaystyle R(t,s) = \frac{\beta}{2\gamma} \frac{{\rm d}}{{\rm d} s} \left\langle p_s\cdot \nabla_q V(s)A(t) \right\rangle - \frac{\beta}{2\gamma} \left\langle L(p_s\cdot \nabla_q V)(s)\,A(t) \right\rangle $} \end{equation} \tag{ 45 }$$

in which the generator L for (43) now reads

$$\begin{equation} Lf(q,p) = p\cdot\nabla_q f + (F-\gamma\,p)\cdot\nabla_p f + D \Delta_p f \end{equation} \tag{ 46 }$$

and is the underdamped version of (3). Note however that the dissipative last term DΔ_pf does not contribute in (45) and that these formulae do not work in case Dⁱ = 0 for an i with ∂V/∂qⁱ ≠ 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].

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–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].

$$\begin{equation} \left.\frac{\partial^2}{\partial h^2}\langle A(t)\rangle^h\right|_{h=0} = \beta\left.\frac{\partial}{\partial h} \langle [A(t)-\bar{A}(0)][V(t)-V(0)]\rangle^h\right|_{h=0}. \end{equation}\noindent \tag{ 47 }$$

$\bar {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 〈·〉⁰). The derivation of (47) uses linear response around nonequilibrium for non-state functions (i.e. the correlation function in the right-hand side).

In the case of observables that depend explicitly on time or are functions of the trajectory over several times (such as products A₁(t₁)A₂(t₂)···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.

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].

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.

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:

$$\begin{equation} \langle V(t)\rangle^h - \langle V\rangle = - \langle V(0)[V(t)- V(0)]\rangle = \textstyle\frac{1}{2}\langle [V(t)- V(0)]^2\rangle \geqslant 0, \end{equation} \tag{ 48 }$$

which is the variance (or the 'diffusion') of the displacement V (x_t) − V (x₀). 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–77].

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.

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 path-space 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.