Exact Response Theory for Time-Dependent and Stochastic Perturbations

1. Introduction

Linear response theory successfully describes the behaviour of macroscopic systems not too far from thermodynamic equilibrium, obtaining the corresponding transport coefficients solely using the equilibrium correlation functions of the microscopic fluctuating currents, computed with the equilibrium dynamics, cf. Ref. [1,2,3]. While the range of applicability of the linear theory covers most phenomena occurring on the scale of our daily life, extending well beyond that scale, contemporary science and technology concern scales that are much smaller or that involve extremely large driving forces, which often exceed the linear regime. Indeed, anomalous phenomena, for which the linear transport coefficients do not exist or vanish, are common at the nanometric scale. Strong drivings, such as high voltages, produce large currents that alter the physical properties of the system at hand, leading to non-linear effects. Under external drivings, microscopic motions may be impaired to the point of producing localization or, in any case, drastic changes of states, like phase transitions. In many such situations, dynamical correlations persist in time, and give rise to behaviours that are not fully understood.

Within the field of nonequilibrum molecular dynamics [4] and, in particular, following to the discovery of the Fluctuation Relation, response theory has been generalized so that small systems as well as large drivings can be treated. A general exact (non-perturbative) response theory, also known as TTCF (Transient Time Correlation Function) has been derived [4], which has proven very effective. In particular, recently it has allowed to treat hard nonequilibrium problems, at low drivings [5,6], drastically improving the signal to noise ratio. Moreover, this theory has allowed the derivation of a host of relations concerning nonequilibrium systems, cf. e.g. Ref. [7]. The key ingredient of this exact response theory is known as the dissipation function, first explicitly defined by Evans and Searles, as the energy dissipation rate which verifies the fluctuation relation [8,9]. Through that, the quantities till then singularly used in specific instances of nonequilibrium molecular dynamics studies have been unified in a single general concept, that can be used in the analysis of any dynamical system. The original result, called dissipation theorem in Ref. [10], concerned systems subjected to time independent perturbations, such as an external constant field. An extension to time dependent perturbations, that can be expanded as Fourier series has been developed in Ref. [11], enlarging the phase space, to make the overall equations of motion autonomous. This step is necessary, in general, to preserve the notion of time reversibility, which plays a fundamental role in many statistical mechanics calculations, including those associated with linear response, such as the Onsager Reciprocal Relations and the Fluctuation Relations [2,12,13].

In the present paper, we take that approach one step forward, allowing stochastic time dependent perturbations, that can be represented as Karhunen-Loève expansions [14]. In Section 2, we recall the Karhunen-Loève theorem, from which we obtain the particular representation of the stochastic perturbation of our dynamical systems. Section 3 briefly illustrates linear response theory, highlighting the role of the initial perturbation, and then we also illustrate the time independent exact response formalism. Section 4 starts from the equations of motion of the deterministic system of interest, to which we add the time dependent stochastic perturbation discussed in Section 2. The first step is to make our system autonomous by enlarging the phase space, to accommodate a (Hamiltonian) harmonic oscillator, which embeds in itself the time dependence of the vector field. Then, we redefine the quantities of interest in the enlarged phase space and we distinguish three cases: a) the case in which a single realization of the stochastic process is fixed, which yields a deterministic process in which only the initial conditions are random; b) the stochastic case in which the initial condition is fixed and one should only average over the realizations of the stochastic coefficients of the Kahrunen-Loève expansion; c) the stochastic case in which averages are taken over both the initial conditions in phase space and the stochastic coefficients of the process. Finally, we compute the generalized exact response formula for any observable of our system. Section 5 illustrates the theory using simple examples, such as the perturbed harmonic oscillator, and we compare the linear and the exact response results. In Section 6 we summarize our work, discussing our results, and noting how even the standard linear response theory may benefit from the extended phase space approach. Indeed, that allows to properly treat the initial condition of the time dependent perturbation, if one only considers finite evolution times, as required in many modern investigations of fluctuating observables.

2. The Karhunen-Loève theorem

The Karhunen-Loève theorem allows a stochastic process to be represented as an infinite linear combination of orthogonal functions, analogous to a Fourier series, with stochastic rather than deterministic coefficients. If ${\left\{X_t\right\}}_{t\in [a,b]}$ is a centered stochastic process satisfying certain continuity conditions, one can decompose it as a sum of pairwise orthogonal functions multiplied by random coefficients that are pairwise uncorrelated, hence orthogonal in probability space. More precisely, the following holds [14]:

Theorem 1 

(Karhunen-Loève). Let $X_t$, $t\in [a,b]$, be a square-integrable stochastic process with zero-mean, defined on a probability space ($\Theta ,F,P)$, with continuous covariance function $K_X(s,t)$. Letting $e_k$ be an orthonormal basis on ${\mathcal{L}}^2\left([a,b]\right)$ formed by the eigenfunctions of $T_{K_X}$ with respective eigenvalues ${\lambda }_k$, $X_t$ admits the following representation:

$$X_t=\sum _{k=1}^{\infty }{Z}_k{e}_k\left(t\right),$$

(1)

where the convergence is in ${\mathcal{L}}^2$, uniform in t and

$$Z_k={\int }_a^b{X}_t{e}_k\left(t\right)dt.$$

(2)

Furthermore, the random variables $Z_k$ have zero-mean, are uncorrelated and have variance ${\lambda }_k$

$$\langle Z_k\rangle =0,\forall k\in \mathbb{N}and\langle Z_i{Z}_j\rangle ={\delta }_{ij}{\lambda }_j,\forall i,j\in \mathbb{N}.$$

(3)

As an example, consider the particular case of a Wiener process, $w_t$ with $t\in [0,T]$ and covariance function $K_w(t,s)=min(s,t)$. To find the corresponding eigenvalues and eigenvectors, we need to solve the following integral equation:

$${\int }_a^b{K}_w(s,t)e\left(s\right)ds={\int }_0^{T}min(s,t)e\left(s\right)ds={\int }_0^{t}se\left(s\right)ds+t{\int }_t^{T}e\left(s\right)ds=\lambda e\left(t\right),0\le t\le T$$

(4)

Differentiating once with respect to t yields:

$${\int }_t^{T}e\left(s\right)ds=\lambda e^{\prime }\left(t\right)$$

(5)

A second differentiation produces:

$$-e\left(t\right)=\lambda e^{\prime \prime }\left(t\right).$$

(6)

whose general solution takes the form:

$$e\left(t\right)=Asin\left(\frac{t}{\sqrt{\lambda }}\right)+Bcos\left(\frac{t}{\sqrt{\lambda }}\right),$$

(7)

where A and B are determined by the boundary conditions. Setting $t=0$ in the integral equation gives $e\left(0\right)=0$, which implies $B=0$; setting $t=T$ in Eq.(5) yields $e^{\prime }\left(T\right)=0$, whence:

$${\lambda }_k={\left(\frac{T}{(k-\frac{1}{2})\pi }\right)}^2,k\ge 1.$$

(8)

The corresponding eigenfunctions are thus given by:

$$e_k\left(t\right)=Asin\left[(k-\frac{1}{2})\pi \frac{t}{T}\right],k\ge 1$$

(9)

and A is finally determined by the normalization condition:

$${\int }_0^T{e}_k^2\left(t\right)dt=1whichleadstoA=\sqrt{\frac{2}{T}}$$

(10)

Finally, we get following representation of the Wiener process:

$$w_t=\frac{\sqrt{2T}}{\pi }\sum _{k=1}^{\infty }{Z}_k\frac{sin\left[(k-\frac{1}{2})\pi \frac{t}{T}\right]}{(k-\frac{1}{2})}$$

(11)

where ${\left\{Z_k\right\}}_{k=1}^{\infty }$ is a sequence of independent Gaussian random variables with zero mean and unit variance. This representation is valid for $t\in [0,T]$ for any $T>0$ and, as stated in the theorem, convergence is in the ${\mathcal{L}}^2$ norm and uniform in t.

3. Response theory

Linear response theory is suitable to describe the evolution of observables of systems subjected to small perturbations, while the exact response applies regardless of the magnitude of the perturbation. To set our notation, let us recall the main aspects of the time dependent linear theory [1]. Given a Hamiltonian dynamical system, with Hamiltonian $H_0$, consider a Hamiltonian perturbation $H_p(\Gamma ,t)=-\mathcal{F}\left(t\right)A(\Gamma )$, whose small intensity is $\mathcal{F}$, and A has the dimensions of energy. The new perturbed Hamiltonian is given by:

$$H(\Gamma ,t)=H_0(\Gamma )+H_p(\Gamma ,t)=H_0(\Gamma )-\mathcal{F}\left(t\right)A(\Gamma )$$

(12)

and the equations of motion take the form:

$$\begin{array}{ccc}& & \frac{d{q}_j}{dt}=\frac{\partial H_0}{\partial p_j}-\mathcal{F}\left(t\right)\frac{\partial A}{\partial p_j}\equiv \frac{\partial H_0}{\partial p_j}-\mathcal{F}\left(t\right)k_j^q,with{k}_j^q=\frac{\partial A}{\partial p_j}\hfill \end{array}$$

(13)

$$\begin{array}{ccc}& & \frac{d{p}_j}{dt}=-\frac{\partial H_0}{\partial q_j}+\mathcal{F}\left(t\right)\frac{\partial A}{\partial q_j}\equiv -\frac{\partial H_0}{\partial q_j}-\mathcal{F}\left(t\right)k_j^p,with{k}_j^p=-\frac{\partial A}{\partial q_j}\hfill \end{array}$$

(14)

Calling G the corresponding vector field, we have:

$$G(\Gamma ,t)=\left(\begin{array}{cc}& \frac{\partial H_0}{\partial p_1}\\ & ⋮\\ & -\frac{\partial H_0}{\partial q_N}\end{array}\right)-\mathcal{F}\left(t\right)\left(\begin{array}{cc}& k_1^q\\ & ⋮\\ & k_N^p\end{array}\right)=G_0(\Gamma )+G_{ext}(\Gamma ,t)$$

(15)

and the Liouville operator

$$\begin{array}{cc}\hfill -i\mathcal{L}& =-\left(G·{\nabla }_{\Gamma }+{\mathrm{div}}_{\Gamma }G\right)=-\left(G_0·{\nabla }_{\Gamma }+{\mathrm{div}}_{\Gamma }{G}_0\right)-\left(G_{ext}·{\nabla }_{\Gamma }+{\mathrm{div}}_{\Gamma }{G}_{ext}\right)\hfill \\ \hfill & =-i\left({\mathcal{L}}_0+{\mathcal{L}}_{ext}\right).\hfill \end{array}$$

(16)

allows us to write the continuity equation for the probability distribution in phase space as:

$${\partial }_{t}f=-i\mathcal{L}f=-i\left({\mathcal{L}}_0+{\mathcal{L}}_{ext}\left(t\right)\right)f$$

(17)

where

$$i{\mathcal{L}}_{0}f=\{f,H_0\}=\sum _{j=1}^N\left(\frac{\partial H_0}{\partial p_j}\frac{\partial }{\partial q_j}-\frac{\partial H_0}{\partial q_j}\frac{\partial }{\partial p_j}\right)f$$

(18)

and

$$i{\mathcal{L}}_{ext}f=-\mathcal{F}\left(t\right)\{f,A\}=-\mathcal{F}\left(t\right)\sum _{j=1}^N\left(k_j^q\frac{\partial }{\partial q_j}+k_j^p\frac{\partial }{\partial p_j}\right)f$$

(19)

whose solution can be expressed as:

$$f_t(\Gamma )=e^{-it{\mathcal{L}}_0}{f}_0(\Gamma )-i{\int }_0^{t}d{t}^{\prime }{e}^{-i(t-t^{\prime }){\mathcal{L}}_0}{\mathcal{L}}_{ext}\left(t^{\prime }\right)f_{t^{\prime }}(\Gamma )$$

(20)

As well known, this is an exact, but not useful, expression. However, the fact that ${\mathcal{L}}_{ext}$ is proportional to the small intensity of the perturbation $\mathcal{F}$, justifies the following linear approximation of $f_t$:

$$f_t(\Gamma )=f_{eq}(\Gamma )-i{\int }_0^{t}d{t}^{\prime }{e}^{-i(t-t^{\prime }){\mathcal{L}}_0}{\mathcal{L}}_{ext}\left(t^{\prime }\right)f_{eq}(\Gamma )+H.O.$$

(21)

where $H.O.$ stands for negligible higher order terms in $\mathcal{F}$. More than the probabilities, though, we are interested in the evolution of observables, which are identified with their evolving ensemble averages. Denote one observable by $\mathcal{O}:\mathcal{M}\to \mathbb{R}$, and its ensemble average at time t by:

$${\langle \mathcal{O}\rangle }_t\equiv {\int }_{\mathcal{M}}\mathcal{O}(\Gamma )f_t(\Gamma )d\Gamma$$

(22)

In the linear regime, we can then write:

$${\langle \mathcal{O}\rangle }_t-{\langle \mathcal{O}\rangle }_0={\int }_0^{t}d{t}^{\prime }R(t-t^{\prime })\mathcal{F}\left(t^{\prime }\right)$$

(23)

where

$$R\left(t\right)=\beta {\int }_{\mathcal{M}}d\Gamma f_0(\Gamma )\dot{A}(\Gamma )e^{it{\mathcal{L}}_0}\mathcal{O}(\Gamma )=\beta {〈\dot{A}(\mathcal{O}\circ S_0^t)〉}_0$$

(24)

is called response function and the subscript 0 in $S_0^t$ indicates the unperturbed evolution. This is a very important result: the response to a small perturbation is determined by the equilibrium dynamics and by the time correlation function of perturbation and observable of interest computed with respect to the known equilibrium ensemble. For thermodynamic systems, this description is fully satisfactory, because the observables of interest in such systems only negligibly fluctuate, and the observed signal concerning each single object practically equals the ensemble averaged signal. This is not more guaranteed when non-thermodynamic systems, or large time dependent perturbations are considered. Therefore, an extension of the linear response theory to small systems or large perturbations is required to address these two issues.

The (time independent) exact response theory initially developed in Ref. [10], may be adapted to describe a generic dynamical system, $S^t:\mathcal{M}\to \mathcal{M}$, on a phase space $\mathcal{M}$, where the time t may be continuous or discrete [15]. Let us focus on the case in which $S^t\Gamma$, with $\Gamma \in \mathcal{M}$ represents the solution at time t of the differential equation

$$\dot{\Gamma }=G(\Gamma )withinitialcondition\Gamma \in \mathcal{M}$$

(25)

For any such system, the phase space variation rate, $\Lambda :\mathcal{M}\to \mathbb{R}$, is the divergence in $\mathcal{M}$ of the vector field:

$$\Lambda ={\nabla }_{\Gamma }·G$$

(26)

We denote its time integral along the phase space trajectory segment delimited by the points reached at time s and time t by:

$${\Lambda }_{s,t}(\Gamma )={\int }_s^t\Lambda (S^r\Gamma )dr$$

(27)

which gives the variation factor of the phase space volume element in going from $S^s\Gamma$ to $S^t\Gamma$.

Assume $\mathcal{M}$ is endowed with a probability measure $d\mu (\Gamma )=f_0(\Gamma )d\Gamma$, of density $f_0$, which evolves as prescribed by the Liouville equation

$$\frac{\partial f_t}{\partial t}(\Gamma )=-G(\Gamma )·{\nabla }_{\Gamma }{f}_t(\Gamma )-f_t(\Gamma ){\nabla }_{\Gamma }·G(\Gamma )$$

(28)

and turns $f_t$ after a time t. This can be rewritten in terms of the dissipation function [16]:

$${\Omega }^{f_t}(\Gamma )=-G(\Gamma )·{\nabla }_{\Gamma }ln{f}_t(\Gamma )-\Lambda (\Gamma )$$

(29)

as

$$\frac{\partial f_t}{\partial t}(\Gamma )={\Omega }^{f_t}(\Gamma )f_t(\Gamma )$$

(30)

Starting at a point $\Gamma \in \mathcal{M}$, the corresponding integral of ${\Omega }^{f_0}$ along the trajectory segment delimited by $S^s\Gamma$ and $S^t\Gamma$ is given by:

$${\Omega }_{s,t}^{f_0}(\Gamma )={\int }_s^t{\Omega }^{f_0}(S^s\Gamma )ds=ln\frac{f_0(S^s\Gamma )}{f_0(S^t\Gamma )}-{\int }_s^t\Lambda (S^s\Gamma )ds,$$

(31)

and the solution of the Liouville equation can be equivalently expressed as:

$$f_{s+t}(S^t\Gamma )=exp\{-{\Lambda }_{0,t}(\Gamma )\}{f}_s(\Gamma ),$$

(32)

and as:

$$f_{s+t}(\Gamma )=exp\left\{{\Omega }_{-t,0}^{f_s}(\Gamma )\right\}{f}_s(\Gamma )$$

(33)

Note that the dynamics is assumed to be invertible, which does not mean that it has to be time-reversal invariant; it suffices that each point in a trajectory has a unique preimage under the dynamics. The ensemble average of an observable $\mathcal{O}$ at time t,

$${\langle \mathcal{O}\rangle }_t=\int \mathcal{O}(\Gamma )f_t(\Gamma )d\Gamma$$

(34)

can eventually be expressed as:

$${\langle \mathcal{O}\rangle }_t={\langle \mathcal{O}\rangle }_0+{\int }_0^t{〈{\Omega }^{f_0}(\mathcal{O}\circ S^s)〉}_{0}ds$$

(35)

which is similar to Eq.(23) with Eq.(24). For properly chosen $f_0$, ${\Omega }^{f_0}$ represents the dissipative flux, like $\dot{A}=\left\{A,H\right\}=\left\{A,H_0\right\}$ which appears in Eq.(24). In particular, nonequilibrium molecular dynamics models require $f_0$ to be the equilibrium distribution of the unperturbed dynamics (drivings set to 0) subjected to the same constraints of the perturbed dynamics. For instance, if the perturbed dynamics preserves the kinetic energy thanks to iso-kinetic thermostats, the unperturbed dynamics must also preserve the kinetic energy, and $f_0$ must be invariant under the resulting (generally non-Hamiltonian) dynamics. In this case, analogously to the linear formula (23), Eq.(35) concerns the correlation function of the evolving observable of interest with the dissipative flux and such correlation function is an average computed with respect to the initial distribution $f_0$. The difference is that the dynamics followed by the observable is the exact perturbed dynamics and not the approximate equilibrium one. However, the formalism can now be used quite more generally, without need for $f_0$ to be the equilibrium distribution, or for the system of interest to be a particle system.

4. Exact response of stochastic processes

In this section, we first adapt the exact response theory illustrated above to a generic abstract dynamical system perturbed by a stochastic, time dependent, vector field, extending the approach for deterministic time dependent perturbations developed in Ref. [11]. Then, we distinguish three cases: a) the deterministic case in which a single realization of the perturbation is fixed (as in e.g. the Green-Kubo linear theory), and so observables are only averaged with respect to the initial conditions of phase space trajectories; b) the stochastic case in which the initial condition in phase space is fixed, the realizations of the perturbation vary, hence observables are averaged over the distribution of such realizations; c) the stochastic case in which both phase space initial conditions and perturbation realizations vary, and observables are averaged over both.

4.1. Wiener process perturbation

Given a dynamical system $\dot{\Gamma }=G(\Gamma )$ on a phase space $\mathcal{M}$, consider the following stochastic equation

$$\dot{\Gamma }=\widehat{G}(\Gamma ,t)=G(\Gamma )+\mathcal{F}\widehat{w}\left(t\right),fort\in [0,T]$$

(36)

where $\mathcal{F}$ is a constant that gives the strength of the perturbation, $T>0$ can be much larger than the scale of observation as in [1], and $w\left(t\right)$ is a Wiener process of the same dimension as $\mathcal{M}$. This can be represented by a Karhunen-Loève expansion as follows [14]:

$$\widehat{w}\left(t\right)=\sum _{n=1}^{\infty }\frac{\sqrt{2T}}{\pi }\frac{Z_n}{(n-1/2)}sin\left(\frac{(n-1/2)\pi }{T}t\right)=\sum _{n=-\infty }^{\infty }{\alpha }_n{e}^{i{\beta }_{n}t}$$

(37)

where the details of the process are determined by the vectors $Z_n=(Z_{1n},...,Z_{Nn})$, or equivalently ${\alpha }_n=({\alpha }_{1n},...,{\alpha }_{Nn})$, with:

$${\alpha }_{kn}=\left\{\begin{array}{cc}\hfill \frac{1}{2i}\frac{\sqrt{2T}}{\pi }\frac{Z_{kn}}{c_n}& ifn>0\hfill \\ \hfill \frac{1}{2i}\frac{\sqrt{2T}}{\pi }\frac{Z_{-kn}}{c_n}& ifn<0\hfill \end{array}\right.c_n=\left\{\begin{array}{cc}\hfill & n-\frac{1}{2}ifn>0,\hfill \\ \hfill & n+\frac{1}{2}ifn<0\hfill \end{array}\right.and{\beta }_n=\frac{c_n\pi }{T}$$

(38)

We recall that this expression holds for all perturbations in ${\mathcal{L}}^2[0,T]$, for any $T>0$, so it is a very general result. As this a zero mean process, it does not include a net force acting on the system. However, such a force, constituting a systematic action about which the Wiener process fluctuates, can be included as a deterministic term in G.

Rather than directly computing the evolution of probability densities, we follow Ref. [11] in order to obtain a more flexible tool, and to preserve the notion of time reversal invariance. Therefore, we first eliminate the time dependence of the perturbation, introducing two new variables $\theta$ and $\varphi$ that evolve as follows:

$$\left\{\begin{array}{c}\dot{\theta }=\varphi \hfill \\ \dot{\varphi }=-\theta \hfill \end{array}\right.$$

(39)

These new variables live on an ellipse, whose axes depend on the initial independently chosen conditions ${\theta }_0$ and ${\varphi }_0$. Then, we can write:

$$\mathcal{F}\widehat{w}\left(t\right)=\mathcal{F}\sum _{n=-\infty }^{\infty }{\alpha }_n{e}^{i{\beta }_{n}t}=\sum _{n=-\infty }^{\infty }{\alpha }_n{w}_n(\theta ,\varphi )=w(\theta ,\varphi )$$

(40)

where the deterministic parts of the perturbation are given by the functions $w_n(\theta ,\varphi )$, and the stochastic parts by the coefficients ${\alpha }_n$. When $\mathcal{F}$ is given, we can write:

$$w_n(\theta ,\varphi )={\mathcal{F}}^{1-{\beta }_n}{\left(\theta \left(t\right)-i\varphi \left(t\right)\right)}^{{\beta }_n}withn\in \mathbb{Z},and\theta \left(t\right)-i\varphi \left(t\right)=\mathcal{F}exp\left(it\right)$$

(41)

Alternatively, the initial condition determines $\mathcal{F}$ through the equality $\theta \left(t\right)-i\varphi \left(t\right)=\mathcal{F}exp\left(it\right)$. The difference between the two situations depends on the initial distribution of the points $(\theta ,\varphi )$, which we are now free to choose as the problem of interest requires.

In both cases, the explicit time dependence of the perturbation w has been turned into an implicit dependence, mediated by $\theta$ and $\varphi$. The new dynamical system writes:

$$\dot{\tilde{\Gamma }}=\tilde{G}(\Gamma ,\theta ,\varphi )=\left(\begin{array}{c}G(\Gamma )+w(\theta ,\varphi )\\ \varphi \\ -\theta \end{array}\right)$$

(42)

which is an autonomous system of differential equations, whose phases $\tilde{\Gamma }=(\Gamma ,\theta ,\varphi )$ live in the new phase space

$$\tilde{\mathcal{M}}={\mathcal{M}}_{\Gamma }\times {\mathcal{M}}_{\theta \varphi }$$

(43)

where ${\mathcal{M}}_{\Gamma }$ coincides with the original phase space of the unperturbed dynamics, $\mathcal{M}$, while ${\mathcal{M}}_{\theta \varphi }\subset {\mathbb{R}}^2$ is the space of the points $(\theta ,\varphi )$, which depends on the problem at hand, and may be bounded or not. In particular, a given initial condition confines $(\theta ,\varphi )$ to a circle of given radius $\mathcal{F}$, as for a harmonic oscillator in one dimension, whose energy is fixed by the initial conditions. The corresponding new phase space variation rate equals the old one:

$$\tilde{\Lambda }=\tilde{\Lambda }(\Gamma ,\theta ,\varphi )=\Lambda (\Gamma )$$

(44)

because the new coordinates describe a Hamiltonian system. Introducing a new probability density, $\widetilde{f_0}=\widetilde{f_0}(\Gamma ,\theta ,\varphi )$, a new dissipation function, ${\tilde{\Omega }}^{f_0}={\tilde{\Omega }}^{f_0}(\Gamma ,\theta ,\varphi )$, can also be obtained, applying the general rule (29) to the quantities having a tilde. The density $\widetilde{f_0}$ can be given in various guises. However, the distribution of $\theta$ and $\varphi$ is not affected by the phase space coordinates $\Gamma$, and should not affect the distribution of $\Gamma$, therefore it can be factorized as:

$${\tilde{f}}_0(\Gamma ,\theta ,\varphi )=f_0(\Gamma )g_0(\theta ,\varphi )$$

(45)

In the case the realization of the perturbation is fixed, and its initial value is given by $w_0=w({\theta }_0,{\varphi }_0)$, one could take:

$${\tilde{f}}_0(\Gamma ,\theta ,\varphi )=f_0(\Gamma )\delta (\theta -{\theta }_0)\delta (\varphi -{\varphi }_0)$$

(46)

where $f_0$ is invariant under the unperturbed dynamics, as in the usual molecular dynamics applications of the exact response formula, and $\delta$ is the Dirac delta function. On the other hand, finite accuracy in setting the initial value of w may be described replacing the $\delta$-function with a smooth non-negative and normalized function, centered on ${\theta }_0$ and ${\varphi }_0$, that vanishes in a suitably narrow or wide interval about 0. Physically, this is more meaningful, and mathematically it is useful since it preserves the differentiability of ${\tilde{f}}_0$. In this case, we would write:

$${\tilde{f}}_0(\Gamma ,\theta ,\varphi )=f_0(\Gamma )g_0(\theta -{\theta }_0)g_0(\varphi -{\varphi }_0)$$

(47)

where the relation to ${\theta }_0$ and ${\varphi }_0$ has been stressed. If, on the other hand, the magnitude of the perturbation is fixed, one may take ${\mathcal{M}}_{\theta \varphi }={\mathcal{S}}^{\mathcal{F}}$, the circle of radius $\mathcal{F}$. Then, because this circle is traversed with uniform speed by the point $\left(\theta \right(t),\varphi (t\left)\right)$ and the origin of times is arbitrary, a natural form of the phase space distribution is:

$$\widetilde{f_0}(\Gamma ,\theta ,\varphi )=\frac{1}{2\pi \mathcal{F}}{f}_0(\Gamma )$$

(48)

Let us rewrite the equations of motion as:

$$\dot{\tilde{\Gamma }}(\Gamma ,\theta ,\varphi )=\left(\begin{array}{c}{\dot{\Gamma }}_1\\ {\dot{\Gamma }}_2\\ ⋮\\ {\dot{\Gamma }}_N\\ \dot{\theta }\\ \dot{\varphi }\end{array}\right)=\left(\begin{array}{c}{G}_1(\Gamma )+\sum _{n=-\infty }^{\infty }{\alpha }_{n1}{w}_n(\theta ,\varphi )\hfill \\ G_2(\Gamma )+\sum _{n=-\infty }^{\infty }{\alpha }_{n2}{w}_n(\theta ,\varphi )\hfill \\ ⋮\hfill \\ G_N(\Gamma )+\sum _{n=-\infty }^{\infty }{\alpha }_{nN}{w}_n(\theta ,\varphi )\hfill \\ \varphi \hfill \\ -\theta \hfill \end{array}\right)$$

(49)

where $({\Gamma }_1,...,{\Gamma }_N)=\Gamma \in \mathcal{M}$, and $(G_1,...,G_N)=G$. Given the initial distribution, the new dissipation function is derived as follows:

$$\begin{array}{ccc}& & {\tilde{\Omega }}_{{\tilde{f}}_0}(\Gamma ,\theta ,\varphi )=-\tilde{\Lambda }(\Gamma ,\theta ,\varphi )-\frac{d}{dt}ln\widetilde{f_0}(\Gamma ,\theta ,\varphi )\hfill \\ & & =-\Lambda (\Gamma )-\frac{1}{f_0(\Gamma )g_0(\theta ,\varphi )}[g_0(\theta ,\varphi )\left(G(\Gamma )+w(\theta ,\varphi )\right)·{\nabla }_{\Gamma }{f}_0(\Gamma )+\hfill \\ & & +f_0(\Gamma )\left(\varphi \frac{\partial g_0}{\partial \theta }(\theta ,\varphi )-\theta \frac{\partial g_0}{\partial \varphi }(\theta ,\varphi )\right)]\hfill \\ & & =-\Lambda (\Gamma )-\frac{1}{f_0(\Gamma )}\sum _{k=1}^N\left[G_k(\Gamma )+\sum _{n=-\infty }^{\infty }{\alpha }_{nk}{w}_n(\theta ,\varphi )\right]\frac{\partial f_0}{\partial {\Gamma }_k}(\Gamma )\hfill \\ & & -\frac{1}{g_0(\theta ,\varphi )}\left(\varphi \frac{\partial g_0}{\partial \theta }(\theta ,\varphi )-\theta \frac{\partial g_0}{\partial \varphi }(\theta ,\varphi )\right)\hfill \\ & & ={\Omega }_{f_0}(\Gamma )-\frac{1}{f_0(\Gamma )}\sum _{n=-\infty }^{\infty }{w}_n(\theta ,\varphi )\sum _{k=1}^N{\alpha }_{nk}\frac{\partial f_0}{\partial {\Gamma }_k}(\Gamma )-\frac{1}{g_0(\theta ,\varphi )}\left(\varphi \frac{\partial g_0}{\partial \theta }(\theta ,\varphi )-\theta \frac{\partial g_0}{\partial \varphi }(\theta ,\varphi )\right)\hfill \end{array}$$

(50)

Consequently, the correlation of ${\tilde{\Omega }}_{{\tilde{f}}_0}$ with the time evolution of one observable $\mathcal{O}$ which, by definition, does not depend on $\theta$ and $\varphi$, can be written as:

$$\begin{array}{ccc}\hfill {〈{\tilde{\Omega }}_{{\tilde{f}}_0}\left(\mathcal{O}\circ {\tilde{S}}^s\right)〉}_{{\tilde{f}}_0}& =& {\int }_{\mathcal{M}\times {\mathcal{M}}_{\theta \varphi }}{\tilde{\Omega }}_{{\tilde{f}}_0}(\Gamma ,\theta ,\varphi )\mathcal{O}\left({\tilde{S}}^s\Gamma \right)f_0(\Gamma )g_0(\theta ,\varphi )d\Gamma d\theta d\varphi \hfill \end{array}$$

(51)

$$\begin{array}{ccc}& =& {\int }_{{\mathcal{M}}_{\theta \varphi }}d\theta d\varphi g_0(\theta ,\varphi ){\int }_{\mathcal{M}}d\Gamma {\tilde{\Omega }}_{{\tilde{f}}_0}(\Gamma ,\theta ,\varphi )\mathcal{O}\left(S^s\Gamma \right)f_0(\Gamma )\hfill \end{array}$$

(52)

because the evolution of $\Gamma$ in the enlarged phase space, ${\tilde{S}}^s$, equals the one in $\mathcal{M}$. Let us consider first the third addend of ${\tilde{\Omega }}_{{\tilde{f}}_0}$. Since it does not depend on $\Gamma$, it yields:

$$\begin{array}{ccc}& & -{\int }_{{\mathcal{M}}_{\theta \varphi }}d\theta d\varphi g_0(\theta ,\varphi )\frac{1}{g_0(\theta ,\varphi )}\left(\varphi \frac{\partial g_0}{\partial \theta }(\theta ,\varphi )-\theta \frac{\partial g_0}{\partial \varphi }(\theta ,\varphi )\right){\int }_{\mathcal{M}}d\Gamma \mathcal{O}\left(S^s\Gamma \right)f_0(\Gamma )\hfill \end{array}$$

(53)

$$\begin{array}{ccc}& & =-{〈\mathcal{O}\circ S^s\Gamma 〉}_0{\int }_{{\mathcal{M}}_{\theta \varphi }}d\theta d\varphi \left(\varphi \frac{\partial g_0}{\partial \theta }(\theta ,\varphi )-\theta \frac{\partial g_0}{\partial \varphi }(\theta ,\varphi )\right)\hfill \end{array}$$

(54)

Because the variables $\theta$ and $\varphi$ are only introduced as auxiliary quantities, and have no effect on the evolution in $\mathcal{M}$ and on the evolution of observables, their distribution $g_0$ can be chosen quite freely. Moreover, their roles are fully interchangeable, so we can take distributions such as $g_0(\theta ,\varphi )=g_0(\varphi ,\theta )$, which makes the integral in expression (54) vanish.

The first addend in the expression of ${\tilde{\Omega }}_{{\tilde{f}}_0}$ yields:

$$\begin{array}{ccc}& & {\int }_{{\mathcal{M}}_{\theta \varphi }}d\theta d\varphi g_0(\theta ,\varphi ){\int }_{\mathcal{M}}d\Gamma {\Omega }_{f_0}(\Gamma )\mathcal{O}\left(S^s\Gamma \right)f_0(\Gamma )\hfill \end{array}$$

(55)

$$\begin{array}{ccc}& & ={\int }_{\mathcal{M}}d\Gamma {\Omega }_{f_0}(\Gamma )\mathcal{O}\left(S^s\Gamma \right)f_0(\Gamma )={〈{\Omega }_{f_0}\left(\mathcal{O}\circ S^s\Gamma \right)〉}_0\hfill \end{array}$$

(56)

because $g_0$ is normalized. The remaining term, concerning the time dependent perturbation finally yields:

$$\begin{array}{c}\hfill -{\int }_{{\mathcal{M}}_{\theta \varphi }}d\theta d\varphi g_0(\theta ,\varphi ){\int }_{\mathcal{M}}d\Gamma \frac{1}{f_0(\Gamma )}\sum _{n=-\infty }^{\infty }{w}_n(\theta ,\varphi )\sum _{k=1}^N{\alpha }_{nk}\frac{\partial f_0}{\partial {\Gamma }_k}(\Gamma )\mathcal{O}\left(S^s\Gamma \right)f_0(\Gamma )\end{array}$$

(57)

$$\begin{array}{c}\hfill -\sum _{n=-\infty }^{\infty }\sum _{k=1}^N{\alpha }_{nk}{\int }_{{\mathcal{M}}_{\theta \varphi }}d\theta d\varphi w_n(\theta ,\varphi )g_0(\theta ,\varphi ){\int }_{\mathcal{M}}d\Gamma \frac{\partial f_0}{\partial {\Gamma }_k}(\Gamma )\mathcal{O}\left(S^s\Gamma \right)\end{array}$$

(58)

At this point, it remains to decide how to deal with the different realizations of the perturbation, namely how to treat the stochastic vectors ${\alpha }_n$.

4.2. Single perturbation realization

Here, we treat our system as deterministic, subjected to the time-dependent perturbation that corresponds to a single realization of the stochastic perturbation. Let us label by $\left(j\right)$ this realization. In this case, we can write:

$${\dot{\tilde{\Gamma }}}^{\left(j\right)}={\tilde{G}}^{\left(j\right)}(\Gamma ,\theta ,\varphi )=G(\Gamma )+\sum _{n=-\infty }^{\infty }{\alpha }_n^{\left(j\right)}{w}_n(\theta ,\varphi )$$

(59)

where the expansion coefficients, denoted by ${\alpha }_n^{\left(j\right)}=({\alpha }_{n1}^{\left(j\right)},...{\alpha }_{nN}^{\left(j\right)})$, characterize the chosen perturbation $\left(j\right)$, and the initial condition as well as the magnitude of the perturbation are fixed by $\theta$ and $\varphi$. The phase space expansion rate is given by:

$$\tilde{\Lambda }(\Gamma ,\theta ,\varphi )={\nabla }_{\tilde{\Gamma }}·{\tilde{G}}^{\left(j\right)}={\nabla }_{\Gamma }·G=\Lambda (\Gamma )$$

(60)

where $\Lambda$ denotes the unperturbed phase space volumes variation rate. The new dissipation function is expressed by:

$$\begin{array}{ccc}\hfill & & {\tilde{\Omega }}_{{\tilde{f}}_0}^{\left(j\right)}(\Gamma ,\theta ,\varphi )=\hfill \\ & & {\Omega }_{f_0}(\Gamma )-\frac{1}{f_0(\Gamma )}\sum _{n=-\infty }^{\infty }{w}_n(\theta ,\varphi )\sum _{k=1}^N{\alpha }_{nk}^{\left(j\right)}\frac{\partial f_0}{\partial {\Gamma }_k}(\Gamma )-\frac{1}{g_0(\theta ,\varphi )}\left(\varphi \frac{\partial g_0}{\partial \theta }(\theta ,\varphi )-\theta \frac{\partial g_0}{\partial \varphi }(\theta ,\varphi )\right)\hfill \end{array}$$

(61)

where $g_0$ is the distribution of the initial condition $({\theta }_0,{\varphi }_0)$ of the auxiliary variables. Choosing $g_0$ as in the previous section, the last addend of ${\tilde{\Omega }}_{{\tilde{f}}_0}^{\left(j\right)}$ does not affect the response of observables. Therefore, given an observable $\mathcal{O}$ (which depends on $\Gamma$ only), and denoting by $S^{{\left(j\right)}^s}$ the evolution operator of the perturbed dynamics in $\mathcal{M}$, we can write:

$$\begin{array}{ccc}{\langle \mathcal{O}\rangle }_t\hfill & =& {\langle \mathcal{O}\rangle }_0+{\int }_0^t{〈\left(\mathcal{O}\circ S^{{\left(j\right)}^s}\right){\tilde{\Omega }}_{{\tilde{f}}_0}^{\left(j\right)}〉}_{0}ds\hfill \end{array}$$

(62)

$$\begin{array}{ccc}& =& {\langle \mathcal{O}\rangle }_0+{\int }_0^t{〈\left(\mathcal{O}\circ S^{{\left(j\right)}^s}\right){\Omega }_{f_0}〉}_{0}ds\hfill \\ & & -{\int }_0^{t}ds\sum _{n=-\infty }^{\infty }\sum _{k=1}^N{\alpha }_{nk}^{\left(j\right)}{\int }_{{\mathcal{M}}_{\theta \varphi }}d\theta d\varphi w_n(\theta ,\varphi )g_0(\theta ,\varphi ){\int }_{\mathcal{M}}d\Gamma \frac{\partial f_0}{\partial {\Gamma }_k}(\Gamma )\mathcal{O}\left(S^{{\left(j\right)}^s}\Gamma \right)\hfill \end{array}$$

(63)

$$\begin{array}{ccc}& =& {\langle \mathcal{O}\rangle }_0+{\int }_0^t{〈\left(\mathcal{O}\circ S^{{\left(j\right)}^s}\right){\Omega }_{f_0}〉}_{0}ds\hfill \\ & & -\mathcal{F}\sum _{n=-\infty }^{\infty }\sum _{k=1}^N{\alpha }_{nk}^{\left(j\right)}{\int }_{{\mathcal{M}}_{\theta \varphi }}d\theta d\varphi {\mathcal{F}}^{-{\beta }_n}{\left(\theta -i\varphi \right)}^{{\beta }_n}{g}_0(\theta ,\varphi ){\int }_0^{t}ds{\int }_{\mathcal{M}}d\Gamma \frac{\partial f_0}{\partial {\Gamma }_k}(\Gamma )\mathcal{O}\left(S^{{\left(j\right)}^s}\Gamma \right)\hfill \end{array}$$

(64)

where ${\mathcal{F}}^{-{\beta }_n}{\left(\theta -i\varphi \right)}^{{\beta }_n}$ does not depend on $\mathcal{F}$, cf. Eq.(41). Here, the integral concerning the auxiliary variables gives the mean of $w_n$ with respect to the initial distribution of $(\theta ,\varphi )$. If this is fixed, then this integral simply gives

$$w_n({\theta }_0,{\varphi }_0)=\mathcal{F}{K}_{{\theta }_0{\varphi }_0}^{(n\pm 1/2)\pi /T}$$

(65)

with $K_{{\theta }_0{\varphi }_0}$ a constant that depends on the initial condition, raised to the power $(n-1/2)\pi /T$ if $n>0$ and to the power $(n+1/2)\pi /T$ if $n<0$. But the situation is analogous, if the initial values of $\theta$ and $\varphi$ are distributed with any other density $g_0$. Naturally, letting $\mathcal{F}$ vanish makes the stochastic contribution vanish as well, thus returning to the response formula for time independent perturbations. So far, we have considered the coefficients ${\alpha }_{nk}^{\left(j\right)}$ as fixed, which makes the dynamics look deterministic, and only averages over the initial conditions in the phase space.

4.3. Stochastic coefficients with fixed initial condition

Let us now average over the distribution of the stochastic coefficients, considering first and second cumulant of the fluctuations generated by the stochastic perturbation, about the response due to the deterministic term G. The higher order cumulants vanish, since the stochastic coefficients are Gaussian random variables. Keeping the initial condition $(\Gamma ,\theta ,\varphi )\in \widehat{\mathcal{M}}$ fixed, we denote by ${\langle ·\rangle }^{\left(st\right)}$ the averages made with respect to all realizations of the stochastic perturbation. Consider the equations of motion (42) and its first cumulant:

$${〈\dot{\tilde{\Gamma }}〉}^{\left(st\right)}=\tilde{G}+\sum _{n=-\infty }^{\infty }{〈{\alpha }_n〉}^{\left(st\right)}{w}_n=\tilde{G}$$

(66)

As anticipated, the first cumulant equals the cumulant of the deterministic part of the equations of motion. Moreover, integrating in time, we obtain

$${〈\tilde{\Gamma }\left(t\right)-\tilde{\Gamma }\left(0\right)〉}^{\left(st\right)}={\int }_0^{t}ds{〈\dot{\tilde{\Gamma }}\left(s\right)〉}^{\left(st\right)}={\int }_0^{t}ds\tilde{G}(\Gamma \left(s\right))=\left(\begin{array}{c}\Gamma \left(t\right)\\ \theta \left(t\right)\\ \theta \left(t\right)\end{array}\right)-\left(\begin{array}{c}{\Gamma }_0\\ {\theta }_0\\ {\theta }_0\end{array}\right)$$

(67)

For the second cumulant, let us consider for instance the component ${\Gamma }_l$ of $\Gamma$. We find:

$${〈{\dot{\Gamma }}_l^2〉}^{\left(st\right)}-{〈{\dot{\Gamma }}_l〉}^{{\left(st\right)}^2}={〈{\dot{\Gamma }}_l^2〉}^{\left(st\right)}-G_l^2$$

(68)

Then, let us compute the following:

$$\begin{array}{ccc}& & {〈{\dot{\Gamma }}_l\left(t\right){\dot{\Gamma }}_l\left(t^{\prime }\right)〉}^{\left(st\right)}={〈\left[G_l(\Gamma \left(t\right))+\sum _{n=-\infty }^{\infty }{\alpha }_{nl}{w}_n(\theta \left(t\right),\varphi \left(t\right))\right]\left[G_l(\Gamma (t’))+\sum _{n=-\infty }^{\infty }{\alpha }_{nl}{w}_n(\theta (t’),\varphi (t’))\right]〉}^{\left(st\right)}\hfill \\ & & ={〈G_l(\Gamma \left(t\right))G_l(\Gamma \left(t^{\prime }\right))〉}^{\left(st\right)}+G_l(\Gamma \left(t\right))\sum _{n=-\infty }^{\infty }{〈{\alpha }_{nl}〉}^{\left(st\right)}{w}_n(\theta \left(t^{\prime }\right),\varphi \left(t^{\prime }\right))\hfill \\ & & +G_l(\Gamma \left(t^{\prime }\right))\sum _{n=-\infty }^{\infty }{〈{\alpha }_{nl}〉}^{\left(st\right)}{w}_n(\theta \left(t\right),\varphi \left(t\right))\hfill \\ & & +\sum _{n=-\infty }^{\infty }\sum _{k=-\infty }^{\infty }{〈{\alpha }_{nl}{\alpha }_{kl}〉}^{\left(st\right)}{w}_n\left(\theta \left(t\right),\varphi \left(t\right)\right)w_k\left(\theta (t’),\varphi (t’)\right)\hfill \end{array}$$

(69)

$$\begin{array}{ccc}& & =G_l\left(\Gamma \left(t\right)\right)G_l\left(\Gamma \left(t^{\prime }\right)\right)+\sum _{n=-\infty }^{\infty }{〈{\alpha }_{nl}^2〉}^{\left(st\right)}{w}_n\left(\theta \left(t\right),\varphi \left(t\right)\right)w_n\left(\theta (t’),\varphi (t’)\right)\hfill \end{array}$$

(70)

Indeed, recalling the fact that the random coefficients $Z_k$ are delta correlated, we have

$${〈{\alpha }_{kl}{\alpha }_{n-k,l}〉}^{\left(st\right)}=\left\{\begin{array}{cc}\frac{1}{4i^2}\frac{T}{{\pi }^2}\frac{{〈Z_{kl}{Z}_{n-k,l}〉}^{\left(st\right)}}{c_k{c}_{n-k}}=-\frac{T}{{\left(2\pi \right)}^2}\frac{{\delta }_{k,n-k}}{c_k{c}_{n-k}}\hfill & ifk,n>0\hfill \\ \frac{1}{4i^2}\frac{2T}{{\pi }^2}\frac{{〈Z_{-kl}{Z}_{-n+k,l}〉}^{\left(st\right)}}{c_k{c}_{n-k}}=-\frac{T}{{\left(2\pi \right)}^2}\frac{{\delta }_{-k,-n+k}}{c_k{c}_{n-k}}\hfill & ifk,n<0\hfill \end{array}\right.$$

(71)

so ${\delta }_{k,n-k}=1$ and ${\delta }_{-k,-n+k}=1$ if and only if $n=2k$. Consequently

$${〈{\alpha }_{kl}^2〉}^{\left(st\right)}=-\frac{T}{{\left(2\pi \right)}^2}\frac{1}{c_k^2}.$$

(72)

Now, setting $t=t^{\prime }$, we obtain:

$${〈{\dot{\Gamma }}_l{\left(t\right)}^2〉}^{\left(st\right)}=G_l{(\Gamma \left(t\right))}^2+\sum _{n=-\infty }^{\infty }{〈{\alpha }_{nl}^2〉}^{\left(st\right)}{w}_n{\left(\theta \left(t\right),\varphi \left(t\right)\right)}^2$$

(73)

and, as a result, the second cumulant takes the following simple form:

$$\begin{array}{ccc}\hfill {〈{\dot{\Gamma }}_l^2\left(t\right)〉}^{\left(st\right)}-{〈{\dot{\Gamma }}_l\left(t\right)〉}^{{\left(st\right)}^2}& =& G_l{\left(\Gamma \left(t\right)\right)}^2+\sum _{k=-\infty }^{\infty }{〈{\alpha }_{kl}^2〉}^{\left(st\right)}{w}_k{\left(\theta \left(t\right),\varphi \left(t\right)\right)}^2-G_l{\left(\Gamma \left(t\right)\right)}^2\hfill \end{array}$$

(74)

$$\begin{array}{ccc}& =& \sum _{k=-\infty }^{\infty }{〈{\alpha }_k^2〉}^{\left(st\right)}{w}_k{\left(\theta \left(t\right),\varphi \left(t\right)\right)}^2\hfill \end{array}$$

(75)

Integrating in time expression (70) we get:

$$\begin{array}{ccc}& & {〈{\int }_0^{t}ds{\dot{\Gamma }}_l\left(s\right){\int }_0^{t^{\prime }}d{s}^{\prime }{\dot{\Gamma }}_l(s’)〉}^{\left(st\right)}={〈\left({\Gamma }_l\left(t\right)-{\Gamma }_l\left(0\right)\right)\left({\Gamma }_l(t’)-{\Gamma }_l\left(0\right)\right)〉}^{\left(st\right)}\hfill \\ & & ={\int }_0^{t}ds{\int }_0^{t^{\prime }}d{s}^{\prime }\left[G_l(\Gamma \left(s\right))G_l(\Gamma \left(s^{\prime }\right))+\sum _{k=-\infty }^{\infty }{〈{\alpha }_{kl}^2〉}^{\left(st\right)}{w}_k\left(\theta \left(s\right),\varphi \left(s\right)\right)w_k\left(\theta (s’),\varphi (s’)\right)\right]\hfill \end{array}$$

(76)

and note that Eqs.(40) and (41) imply:

$$\frac{d}{ds}{w}_n\left(\theta \left(s\right),\varphi \left(s\right)\right)=i{\beta }_n{w}_n\left(\theta \left(s\right),\varphi \left(s\right)\right)$$

(77)

so that we can write

$${\int }_0^t{w}_n\left(\theta \left(s\right),\varphi \left(s\right)\right)ds=\frac{1}{i{\beta }_n}\left[w_n\left(\theta \left(t\right),\varphi \left(t\right)\right)-w_n\left({\theta }_0,{\varphi }_0\right)\right]$$

(78)

Then, we obtain:

$$\begin{array}{ccc}& & {〈\left({\Gamma }_l\left(t\right)-{\Gamma }_l\left(0\right)\right)\left({\Gamma }_l\left(t^{\prime }\right)-{\Gamma }_l\left(0\right)\right)〉}^{\left(st\right)}=\hfill \\ & & {\int }_0^{t}ds{\int }_0^{t^{\prime }}d{s}^{\prime }{G}_l\left(\Gamma \left(s\right)\right)G_l\left(\Gamma \left(s^{\prime }\right)\right)+\sum _{k=-\infty }^{\infty }\frac{{〈{\alpha }_{kl}^2〉}^{\left(st\right)}}{{\left(i{\beta }_k\right)}^2}\left[w_k\left(\theta \left(t\right),\varphi \left(t\right)\right)-w_k\left({\theta }_0,{\varphi }_0\right)\right]\hfill \\ & & \times \left[w_k\left(\theta \left(t^{\prime }\right),\varphi \left(t^{\prime }\right)\right)-w_k\left({\theta }_0,{\varphi }_0\right)\right]\hfill \end{array}$$

(79)

Consequently, observing that

$$\begin{array}{ccc}& & {〈\left({\Gamma }_l\left(t\right)-\Gamma \left(0\right)\right)\left({\Gamma }_l\left(t^{\prime }\right)-{\Gamma }_l\left(0\right)\right)〉}^{\left(st\right)}-{〈{\Gamma }_l\left(t\right)-{\Gamma }_l\left(0\right)〉}^{\left(st\right)}{〈{\Gamma }_l\left(t^{\prime }\right)-{\Gamma }_l\left(0\right)〉}^{\left(st\right)}\hfill \end{array}$$

(80)

$$\begin{array}{ccc}& & ={〈{\Gamma }_l\left(t\right){\Gamma }_l\left(t^{\prime }\right)〉}^{\left(st\right)}-{〈{\Gamma }_l\left(t\right)〉}^{\left(st\right)}{〈{\Gamma }_l\left(t^{\prime }\right)〉}^{\left(st\right)}\hfill \end{array}$$

(81)

the auto-correlation function is expressed by:

$$\begin{array}{cc}\hfill & {〈{\Gamma }_l\left(t\right){\Gamma }_l\left(t^{\prime }\right)〉}^{\left(st\right)}-{〈{\Gamma }_l\left(t\right)〉}^{\left(st\right)}{〈{\Gamma }_l\left(t^{\prime }\right)〉}^{\left(st\right)}=\hfill \\ \hfill & \sum _{k=-\infty }^{\infty }\frac{{〈{\alpha }_{kl}^2〉}^{\left(st\right)}}{{\left(i{\beta }_k\right)}^2}\left[w_k\left(\theta \left(t\right),\varphi \left(t\right)\right)-w_k\left({\theta }_0,{\varphi }_0\right)\right]\left[w_k\left(\theta (t’),\varphi (t’)\right)-w_k\left({\theta }_0,{\varphi }_0\right)\right]\hfill \\ \hfill & +{\int }_0^{t}ds{\int }_0^{t^{\prime }}d{s}^{\prime }\left[G_l(\Gamma \left(s\right))G_l(\Gamma \left(s^{\prime }\right))-G_l(\Gamma \left(s\right))G_l(\Gamma \left(s^{\prime }\right))\right]\hfill \\ \hfill & =\sum _{k=-\infty }^{\infty }\frac{{〈{\alpha }_{kl}^2〉}^{\left(st\right)}}{{\left(i{\beta }_k\right)}^2}\left[w_k\left(\theta \left(t\right),\varphi \left(t\right)\right)-w_k\left({\theta }_0,{\varphi }_0\right)\right]\left[w_k\left(\theta (t’),\varphi (t’)\right)-w_k\left({\theta }_0,{\varphi }_0\right)\right]\hfill \end{array}$$

(82)

and $t=t^{\prime }$ yields:

$$\begin{array}{cc}\hfill {〈\Gamma {\left(t\right)}^2〉}^{\left(st\right)}-{〈\Gamma \left(t\right)〉}^{{\left(st\right)}^2}& =\sum _{k=-\infty }^{\infty }\frac{{〈{\alpha }_{kl}^2〉}^{\left(st\right)}}{{\left(i{\beta }_k\right)}^2}{\left[w_k\left(\theta \left(t\right),\varphi \left(t\right)\right)-w_k\left({\theta }_0,{\varphi }_0\right)\right]}^2\hfill \end{array}$$

(83)

For the dissipation function, the first cumulant is trivial since the stochastic coefficients have zero mean:

$${〈{\tilde{\Omega }}_{f_0}〉}^{\left(st\right)}={〈{\Omega }_{f_0}〉}^{\left(st\right)}$$

(84)

For the second cumulant, one has:

$$\begin{array}{ccc}\hfill {〈{\left({\tilde{\Omega }}_{f_0}-{〈{\tilde{\Omega }}_{f_0}〉}^{\left(st\right)}\right)}^2〉}^{\left(st\right)}& =& {〈{\tilde{\Omega }}_{f_0}^2〉}^{\left(st\right)}-{〈{\tilde{\Omega }}_{f_0}〉}^{{\left(st\right)}^2}\hfill \end{array}$$

(85)

$$\begin{array}{ccc}& =& {〈\left({\Omega }_{f_0}w\frac{\partial ln{f}_0}{\partial \Gamma }\right)\left({\Omega }_{f_0}w\frac{\partial ln{f}_0}{\partial \Gamma }\right)〉}^{\left(st\right)}-{〈{\Omega }_{f_0}〉}^{{\left(st\right)}^2}\hfill \end{array}$$

(86)

$$\begin{array}{ccc}& =& {〈{\Omega }_{f_0}^2〉}^{\left(st\right)}-{〈{\Omega }_{f_0}〉}^2+\sum _{k=-\infty }^{\infty }{〈{\alpha }_k^2〉}^{\left(st\right)}{w}_k^2{\left(\frac{\partial ln{f}_0}{\partial \Gamma }\right)}^2\hfill \end{array}$$

(87)

thanks to the statistical properties of stochastic coefficients.

4.4. Averaging over initial conditions and stochastic coefficients

Given a single realization of the stochastic perturbation, the response formula for an observable $\mathcal{O}$ takes the form:

$${〈\mathcal{O}〉}_t={〈\mathcal{O}〉}_0+{\int }_0^t{〈(\mathcal{O}\circ S^{{\left(j\right)}^s}){\Omega }^{f_0}〉}_{0}ds$$

(88)

If we further average over the realizations of the stochastic process, we use the following notation:

$${〈{〈\mathcal{O}〉}_t〉}^{\left(st\right)}$$

(89)

and we can write:

$$\begin{array}{ccc}\hfill {〈{〈\mathcal{O}〉}_{{\tilde{f}}_t}〉}^{\left(st\right)}& =& {〈{〈\mathcal{O}〉}_{{\tilde{f}}_0}〉}^{\left(st\right)}+{\int }_0^t{〈{〈(\mathcal{O}\circ S^{{\left(j\right)}^s}){\tilde{\Omega }}_{f_0}〉}_{{\tilde{f}}_0}〉}^{\left(st\right)}ds\hfill \\ & =& {〈\mathcal{O}〉}_{{\tilde{f}}_0}+{\int }_0^t{〈{〈(\mathcal{O}\circ S^{{\left(j\right)}^s})\left[{\Omega }_{f_0}-w\frac{\partial ln{f}_0(\Gamma )}{\partial \Gamma }\right]〉}_{{\tilde{f}}_0}^{\left(st\right)}〉}^{\left(st\right)}ds\hfill \\ & =& {〈\mathcal{O}〉}_t-{〈\sum _{n=-\infty }^{+\infty }{\alpha }_n{\int }_0^{t}ds{〈(\mathcal{O}\circ S^s)w_n\frac{\partial ln{f}_0(\Gamma )}{\partial \Gamma }〉}_{{\tilde{f}}_0}〉}^{\left(st\right)}\hfill \end{array}$$

(90)

5. Harmonic oscillator

To illustrate the theory developed above, let us consider some simple example, which can be treated analytically. More complex situations can be dealt with performing numerical integration. Let us consider a harmonic oscillator in one dimension, that is perturbed with a time dependent force. We will compare the the linear and the exact response. Let $H_0=\frac{p^2}{2}+\frac{{\omega }^2{q}^2}{2}$ be the unperturbed Hamiltonian, and consider the following perturbation:

$$H_p=-w\left(t\right)q=-ϵZsin\left(\gamma t\right)q$$

(91)

where $ϵ$ is a pure number, and Z has the dimension of force, that may be deterministic or stochastic. Then, the perturbed energy is given by

$$H(q,p,t)=H_0+H_p=\frac{p^2+{\omega }^2{q}^2}{2}-w\left(t\right)q$$

(92)

and the equations of motion are expressed by:

$$G(q,p,t)=\left(\begin{array}{c}\dot{q}\\ \dot{p}\end{array}\right)=\left(\begin{array}{c}p\\ w\left(t\right)-{\omega }^{2}q\end{array}\right)$$

(93)

We make the system autonomous introducing two new variables $\theta$ and $\varphi$, so that the dynamics is given by:

$$\tilde{G}(\Gamma ,\theta ,\varphi )=\left(\begin{array}{c}\dot{q}\\ \dot{p}\\ \dot{\theta }\\ \dot{\varphi }\end{array}\right)=\left(\begin{array}{c}p\\ w(\theta ,\varphi )-{\omega }^{2}q\\ \varphi \\ -{\gamma }^2\theta \end{array}\right)$$

(94)

and the perturbation can be written as:

$$w\left(t\right)=w\left(\theta \right(t),\varphi (t\left)\right)=-ϵZ\varphi \left(t\right)$$

(95)

The motion is then given by

$${\tilde{S}}_{\omega }^t\left(\begin{array}{c}q\\ p\\ \theta \\ \varphi \end{array}\right)=\left(\begin{array}{c}{C}_{1}cos\left(\omega t\right)+\frac{C_2}{\omega }sin\left(\omega t\right)+\frac{ϵZ}{({\omega }^2-{\gamma }^2)}\left[{\theta }_{0}sin\left(\gamma t\right)-\frac{{\varphi }_0}{\gamma }cos\left(\gamma t\right)\right]\\ -\omega C_{1}sin\left(\omega t\right)+C_{2}cos\left(\omega t\right)+\frac{ϵZ}{({\omega }^2-{\gamma }^2)}\left[{\theta }_0\gamma cos\left(\gamma t\right)+{\varphi }_{0}sin\left(\gamma t\right)\right]\\ {\theta }_{0}cos\left(\gamma t\right)+\frac{{\varphi }_0}{\gamma }sin\left(\gamma t\right)\\ -\gamma {\theta }_{0}sin\left(\gamma t\right)+{\varphi }_{0}cos\left(\gamma t\right)\end{array}\right)$$

(96)

where

$$C_1=\left(q_0+ϵZ\frac{{\varphi }_0}{\gamma ({\omega }^2-{\gamma }^2)}\right),C_2=\left(p_0-ϵZ\frac{\gamma {\theta }_0}{({\omega }^2-{\gamma }^2)}\right)$$

(97)

Let us compare the linear and the exact response of our system to the perturbation. First, suppose that $\widetilde{f_0}$ takes the following form:

$$\widetilde{f_0}(\Gamma ,\theta ,\varphi )=\frac{f_0(\Gamma )}{\sigma }=\frac{1}{\sigma }\frac{e^{-\beta H_0}}{Z_0},$$

(98)

where $\sigma$ is the circumference of the variables $\theta$ and $\varphi$, whose radius depends on the initial condition $\mathcal{F}=\sqrt{{\theta }_0^2+{\varphi }_0^2}$. Then,

$$\tilde{\Lambda }(\Gamma ,\theta ,\varphi )={\nabla }_{\tilde{\Gamma }}·\tilde{G}(\Gamma ,\theta ,\varphi )=\frac{\partial }{\partial q}\dot{q}+\frac{\partial }{\partial p}\dot{p}+\frac{\partial }{\partial \theta }\dot{\theta }+\frac{\partial }{\partial \varphi }\dot{\varphi }=0,$$

(99)

$$\begin{array}{cc}\hfill {\tilde{\Omega }}^{f_0}(\Gamma ,\theta ,\varphi )& =-\tilde{\Lambda }(\Gamma ,\theta ,\varphi )-\tilde{G}(\Gamma ,\theta ,\varphi )·{\nabla }_{\tilde{\Gamma }}ln{\tilde{f}}_0(\Gamma ,\theta ,\varphi )\hfill \\ \hfill & =-\left(p,w(\theta ,\varphi )-{\omega }^{2}q,\varphi ,-{\gamma }^2\theta \right)·\left(\begin{array}{cc}& -\beta {\omega }^{2}q\\ & -\beta p\\ & 0\\ & 0\end{array}\right)\hfill \\ \hfill & =\beta {\omega }^{2}pq+\beta w(\theta ,\varphi )p-\beta {\omega }^{2}qp=\beta w(\theta ,\varphi )p\hfill \end{array}$$

(100)

The exact response theory yields:

$${〈\mathcal{O}〉}_{{\tilde{f}}_t}^{\left(st\right)}={〈\mathcal{O}〉}_{{\tilde{f}}_0}^{\left(st\right)}+{\int }_0^t{〈(\mathcal{O}\circ S_w^s)·{\tilde{\Omega }}^{f_0}〉}_{{\tilde{f}}_0}^{\left(st\right)}ds,$$

(101)

where $S_w^s$ denotes the perturbed dynamics and the ensemble average is taken over the augmented phase space. The linear response theory yields:

$${〈\mathcal{O}〉}_{f_t}^{\left(st\right)}={〈\mathcal{O}〉}_{f_0}^{\left(st\right)}+\beta {\int }_0^t\mathcal{F}{〈(\mathcal{O}\circ S_0^s)·\dot{A}〉}_{f_0}^{\left(st\right)}ds$$

(102)

where $S_0^s$ denote the unperturbed dynamics.

Let us consider, for instance, $\mathcal{O}=q$. To proceed, we need some preliminary results.

$$\begin{array}{cc}\hfill {〈q〉}_{f_0}^{\left(st\right)}& =\int \int dqdp\frac{e^{-\beta H_0}}{Z_0}q=\frac{{\int }_{-\infty }^{\infty }{e}^{-\frac{\beta p^2}{2}}dp{\int }_{-\infty }^{\infty }q{e}^{-\frac{\beta {\omega }^2{q}^2}{2}}dq}{\int \int e^{-\beta H_0}dqdp}\hfill \\ \hfill & =\frac{1}{{\int }_{-\infty }^{\infty }{e}^{-\frac{\beta {\omega }^2{q}^2}{2}}dq}\frac{-e^{-\frac{\beta {\omega }^2{q}^2}{2}}}{\beta {\omega }^2}{|}_{-\infty }^{\infty }=\frac{1}{\sqrt{2\pi /\beta {\omega }^2}}\frac{-e^{-\frac{\beta {\omega }^2{q}^2}{2}}}{\beta {\omega }^2}{|}_{-\infty }^{\infty }=0,\hfill \end{array}$$

(103)

$$\begin{array}{cc}\hfill {〈p〉}_{f_0}^{\left(st\right)}& =\int \int dqdp\frac{e^{-\beta H_0}}{Z_0}p=\frac{{\int }_{-\infty }^{\infty }p{e}^{-\frac{\beta p^2}{2}}dp{\int }_{-\infty }^{\infty }{e}^{-\frac{\beta {\omega }^2{q}^2}{2}}dq}{\int \int e^{-\beta H_0}dqdp}\hfill \\ \hfill & =\frac{1}{{\int }_{-\infty }^{\infty }{e}^{-\frac{\beta p^2}{2}}dp}\frac{-e^{-\frac{\beta p^2}{2}}}{\beta }{|}_{-\infty }^{\infty }=\frac{1}{\sqrt{2\pi /\beta }}\frac{-e^{-\frac{\beta p^2}{2}}}{\beta }{|}_{-\infty }^{\infty }=0,\hfill \end{array}$$

(104)

$$\begin{array}{cc}\hfill {〈p^2〉}_{f_0}^{\left(st\right)}& =\int \int dqdp\frac{e^{-\beta H_0}}{Z_0}{p}^2=\frac{{\int }_{-\infty }^{\infty }{p}^2{e}^{-\frac{\beta p^2}{2}}dp}{{\int }_{-\infty }^{\infty }{e}^{-\frac{\beta p^2}{2}}dp}\hfill \\ \hfill & =\frac{-\frac{p}{\beta }{e}^{-\frac{\beta p^2}{2}}{|}_{-\infty }^{\infty }+{\int }_{-\infty }^{\infty }\frac{1}{\beta }{e}^{-\frac{\beta p^2}{2}}}{\sqrt{2\pi /\beta }}=\frac{\frac{1}{\beta }\sqrt{2\pi /\beta }}{\sqrt{2\pi /\beta }}=k_{b}T,\hfill \end{array}$$

(105)

$$\begin{array}{cc}\hfill {〈q^2〉}_{f_0}^{\left(st\right)}& =\int \int dqdp\frac{e^{-\beta H_0}}{Z_0}{q}^2=\frac{{\int }_{-\infty }^{\infty }{q}^2{e}^{-\frac{\beta {\omega }^2{q}^2}{2}}dq}{{\int }_{-\infty }^{\infty }{e}^{-\frac{\beta {\omega }^2{q}^2}{2}}dq}\hfill \\ \hfill & =\frac{-\frac{q}{\beta {\omega }^2}{e}^{-\frac{\beta {\omega }^2{q}^2}{2}}{|}_{-\infty }^{\infty }+{\int }_{-\infty }^{\infty }\frac{1}{\beta {\omega }^2}{e}^{-\frac{\beta {\omega }^2{q}^2}{2}}}{\sqrt{2\pi /\beta {\omega }^2}}=\frac{\frac{1}{\beta {\omega }^2}\sqrt{2\pi /\beta {\omega }^2}}{\sqrt{2\pi /\beta {\omega }^2}}=\frac{k_{b}T}{{\omega }^2}.\hfill \end{array}$$

(106)

Then, the linear response for a single realization of the perturbation, i.e. for a given value Z, yields:

$$\begin{array}{cc}\hfill {〈q〉}_t^{\left(st\right)}& ={〈q〉}_0^{\left(st\right)}+\beta {\int }_0^t\mathcal{F}{〈(q\circ S_0^s)·\dot{A}〉}_0^{\left(st\right)}ds\hfill \\ \hfill & =\frac{\mathcal{F}Z}{2\omega }\left(\frac{sin\left[\right(\omega -\gamma \left)t\right]}{(\omega -\gamma )}-\frac{sin\left[\right(\omega +\gamma \left)t\right]}{(\omega +\gamma )}\right)\hfill \end{array}$$

(107)

On the other hand, the exact response is expressed by:

$$\begin{array}{cc}\hfill {〈q〉}_{{\tilde{f}}_t}^{\left(st\right)}& ={〈q〉}_{{\tilde{f}}_0}^{\left(st\right)}+{\int }_0^t{〈(q\circ S_w^s)·{\tilde{\Omega }}^{f_0}〉}_{{\tilde{f}}_0}^{\left(st\right)}ds={\int }_0^t{〈(q\circ S_w^s)·\frac{p}{k_{b}T}w(\theta ,\varphi )〉}_{{\tilde{f}}_0}^{\left(st\right)}ds\hfill \\ \hfill & ={\int }_0^{t}ds\int d\theta d\varphi dqdpw(\theta ,\varphi )\frac{1}{\sigma }\frac{e^{-H_0/\left(k_{b}T\right)}}{Z_0}\frac{p}{k_{b}T}[\left(q+\mathcal{F}Z\frac{\varphi }{\gamma ({\omega }^2-{\gamma }^2)}\right)cos\left(\omega s\right)\hfill \\ \hfill & +\frac{1}{\omega }\left(p-\mathcal{F}Z\frac{\gamma \theta }{({\omega }^2-{\gamma }^2)}\right)sin\left(\omega s\right)+\frac{\mathcal{F}Z}{({\omega }^2-{\gamma }^2)}\left(\theta sin\left(\gamma s\right)-\frac{\varphi }{\gamma }cos\left(\gamma s\right)\right)]\hfill \\ \hfill & ={\int }_0^{t}ds\int d\theta d\varphi dqdp\frac{w(\theta ,\varphi )}{k_{b}T}\frac{1}{\sigma }\frac{e^{-H_0/\left(k_{b}T\right)}}{Z_0}\frac{p^2}{\omega }sin\left(\omega s\right)\hfill \\ \hfill & =-\frac{\mathcal{F}Z}{\sigma \gamma }{\int }_0^{t}ds\frac{sin\left(\omega s\right)}{\omega }\int \int d\theta d\varphi (-\gamma \theta sin\left(\gamma s\right)+\varphi cos\left(\gamma s\right))=0\hfill \end{array}$$

(108)

because the integral in $\theta$ and $\varphi$ vanishes (see Appendix for details). Therefore, linear and exact response differ, under the chosen conditions. If, on the other hand, the initial condition of the auxiliary variables is fixed to the single point ${\theta }_0=1$ and ${\varphi }_0=0$, we have $\mathcal{F}=1$, and

$$\widetilde{f_0}(\Gamma ,\theta ,\varphi )=\delta ({\theta }_0-1)\delta \left({\varphi }_0\right)f_0(\Gamma )=\delta ({\theta }_0-1)\delta \left({\varphi }_0\right)\frac{e^{-\beta H_0}}{Z_0}$$

(109)

The exact response is then expressed by:

$$\begin{array}{cc}\hfill {〈q〉}_{{\tilde{f}}_t}^{\left(st\right)}& ={〈q〉}_{{\tilde{f}}_0}^{\left(st\right)}+{\int }_0^t{〈(q\circ S_w^s)·\widetilde{{\Omega }^{f_0}}〉}_{{\tilde{f}}_0}^{\left(st\right)}ds={\int }_0^t{〈(q\circ S_w^s)·\frac{p}{k_{b}T}w(\theta ,\varphi )〉}_{{\tilde{f}}_0}^{\left(st\right)}ds\hfill \\ \hfill & ={\int }_0^{t}ds\int d\theta d\varphi dqdpw(\theta ,\varphi )\delta (\theta -1)\delta \left(\varphi \right)\frac{e^{-H_0/\left(k_{b}T\right)}}{Z_0}\frac{p}{k_{b}T}[\left(q+Z\frac{\varphi }{\gamma ({\omega }^2-{\gamma }^2)}\right)cos\left(\omega s\right)\hfill \\ \hfill & +\frac{1}{\omega }\left(p-Z\frac{\gamma \theta }{({\omega }^2-{\gamma }^2)}\right)sin\left(\omega s\right)+\frac{Z}{({\omega }^2-{\gamma }^2)}\left(\theta sin\left(\gamma s\right)-\frac{\varphi }{\gamma }cos\left(\gamma s\right)\right)]\hfill \\ \hfill & ={\int }_0^{t}ds\int d\theta d\varphi dqdp\frac{w(\theta ,\varphi )}{k_{b}T}\delta (\theta -1)\delta \left(\varphi \right)\frac{e^{-H_0/\left(k_{b}T\right)}}{Z_0}\frac{p^2}{\omega }sin\left(\omega s\right)\hfill \\ \hfill & =-\frac{Z}{\gamma }{\int }_0^{t}ds\frac{sin\left(\omega s\right)}{\omega }\int \int d\theta d\varphi \delta (\theta -1)\delta \left(\varphi \right)(-\gamma \theta sin\left(\gamma s\right)+\varphi cos\left(\gamma s\right))\hfill \\ \hfill & =Z{\int }_0^{t}ds\frac{sin\left(\omega s\right)}{\omega }sin\left(\gamma s\right)=\frac{Z}{2\omega }\left(\frac{sin\left[\right(\omega -\gamma \left)t\right]}{(\omega -\gamma )}-\frac{sin\left[\right(\omega +\gamma \left)t\right]}{(\omega +\gamma )}\right)\hfill \end{array}$$

(110)

where the integral on the curve in $d\theta$ and $d\varphi$ gives the length of the curve that is $\sigma$. This is the same as the linear response with $({\theta }_0,{\varphi }_0)$ uniformly distributed in the unit circle. If we take as an observable $\mathcal{O}=p$ we get something similar. For the linear response we have:

$$\begin{array}{cc}\hfill {〈p〉}_t^{\left(st\right)}& ={〈p〉}_0^{\left(st\right)}+\beta {\int }_0^t\mathcal{F}\left(s\right){〈(p\circ S_0^s)·\dot{A}〉}_0^{\left(st\right)}ds\hfill \\ \hfill & =\mathcal{F}Z{\int }_0^{t}dssin\left(\gamma s\right)\int \int dqdp\frac{e^{-H_0/\left(k_{b}T\right)}}{Z_0}\frac{p}{k_{b}T}\left[-\omega qsin\left(\omega s\right)+pcos\left(\omega s\right)\right]\hfill \\ \hfill & =\mathcal{F}Z{\int }_0^{t}dssin\left(\gamma s\right)\int \int dqdp\frac{e^{-H_0/\left(k_{b}T\right)}}{Z_0}\frac{p^2}{k_{b}T}cos\left(\omega s\right)\hfill \\ \hfill & =\mathcal{F}Z{\int }_0^{t}dssin\left(\gamma s\right)cos\left(\omega s\right)=\frac{ϵZ}{2}{\int }_0^t\left[sin\left[\right(\gamma -\omega \left)s\right]+sin\left[\right(\gamma +\omega \left)s\right]\right]ds\hfill \\ \hfill & =\frac{\mathcal{F}Z}{2}\left(-\frac{cos\left[\right(\gamma -\omega \left)t\right]}{(\gamma -\omega )}-\frac{cos\left[\right(\gamma +\omega \left)t\right]}{(\gamma +\omega )}+2\right)\hfill \end{array}$$

(111)

and for the exact response we have:

$$\begin{array}{cc}\hfill {〈p〉}_{{\tilde{f}}_t}^{\left(st\right)}& ={〈p〉}_{{\tilde{f}}_0}^{\left(st\right)}+{\int }_0^t{〈(p\circ S_w^s)·{\tilde{\Omega }}^{f_0}〉}_{{\tilde{f}}_0}^{\left(st\right)}ds={\int }_0^t{〈(p\circ S_w^s)·\frac{p}{k_{b}T}w(\theta ,\varphi )〉}_{{\tilde{f}}_0}^{\left(st\right)}ds\hfill \\ \hfill & ={\int }_0^{t}ds\int d\theta d\varphi dqdpw(\theta ,\varphi )\frac{1}{\sigma }\frac{e^{-H_0/\left(k_{b}T\right)}}{Z_0}\frac{p}{k_{b}T}[-\omega \left(q+\mathcal{F}Z\frac{\varphi }{\gamma ({\omega }^2-{\gamma }^2)}\right)sin\left(\omega s\right)\hfill \\ \hfill & +\left(p-\mathcal{F}Z\frac{\gamma \theta }{({\omega }^2-{\gamma }^2)}\right)cos\left(\omega s\right)+\frac{\mathcal{F}Z}{({\omega }^2-{\gamma }^2)}\left(\theta \gamma cos\left(\gamma s\right)+\varphi sin\left(\gamma s\right)\right)]\hfill \\ \hfill & ={\int }_0^{t}ds\int d\theta d\varphi dqdp\frac{w(\theta ,\varphi )}{k_{b}T}\frac{1}{\sigma }\frac{e^{-H_0/\left(k_{b}T\right)}}{Z_0}\frac{p^2}{\omega }cos\left(\omega s\right)\hfill \\ \hfill & =-\frac{\mathcal{F}Z}{\sigma \gamma }{\int }_0^{t}ds\frac{cos\left(\omega s\right)}{\omega }\int \int d\theta d\varphi (-\gamma \theta sin\left(\gamma s\right)+\varphi cos\left(\gamma s\right))=0\hfill \end{array}$$

(112)

Again, the same response is obtained only for fixed initial condition $({\theta }_0=1,{\varphi }_0=0)$. This is not strange, it happens when the exact response amounts to a first order perturbation, which is the case in special situations [20]. However, in general, that condition is not verified. For instance, consider the observable $\mathcal{O}=qp$, the linear theory yields:

$$\begin{array}{cc}\hfill {〈qp〉}_{f_t}^{\left(st\right)}& ={〈qp〉}_{f_0}^{\left(st\right)}+\beta {\int }_0^t\mathcal{F}\left(s\right){〈(q\circ S_0^s)(p\circ S_0^s)·\dot{A}〉}_{f_0}^{\left(st\right)}ds\hfill \\ \hfill & =\mathcal{F}Z{\int }_0^{t}dssin\left(\gamma s\right)\int \int dqdp\frac{e^{-H_0/\left(k_{b}T\right)}}{Z_0}\frac{p}{k_{b}T}[\left(qcos\left(\omega s\right)+\frac{p}{\omega }sin\left(\omega s\right)\right)\hfill \\ \hfill & \times \left(-\omega qsin\left(\omega s\right)+pcos\left(\omega s\right)\right)]=0\hfill \end{array}$$

(113)

because integrating the terms $p{q}^2$, $q{p}^2$, $p^3$ over $dqdp$ gives zero. For the exact response, we have:

$$\begin{array}{cc}\hfill {\langle qp\rangle }_{{\tilde{f}}_t}& ={\langle qp\rangle }_{{\tilde{f}}_0}+{\int }_0^{t}ds{\langle w(\theta ,\varphi )\frac{p}{k_{b}T}(S_w^S\circ q)(S_w^s\circ p)\rangle }_{{\tilde{f}}_0}\hfill \\ \hfill & ={\int }_0^{t}ds〈w(\theta ,\varphi )\frac{p}{k_{b}T}\left\{\right[C_{1}cos\left(\omega t\right)+\frac{C_2}{\omega }sin\left(\omega t\right)+\frac{ϵZ}{({\omega }^2-{\gamma }^2)}({\theta }_{0}sin\left(\gamma t\right)+\hfill \\ \hfill & {-\frac{{\varphi }_0}{\gamma }cos\left(\gamma t\right)\left)\right]\left[-\omega C_{1}sin\left(\omega t\right)+C_{2}cos\left(\omega t\right)+\frac{ϵZ}{({\omega }^2-{\gamma }^2)}\left({\theta }_0\gamma cos\left(\gamma t\right)+{\varphi }_{0}sin\left(\gamma t\right)\right)\right]\}〉}_{{\tilde{f}}_0}\hfill \end{array}$$

(114)

Taking into account all the terms that cancel under integration over $dqdp$, and those that cancel by integrating over $d\theta d\varphi$, what remains are elliptic integrals over $d\theta d\varphi$, that can be solved numerically, and integrals over t of sine and cosine terms. The result does not vanish, as can be seen by inspection. In this example, even fixing the initial conditions as $({\theta }_0=1,{\varphi }_0=0)$ and using the Dirac delta distribution, we get again a non-vanishing result:

$$\begin{array}{ccc}\hfill & & {〈qp〉}_{{\tilde{f}}_t}^{\left(st\right)}=\hfill \\ & & -\frac{\mathcal{F}Z}{\gamma }{\int }_0^{t}ds\frac{-2\mathcal{F}Z\gamma }{\omega ({\omega }^2-{\gamma }^2)}sin\left(\omega s\right)cos\left(\omega s\right)\int \int d\theta d\varphi \delta (\theta -1)\delta \left(\varphi \right)\left(-\gamma \theta sin\left(\gamma s\right)+\varphi cos\left(\gamma s\right)\right)\hfill \\ & & -\frac{\mathcal{F}Z}{\gamma }{\int }_0^{t}ds\frac{\mathcal{F}Z\gamma }{\omega ({\omega }^2-{\gamma }^2)}cos\left(\gamma s\right)\int \int d\theta d\varphi \delta (\theta -1)\delta \left(\varphi \right)\left(-\gamma \theta sin\left(\gamma s\right)+\varphi cos\left(\gamma s\right)\right)\hfill \\ & & =-\frac{2{\mathcal{F}}^2{Z}^2}{\omega ({\omega }^2-{\gamma }^2)}\left[{\displaystyle \frac{\left(\gamma -2\omega \right)sin\left(\gamma t+2\omega t\right)}{4{\gamma }^2-16{\omega }^2}}-{\displaystyle \frac{\left(\gamma +2\omega \right)sin\left(\gamma t-2\omega t\right)}{4{\gamma }^2-16{\omega }^2}}+\frac{{sin}^2\gamma t}{4\gamma }\right]\hfill \end{array}$$

(115)

In this case, the linear response yields zero because the exact response is second order in the perturbation magnitude $\mathcal{F}$. When that is sufficiently small, the terms of order $O\left({\mathcal{F}}^2\right)$ can be neglected, but not in general. Furthermore, forced oscillators may undergo resonance phenomena, greatly amplifying the amplitude of their oscillations even if they are driven by very weak forces. This is captured by the denominators of expression (116), which are missing in the linear response. Near the resonances, the linear response is bound to fail, even for small driving, and the exact formula is necessary.

6. Concluding remarks

In this paper we have provided a generalization to stochastic processes of the exact response theory developed for deterministic molecular dynamics models, also known as TTCF in the field of molecular dynamics, where it mostly concerns time independent perturbations. This has been done starting from the previously developed time dependent response theory of Ref. [11], expressing the stochastic perturbation as a Karhunen-Loève expansion, and then adding two auxiliary variables that make autonomous the corresponding time dependent system of ordinary differential equations. The main ingredient of the exact response theory, the dissipation function, has then been derived in the enlarged phase space. This allows us to treat processes whose initial conditions are variously distributed, adding flexibility to the common response theory, that considers fixed initial perturbations. Moreover, it preserves the notion of time reversibility in the equations of motion, which is a key ingredient of numerous statistical mechanics results, notably Onsager Reciprocal Relations and Fluctuation Relations.

We illustrated the theory applying it to a one-dimensional harmonic oscillator, perturbed by a sinusoidal force, whose amplitude can be either deterministic or stochastic. This is enough to compare the exact response with the linear response. Although in special situations the two responses may coincide, up to a certain degree defining the range of applicability of the linear theory, they differ in general, even in such simple situations. The exact formula is necessary for large perturbations, as expected, but also for small perturbations in presence of resonance phenomena. Furthermore, we have shown that different results are obtained for different distributions of the initial perturbation, something that is not obvious within the standard linear response theory. Indeed, enlarging the phase space to make the system autonomous, allows us to choose the distribution of the initial condition of the perturbation, referring to different kinds of experiments, in which the initial perturbation can be either deterministically fixed, as commonly assumed, or distributed in various ways. Note, we explicitly treated only zero mean perturbations, because the possible net perturbation may be included in the autonomous part of the dynamics, for which the exact theory has been developed in the past.

The above comes in addition to the fact that many theoretical results can be derived from the exact response formula, before they can be numerically computed in simulations, e.g. [7]; that linear response cannot treat phase transitions [17]; and that various studies show an impressively better performance of the TTCF formula, compared to other averaging methods for time independent perturbations [6,19]. We expect this to be the case with our approach, when dealing with time dependent (deterministic or stochastic) perturbations. That constitutes the purposes of future works.