On some aspects of the response to stochastic and deterministic forcings

Response theory aims at predicting the how the statistical properties of a system under investigation changes as a result of the application of (weak) external stimuli. In linear approximation, the cornerstone result to this end is the fluctuation-dissipation theorem (FDT) [Kub57, Kub66], which relates the statistics of the unforced fluctuations and the relaxation properties of the system with its response to external forcings. Specifically, the FDT allows for writing the Green function describing the response of the system as a correlation between observables of the unperturbed state. Hence, the FDT classically provides a powerful tool for predicting the properties of forced fluctuations from those of the unforced ones. See [LV07, MPRV08, SV19] for a comprehensive summary of the main results and more modern perspectives, and [LC12] for an extension of the theory to the nonlinear case.

The applicability of the FDT relies on the existence of a smooth (with respect to Lebesgue) invariant measure, which allows for applying an integration by parts in the derivation of the final formulas. In the case of deterministic dynamics, this is viable only for systems with vanishing average phase space contraction, as in the case of physical system at equilibrium. When non-conservative forces come into play, phase space volume contracts, making the invariant measure describing the statistical steady state singular with respect to the Lebesgue measure. In this case, one needs to resort to more general response formulas. Ruelle [Rue98a, Rue09] introduced a response theory able to predict the change in the statistical properties of axiom A systems due to small perturbation to its dynamics. Directly implementing the Ruelle response formulas is particularly nontrivial in the case of deterministic chaotic systems, because of the radically different properties of the tangent space in its unstable and stable directions [AM07]. Fortunately, recent contributions based on adjoint and shadowing methods seem to provide a convincing way forward [Wan13, CW20, Ni20, SW22]. Ruelle's response theory has been then reformulated and extended through the viewpoint of functional analysis by studying via a perturbative expansion the corrections to the transfer operator [Bal00] due to the applied forcing [GL06, BL07, Bal08].

Systems of stochastic differential equations (SDEs) [Kha66, KP92, Arn98] typically possess smooth invariant probability densities, which ensures the applicability of the FDT and the straightforward derivation of linear response formulas [HT82, HM10, Pav14]. See also the recent work by [DG22], where the Feynman–Kac formalism is used to derive expressions for the covariances of general additive dynamical functionals for stochastic dynamical systems, and the comprehensive review by [Mae20], where a trajectory-based approach on the computation of the response focusing on computing and interpreting the time-symmetric and time-antisymmetric contributions is discussed in detail.

Linear response theory has been used to study sensitivity and transport properties of a plethora of physical problems and systems, including stochastic resonance [GHJM98], optical materials [LSPV05], simple toy models of chaotic dynamics [Rei02, CS07, Luc09, LS11], Markov chains [Mit04, Luc16, ADF18, SGLCG21], neural networks [Ces19], turbulence [Kan20], galactic dynamics [BT08], financial markets [PTSSG+21], plasma [Nam76], interacting identical agents [LPZ20, ZLP21, ACCR22], optomechanical systems [MDN21], and the climate system [Lei75, NBH93, LRL17, LLR20, BDH21], just to name a few; see also the many examples discussed in [Ött05], the theoretical developments presented in [Luc08, BM13, BG18, Luc18, SV19, WG19], and the recent special issue edited by Gottwald [Got20].

When considering SDEs, the theory of Markov semigroups and of their spectral decomposition (see [TLLD18, CDN+20]) can be applied to explain the properties of correlation functions. In the formal limit of weak noise, the spectral properties resemble those of the underlying deterministic component (cf [KL99]) but, as mentioned above, in this case formal problems emerge when trying to relate correlations with response formulas. Nevertheless, the importance of the spectral features of Markov semigroups in gauging the smoothness of the response is crucial, even for high-dimensional and physically relevant systems. Indeed, the work of [CNK+14] is the first study to relate the eigenvalues of reduced Markov representations of geophysical flows with the rough/smooth dependence on model parameter uncertainty.

The first goal of this note is to study the linear response of SDEs to perturbations that can affect either the drift term or the stochastic forcing. Our strategy will rely on investigating the impact of applying perturbations to the dynamics in terms of corrections to the Fokker–Planck equation or to the backward-Kolmogorov equation, which describe the evolution of the smooth measures and observables of the system, respectively. This will leads us to computing the Green functions of the system and subsequently express it in terms of Koopman modes and the corresponding eigenvalues [Mez05, BMM12], taking advantage of the spectral theory of Markov semigroups. The basic derivations are revisited here in section 2. The advantage of taking this angle on the problem is that one can relate the relaxation and oscillatory time scales of the Koopman modes, which are intrinsic to the unperturbed system, with the time-dependent response of general observables for general forcings. This also clarifies the close relationship that exists between the presence of slow decay of correlations and, loosely speaking, the strong sensitivity of a system to perturbations, pointing to its proximity to a critical transition [TLD18, TLLD18, GL20, SBBBC09, LHKHL08, Kue11]. Additionally, we derive similarly general explicit expressions for the spectral properties of the frequency-dependent response of the system. Such a connection is particularly illuminating in the case of systems whose Koopman modes possess a clear distinction between the point and the residual spectrum [SGLCG21, CDN+20, theorems 4–6], as discussed in detail below.

We will also explore another aspect of the relationship between the response of a system to perturbations and stochasticity. Stochastic perturbations to deterministic systems have been widely studied in the literature [Kif88, FW98] and the investigation of the weak-noise limit in the case the dynamics of the underlying system is chaotic is of key importance for constructing the so-called physical measure [ER85, Rue89]. The reader is also referred to [Gar09, chapter 7] for an introduction to the impact of stochastic perturbations and to [BK98, KL99] for abstract results concerning the transfer operator. In section 3 we shall, nonetheless, follow previous work based on response theory whereby, taking advantage of Ruelle's formalism, the effects of adding noise to a deterministic system are computed [Luc12]. The same question is tackled in [Abr17], although the author resorts to expanding the perturbed flow in terms of the tangent map leading to an apparently different formula to that found in [Luc12]. The relationship between these two formulas will be discussed below in proposition 3.1, elucidating the importance of considering the correct interpretation for the noise term.

There is a somewhat surprising lack of studies focusing on the systematic investigation of the response of higher-order correlations to forcings. As shown in [LW17], this problem is instead of relevance also for understanding the sensitivity to coarse-grained representations of the system of interest. As additional output of our work, in section 4 we clarify that it is possible to extend the results presented in [LW17] to the case of perturbations acting on either the drift of the noise law of a system, whether its background dynamics is deterministic or stochastic. We then achieve a possibly useful extension of the response formulas able to encompass many cases of interest where one wishes to study the sensitivity of the correlations of a system to perturbations.

Finally, in section 5 we draw our conclusions and present perspectives for future work. Appendix A revisits the homogeneous equation for the response derived originally in [Ken71] where, here, we rewrite such equation in terms of the stochastic analysis language employed throughout this note.

In this section, the construction of the Green function for an SDE is revisited using Fokker–Planck operator expansions. We start by considering the d-dimensional Itô SDE with deterministic drift $\mathbf{F}:{\mathbb{R}}^{d}\to {\mathbb{R}}^{d}$, perturbation field $\mathbf{G}:{\mathbb{R}}^{d}\to {\mathbb{R}}^{d}$, modulated by the bounded function $g:\mathbb{R}\to \mathbb{R}$, and noise law given by the d × p matrix ${\Sigma}\in {\mathbb{R}}^{d\times p}$, p ⩾ 1:

$$\begin{equation}\mathrm{d}\mathbf{x}(t)=\left[\mathbf{F}(\mathbf{x})+\varepsilon g(t)\mathbf{G}(\mathbf{x})\right]\mathrm{d}t+{\Sigma}(\mathbf{x})\mathrm{d}{W}_{t},\end{equation} \tag{ 1 }$$

where W_t is a p-dimensional Wiener process and ɛ is a real number. We assume that the SDE (1) generates a process x(t) in ${\mathbb{R}}^{d}$ for an initial condition at a certain value of time. The SDE is autonomous when ɛ = 0. When ɛ ≠ 0, a time-modulated perturbation vector field G is activated through the function g. The evolution of probability density functions associated with the SDE (1) is given by the Fokker–Planck equation:

$$\begin{equation}{\partial }_{t}\rho =\mathcal{L}(t)\rho =-\nabla \cdot \left(\mathbf{F}\rho \right)-\varepsilon g(t)\nabla \cdot \left(\mathbf{G}\rho \right)+\frac{1}{2}{\nabla }^{2}:\left({\Sigma}{{\Sigma}}^{\top }\rho \right),\end{equation} \tag{ 2 }$$

where ':' denotes the Frobenius inner product, ΣΣ^⊤ is the noise covariance matrix and the operator $\mathcal{L}$ can be expressed in terms of the two operators: ${\mathcal{L}}_{0}:D({\mathcal{L}}_{0})\subseteq \mathcal{B}\to \mathcal{B}$ and ${\mathcal{L}}_{1}:D({\mathcal{L}}_{1})\subseteq \mathcal{B}\to \mathcal{B}$ so that $\mathcal{L}(t)={\mathcal{L}}_{0}+\varepsilon g(t){\mathcal{L}}_{1}$ and their domains, denoted by D(⋅), are densely defined in a suitable Banach space $\mathcal{B}$ and are assumed to satisfy $D(\mathcal{L}(t))=D({\mathcal{L}}_{0})=D({\mathcal{L}}_{1})$. These operators are defined as:

$$\begin{equation}{\mathcal{L}}_{0}\rho =-\nabla \cdot \left(\mathbf{F}\rho \right)+\frac{1}{2}{\nabla }^{2}:\left({\Sigma}{{\Sigma}}^{\top }\rho \right);\end{equation} \tag{ 3a }$$

$$\begin{equation}{\mathcal{L}}_{1}\rho =-\nabla \cdot \left(\mathbf{G}\rho \right),\end{equation} \tag{ 3b }$$

for every ρ in $D({\mathcal{L}}_{0})$. We immediately identify that ${\mathcal{L}}_{1}$ is the operator corresponding to the perturbation introduced in equation (1) and hence, we refer to it as the perturbation operator. Because the right-hand side (rhs) of equation (1) is linear, its solutions when ɛ = 0 are given by the exponential map family ${\left\{{\mathrm{e}}^{t{\mathcal{L}}_{0}}\right\}}_{t\geqslant 0}$, which we call the semigroup generated by ${\mathcal{L}}_{0}$. In fact, ${\mathcal{L}}_{0}$ generates a strongly continuous semigroup so that the properties of the exponential function are inherited [EN00]. It is noted, however, that $\mathcal{L}$ is time-dependent, so the way it generates an evolution semigroup is by means of the time-ordered exponential operator [Gil17].

The regularity of the transition probabilities in equation (2) and encoded in $\mathcal{B}$ is not granted by merely adding noise, but one has to ensure that the latter spreads in all directions. To this end, Hörmander's condition [Hör67] is invoked so that at any point in phase space $\mathcal{X}$, the tangent linear is recovered by the directions generated by the second order differential operator in ${\mathcal{L}}_{0}$ stemming from the noise. Hence, the smoothness and boundedness of the coefficients of equation (1) together with Hörmander's condition make sure that we are dealing with an hypoelliptic diffusion process and ensure the existence of smooth transition probabilities solving equation (2) and an invariant probability density function ρ₀ that solves ${\mathcal{L}}_{0}{\rho }_{0}\equiv 0$ [Pav14].

Since the SDE is non-autonomous, the density functions ρ are dependent on the phase variable x at time t and the initial condition x₀ at time s. If one, instead, assumes that the initial condition x₀ at time s = 0 is distributed according to the unperturbed state ρ₀ (ɛ = 0), we can understand the densities in equation (2) as

$$\begin{equation}\rho (\mathbf{x},t){:=}\int \rho (\mathbf{x},\mathbf{y},t,0){\rho }_{0}(\mathbf{y})\mathrm{d}\mathbf{y},\end{equation} \tag{ 4 }$$

so that ρ(x, t) determines the probability of encountering the process in x at time t for a ρ₀-distributed initial condition; see also [Pav14, chapter 2] for more general considerations.

The goal now is to calculate the distribution function for the perturbed system i.e., when ɛ ≠ 0. To this end, we shall solve equation (2) for successive orders of ɛ. Let us assume that a solution of equation (2) ρ can be written as:

$$\begin{equation}{\rho }_{\varepsilon }={\rho }_{0}+\varepsilon {\rho }_{1}+\mathcal{O}({\varepsilon }^{2}).\end{equation} \tag{ 5 }$$

Such an asymptotic expansion is the starting point of virtually all linear response formulas for statistical mechanical systems. We assume here that our systems of interest comply with the conditions discussed in [HM10]. A discussion of the radius of convergence of the expansion in the context of finite-state Markov processes is presented in [Luc16, SGL20].

While ρ₀ is provided, ρ₁ is obtained by plugging equation (5) into (2) and gathering, first, the ɛ-order terms and applying the variation-of-parameters formula:

$$\begin{equation}{\rho }_{1}(\cdot ,t)={\int }_{0}^{t}{\mathrm{e}}^{(t-s){\mathcal{L}}_{0}}g(s){\mathcal{L}}_{1}{\rho }_{0}\mathrm{d}s,\end{equation} \tag{ 6 }$$

where the time-modulation g appears in the integrand. The kth order element of equation (5) is obtained analogously by gathering ɛ^k -terms [Kub57, Rue98b, Luc08]. The expectation value of a general observable Ψ can be written as:

$$\begin{equation}\left\langle {\Psi},\rho (\cdot ,t)\right\rangle =\sum\limits _{k=0}^{\infty }\,{\delta }^{(k)}\left[{\Psi}\right](t)=\int {\Psi}(\mathbf{x}){\rho }_{0}(\mathbf{x})\mathrm{d}\mathbf{x}+\varepsilon \int {\Psi}(\mathbf{x}){\rho }_{1}(\mathbf{x},t)\mathrm{d}\mathbf{x}+\mathcal{O}\left({\varepsilon }^{2}\right).\end{equation} \tag{ 7 }$$

Hence, after integrating by parts the linear response accounting for the first order corrections is:

$$\begin{equation}{\delta }^{(1)}[{\Psi}](t)=\varepsilon {\int }_{-\mathrm{\infty }}^{\mathrm{\infty }}g(s)\int {\Theta}(t-s){\mathrm{e}}^{(t-s){\mathcal{L}}_{0}^{\ast }}{\Psi}(\mathbf{x}){\mathcal{L}}_{1}{\rho }_{0}(\mathbf{x})\mathrm{d}\mathbf{x}\mathrm{d}s=\varepsilon (g\bullet \mathcal{G})(t),\end{equation} \tag{ 8 }$$

where Θ(t) is the Heaviside distribution, '*' denotes the dual representation, and $\bullet$ indicates the convolution product. Specifically, the generator of the unperturbed backward-Kolmogorov equation ${\mathcal{L}}_{0}^{\ast }$ and the analogous perturbation operator ${\mathcal{L}}_{1}^{\ast }$ act on measurable observables Ψ and can be written as

$$\begin{equation}{\mathcal{L}}_{0}^{\ast }{\Psi}=\mathbf{F}\cdot \nabla {\Psi}+\frac{1}{2}{\Sigma}{{\Sigma}}^{\top }:{\nabla }^{2}{\Psi};\end{equation} \tag{ 9a }$$

$$\begin{equation}{\mathcal{L}}_{1}^{\ast }{\Psi}=\mathbf{G}\cdot \nabla {\Psi}.\end{equation} \tag{ 9b }$$

Furthermore, the Green function $\mathcal{G}(t)={\Theta}(t)\int {\mathrm{e}}^{t{\mathcal{L}}_{0}^{\ast }}{\Psi}(\mathbf{x}){\mathcal{L}}_{1}{\rho }_{0}(\mathbf{x})\mathrm{d}\mathbf{x}$ has a non-negative support because of causality [Rue09, Luc18]. Equation (8) implies that once the Green function is known or computed off-line from the unperturbed regime, the linear response can be readily calculated regardless of the time-modulation g.

Remark 2.1. It is also possible to study the impact of applying a perturbation to the stochastic term in equation (1) as follows:

$$\begin{equation}\mathrm{d}\mathbf{x}(t)=\mathbf{F}(\mathbf{x})\mathrm{d}t+{\Sigma}(\mathbf{x})\mathrm{d}{W}_{t}+\varepsilon g(t){\Gamma}(\mathbf{x})\mathrm{d}{W}_{t},\end{equation} \tag{ 10 }$$

where ${\Gamma}\in {\mathbb{R}}^{d\times p}$ is a d × p matrix defining the correction to the background noise law. Note that we allow for the presence of nontrivial time modulations in the properties of the noise through the bounded function g(t). Up to first order in ɛ, the Fokker–Planck equation reads as follows:

$$\begin{equation}{\partial }_{t}\rho =\mathcal{L}(t)\rho =-\nabla \cdot \left(\mathbf{F}\rho \right)-\varepsilon g(t)\frac{1}{2}{\nabla }^{2}:\left(\left({\Sigma}{{\Gamma}}^{\top }+{\Gamma}{{\Sigma}}^{\top }\right)\rho \right)+\frac{1}{2}{\nabla }^{2}:\left({\Sigma}{{\Sigma}}^{\top }\rho \right).\end{equation} \tag{ 11 }$$

In this case, one has $\mathcal{L}(t)={\mathcal{L}}_{0}+\varepsilon g(t){\mathcal{L}}_{1}$ where ${\mathcal{L}}_{1}$ is:

$$\begin{equation}{\mathcal{L}}_{1}\rho =\frac{1}{2}{\nabla }^{2}:\left(\left({\Sigma}{{\Gamma}}^{\top }+{\Gamma}{{\Sigma}}^{\top }\right)\rho \right).\end{equation} \tag{ 12 }$$

Hence, by plugging this expression of ${\mathcal{L}}_{1}$ into equations (6)–(8) one can compute the first order correction to the statistical properties of the system due to perturbations to the stochastic diffusion. The same formulas—apart from the specific expression of ${\mathcal{L}}_{1}$—can be used to evaluate at leading the impact of either deterministic or stochastic forcing onto the system. Additionally, because of linearity, the impact of deterministic (as in equation (1)) and stochastic perturbations (as in equation (10)) can be treated separately and then added up to account for the response of the system to a combined perturbation.

Note that if ΣΓ^⊤ + Σ^⊤Γ vanishes (and in particular if Σ ≡ 0, i.e. if the background dynamics is deterministic) the leading order response to a stochastic perturbation is proportional to ɛ², see [Luc12, Abr17] and discussion in section 3. Note also that the dual operator ${\mathcal{L}}_{1}^{\ast }$ acting on measurable observables Ψ can be written as:

$$\begin{equation}{\mathcal{L}}_{1}^{\ast }{\Psi}=\frac{1}{2}\left({\Sigma}{{\Gamma}}^{\top }+{\Gamma}{{\Sigma}}^{\top }\right):{\nabla }^{2}{\Psi}.\end{equation} \tag{ 13 }$$

Additionally, note also that if we include in equation (10) extra noise terms of the form ɛf(t)Δ(x)dV_t , where ${\Delta}\in {\mathbb{R}}^{d\times q}$ is a d × q matrix, f is a bounded function, and V_t is a q-dimensional Wiener process with no correlation with the previously defined W_t p-dimensional Wiener process, their leading-order impact on the measure is $\mathcal{O}({\varepsilon }^{2})$, because the first order term of equations (12) and (13) would be zero.

2.1. Spectral decomposition of the linear response

Under mild conditions on the SDE (1) and its associated Fokker–Planck equation (2), the semigroup ${\left\{{\mathrm{e}}^{t{\mathcal{L}}_{0}}\right\}}_{t\geqslant 0}$ enjoys an exponential convergence to equilibrium due to ${\mathcal{L}}_{0}$ having a nontrivial set of dominating eigenvalues [CDN+20, theorems 4–6]. The target of this section is, thus, to decompose the Green function $\mathcal{G}$ in terms of the spectrum of the unforced generator ${\mathcal{L}}_{0}$ of the Fokker–Planck equation in order to extract the relaxation timescales of the response. For this, we shall further require the semigroup ${\left\{{\mathrm{e}}^{t{\mathcal{L}}_{0}}\right\}}_{t\geqslant 0}$ to be quasi-compact, this is, it approaches the space of compact operators as t tends to infinity in the operator norm [EN00]. Quasi-compact semigroups are constituted by operators having a strictly negative essential growth bound r_ess, ⁴ implying that spectral components of the semigroup with modulus larger than ${\mathrm{e}}^{{r}_{\text{ess}}t}$ are eigenvalues of finite algebraic multiplicity, which are related to those of the generator by means of the spectral mapping theorem [EN00, chapter IV]. Namely, if λ_j is an eigenvalue of ${\mathcal{L}}_{0}$ with eigenfunction ψ_j , so is ${\mathrm{e}}^{{\lambda }_{j}t}$ of ${\mathrm{e}}^{t{\mathcal{L}}_{0}}$ relative to the same eigenfunction. Therefore, if ${\left\{{\lambda }_{j}\right\}}_{j=1}^{M}$ are M eigenvalues of finite algebraic multiplicity a_j such that λ_j > r_ess the evolution operator ${\mathrm{e}}^{t{\mathcal{L}}_{0}}$ is decomposed as:

$$\begin{equation}{\mathrm{e}}^{t{\mathcal{L}}_{0}}=\sum\limits _{j=1}^{M}\,{T}_{j}(t)+\mathcal{R}(t),\end{equation} \tag{ 14 }$$

where T_j (t) denotes the contribution relative to the eigenvalue λ_j and $\mathcal{R}$ is the operator accounting for the essential spectrum. The operator T_j (t) is defined as:

$$\begin{equation}{T}_{j}(t)=\sum\limits _{k=0}^{{a}_{j}-1}\,{\mathrm{e}}^{{\lambda }_{j}t}\frac{{t}^{k}}{k!}{\left({\mathcal{L}}_{0}-{\lambda }_{j}\right)}^{k}{{\Pi}}_{j},\end{equation} \tag{ 15 }$$

where Π_j denotes the spectral projector corresponding to the eigenvalue λ_j . It is worth stressing that equation (15) applies also to non-diagonalisable operators. Indeed, in case the eigenvalue λ_j has non-unit multiplicity, i.e. a_j > 1, equation (15) accounts for the generalised eigenfunctions and possible 'Jordan blocks'. For notational convenience, we shall further assume that ${\lambda }_{0}=0 > \mathfrak{R}\mathfrak{e}{\lambda }_{1}\geqslant \mathfrak{R}\mathfrak{e}{\lambda }_{2}\geqslant \dots$. The hope is that, as t tends to infinity, the contributions from the essential spectrum decay. In fact, if $r > \mathrm{sup}\enspace \left\{{r}_{\text{ess}}\right\}\cup \left\{\mathfrak{R}\mathfrak{e}\lambda :\lambda \in \sigma ({\mathcal{L}}_{0}){\backslash}\left\{{\lambda }_{0},\dots ,{\lambda }_{M}\right\}\right\}$ there exists a constant C > 0 such that:

$$\begin{equation}{\Vert}\mathcal{R}(t){\Vert}\leqslant C{\mathrm{e}}^{rt},\end{equation} \tag{ 16 }$$

for all positive values of t. This ensures that, as time goes to infinity, the norm of the residual operator decays to zero. The finiteness of the eigenvalues assumed here implies that all of them, but for λ₀ = 0, are contained in a strip of the complex plane that is at a distance $\gamma =-\mathfrak{R}\mathfrak{e}{\lambda }_{1}$ from the imaginary axis. The number γ is called the spectral gap of ${\mathcal{L}}_{0}$ and it guarantees that correlations decay exponentially fast [CDN+20, Bal00, section 1.3]. It must be noted, though, that the point spectrum need not be, in general, finite, so that it is possible that it accumulates around the leading zero eigenvalue. This pathological case is responsible for subexponential decay rates [Rue86] and a signal of the system approaching a critical point, as is the case of the pitchfork bifurcation [GNPT95] or the critical transitions between competing attractors in atmospheric circulation models [TLLD18], or the rough/smooth parameter dependence in geophysical flows [CNK+14].

The spectral decomposition of quasi-compact operators allows us to naturally link the theory of the Green function with the eigenvalues of the generator of the unperturbed Fokker–Planck semigroup ${\mathcal{L}}_{0}$. We remark that in the case of smooth invariant measures, such eigenvalues agree with those of the generator of the backward-Kolmogorov equation ${\mathcal{L}}_{0}^{\ast }$. If we assume for the sake of simplicity that our observable of interest Ψ projects entirely onto the point spectrum of ${\mathcal{L}}_{0}^{\ast }$, we obtain the following closed formula:

$$\begin{equation}{\delta }^{(1)}[{\Psi}](t)=\varepsilon (g\bullet \mathcal{G})(t)=\varepsilon {\int }_{-\mathrm{\infty }}^{\mathrm{\infty }}g(t-s)\mathcal{G}(s)\mathrm{d}s\end{equation} \tag{ 17a }$$

$$\begin{equation}=\varepsilon {\int }_{-\infty }^{\infty }g(t-s)\int {\Theta}(s){\mathrm{e}}^{s{\mathcal{L}}_{0}^{\ast }}{\Psi}(\mathbf{x}){\mathcal{L}}_{1}{\rho }_{0}(\mathbf{x})\mathrm{d}\mathbf{x}\mathrm{d}s\end{equation} \tag{ 17b }$$

$$\begin{equation}=\varepsilon {\int }_{-\infty }^{\infty }g(t-s)\int \left({\Theta}(s)\sum\limits _{j=1}^{M}{\mathrm{e}}^{{\lambda }_{j}s}{{\Pi}}_{j}^{\ast }{\Psi}(\mathbf{x})\right){\mathcal{L}}_{1}{\rho }_{0}(\mathbf{x})\mathrm{d}s\mathrm{d}\mathbf{x}\end{equation} \tag{ 17c }$$

$$\begin{equation}=\varepsilon \sum\limits _{j=1}^{M}\,\sum\limits _{k=1}^{{a}_{j}-1}\,{\alpha }_{j}^{(k)}\frac{1}{k!}{\int }_{-\infty }^{\infty }g(t-s){\Theta}(s){\mathrm{e}}^{{\lambda }_{j}s}{s}^{k}\mathrm{d}s\end{equation} \tag{ 17d }$$

$$\begin{equation}=\varepsilon \sum\limits _{j=1}^{M}\;\sum\limits _{k=1}^{{a}_{j}-1}\;{\alpha }_{j}^{(k)}\frac{1}{k!}[g\bullet ({\Theta}(\circ ){\mathrm{e}}^{{\lambda }_{j}\circ }{\circ }^{k})](t),\end{equation} \tag{ 17e }$$

where the coefficients ${\alpha }_{j}^{(k)}$ are defined as:

$$\begin{equation}{\alpha }_{j}^{(k)}=\int {({\mathcal{L}}_{0}^{\ast }-{\lambda }_{j})}^{k}{{\Pi}}_{j}^{\ast }{\Psi}(\mathbf{x}){\mathcal{L}}_{1}{\rho }_{0}(\mathbf{x})\mathrm{d}\mathbf{x},\end{equation} \tag{ 18 }$$

and ${{\Pi}}_{j}^{\ast }$ is the projector of the observable Ψ on the jth mode of the generator ${\mathcal{L}}_{0}^{\ast }$. The α coefficients weight the contribution of the various modes of the response for a given combination of observable and forcing.

The calculations presented in equation (17) assume that the observable of interest Ψ projects entirely onto the point-spectrum of ${\mathcal{L}}_{0}$, but in general, Ψ can possess a significant component contained in the residual spectrum. The quasi-compact assumption, however, ensures that after 1/r time units, the residual contribution to the response δ⁽¹⁾[Ψ](t) will be negligible, where r was introduced in equation (16). In essence, the further away the residual spectrum is from the imaginary axis, the less significant will its contribution be to the response as a function of time. Hence, the reason for neglecting the residual spectrum is twofold. First, the calculations become exact and, second, we emphasize the role of the dominating eigenvalues in determining the decay rates of the response function.

As clear from the discussion above, these results apply to either a deterministic or a stochastic perturbation to the dynamics of an SDE or to a combination thereof; one just needs to use the right definition for the operator ${\mathcal{L}}_{1}$. Notice that the contribution of the eigenvalue 0 is not present because it is simple and its associated dual eigenfunction is constant, making α₀ = 0. We understand here that the spectral gap gauges the sensitivity of the system with respect to external perturbations: if 0 < γ ≪ 1 the linear response δ⁽¹⁾[Ψ](t) to an impulse defined by g(t) = δ(t), where δ(t) is Dirac's delta, will decay very slowly to zero for any observable Ψ. In physical terms, this means that the negative feedbacks of the system are very weak.

At a practical level, when a system is weakly forced so that its response is in the linear regime, δ⁽¹⁾[Ψ](t) can be observed from numerical simulations. In particular, if g(t) = δ(t), the perturbation would be equivalent to an ɛ-sized shift in the initial condition of the unperturbed system. This comment raises the question of using empirical measurements of the response ${\delta }^{(1)}\left[{\Psi}\right](t)$ to find the decomposition of the Green function in terms of the point spectrum, which can be thence convoluted against different time-modulations. In effect, equation (17) provides a basis for the empirically observed response ${\delta }^{(1)}\left[{\Psi}\right](t)$ to project onto.

In Fourier domain, the convolution product defining ${\delta }^{(1)}\left[{\Psi}\right](t)$ is transformed into a regular product:

$$\begin{equation}\mathfrak{F}[{\delta }^{(1)}[{\Psi}]](\omega )=\varepsilon \mathfrak{F}[g\bullet \mathcal{G}](\omega )=\varepsilon \mathfrak{F}[g](\omega )\mathfrak{F}[\mathcal{G}](\omega ),\end{equation} \tag{ 19 }$$

where $\mathfrak{F}$ indicates the application of the Fourier transform, ω is a real frequency and $\mathfrak{F}\left[\mathcal{G}\right]$ is the susceptibility function which is analytic in the upper complex plane, by virtue of the causality constraint. Furthermore, if one assumes polynomial decay of $\mathfrak{F}\left[\mathcal{G}\right](\omega )$, Kramers–Kronig relations [LSPV05], can be derived, which link the imaginary and real parts of the susceptibility function via Cauchy integrals [Rue98b, Luc08]. Under suitable integrability assumptions, one can, moreover, link the logarithm of the modulus of the susceptibility with its real part, and then use Kramers–Kronig relations to find the imaginary counterpart [Luc12].

Exploiting the spectral decomposition of quasi-compact semigroups and assuming that a given observable does not project onto the residual spectrum, the Fourier transform of the linear response is written as:

$$\begin{equation}\mathfrak{F}\left[\mathcal{G}\right](\omega )=\sum\limits _{j=1}^{M}\,\sum\limits _{k=1}^{{a}_{j}-1}\,{\alpha }_{j}^{(k)}\frac{1}{k!}\mathfrak{F}\left[{\Theta}{\mathrm{e}}^{{\lambda }_{j}{\circ}}{{\circ}}^{k}\right](\omega )=\sum\limits _{j=1}^{M}\,\sum\limits _{k=1}^{{a}_{j}-1}\,\frac{{\alpha }_{j}^{(k)}}{{\left(i\omega -{\lambda }_{j}\right)}^{k+1}}.\end{equation} \tag{ 20 }$$

The susceptibility function $\mathfrak{F}\left[\mathcal{G}\right]$ can be meromorphically extended to all values of ω in $\mathbb{C}$ such that $\mathfrak{I}\mathfrak{m}\omega > {\omega }_{\text{ess}}$, where the poles [Rue86] correspond to those of the resolvent norm ${\Vert}R(\cdot ,{\mathcal{L}}_{0}){\Vert}$ which are precisely given by the eigenvalues of ${\mathcal{L}}_{0}$ with the exception of that at ω = 0 by the observation made below equation (18). It is then clear that the location of the poles in the susceptibility function depends, uniquely, on the underlying system and not on the time-modulation function g, on the applied forcing, or on the observable in question.

Correspondingly, equation (20) implies the existence of resonances in the response for $\omega =\mathfrak{I}\mathfrak{m}{\lambda }_{j}$, j = 1, ..., M. Neglecting for the moment the possible existence of nonunitary algebraic multiplicities, the visibility of each resonance (the height of the corresponding peak when considering real values for ω) will, instead, depend on the width, determined by $\mathfrak{R}\mathfrak{e}{\lambda }_{j}$ (which depends neither on the observable nor on the applied forcing), and by the weight given by the spectral coefficient ${\alpha }_{j}^{(1)}$ (which depends, instead, on both the observable and the applied forcing). In the limit of a vanishing spectral gap γ, the spectrum will be dominated by the resonance associated to λ₁ unless ${\alpha }_{j}^{(1)}$ vanishes.

Note also that, in general, the assumption that the contribution of the residual spectrum to the response operators is negligible for large times implies that its contribution to the susceptibility is negligible in the low-frequency range.

Remark 2.2. The existence of a smooth invariant density for equation (1) is enough for the fluctuation–dissipation theorem to hold [Pav14]. This implies that the spectral decomposition of correlation functions in stochastic systems, see [CDN+20], would equally yield equations (17) and (20). However, in the limit of small noise, density representations and amenable functional settings for non-conservative can fail to exist. Attention is drawn, instead, to the Green function which is possible to define more generally using Ruelle's response theory for deterministic flows [Rue09, Luc08]. The ultimate goal is therefore to replace ρ₀(x)dx by ρ₀(dx) in equation (8) and hope that the spectral decomposition of equation (17) holds. If such is the case, the relaxation rates would still be encoded in the Green function through the generator ${\mathcal{L}}_{0}^{\ast }$ while correlation functions would cease to explain the response at leading order.

The theory of random dynamical systems regards stochastic flows as transformations of phase space parameterised in time but also in the noise realisation [Arn98]. In this sense, one can view noise as the inhomogeneous component of an equation that makes it, loosely speaking, non-autonomous. Interestingly, when the noise component is relatively weak, it can be regarded as a perturbation to the otherwise unperturbed deterministic system. Consequently, the leading order modification to the system's statistics is expected to be captured through linear response theory, in the spirit of the previous section.

Indeed, such is the case, although two main issues arise. First, the reference dynamics is deterministic, so that the matrix Σ in equation (10) vanishes and hence, the unperturbed system will not, in general, possess an absolutely continuous invariant measure with respect to Lebesgue. This prevents the possibility of using density functions when studying the evolution of the measure. Second, the Green function can be seamlessly applied to derive realisation-dependent response formulas although these do not agree with a perturbation expansion of the Fokker–Planck equation. This section, therefore, aims at (i), determining the leading order correction to the statistics using operator relations and (ii), proving that Green function and operator based response formulas are equal only when the former is taken in the Stratonovich setting.

We proceed by considering a stochastically perturbed system that obeys the following Itô SDE:

$$\begin{equation}\mathrm{d}\mathbf{x}(t)=\mathbf{F}(\mathbf{x})\mathrm{d}t+\varepsilon {\Gamma}(\mathbf{x})\mathrm{d}{W}_{t},\end{equation} \tag{ 21 }$$

where $\mathbf{F}:{\mathbb{R}}^{d}\to {\mathbb{R}}^{d}$ is the drift, ${\Gamma}:{\mathbb{R}}^{d}\to {\mathbb{R}}^{d\times p}$ is a perturbation matrix, and W_t denotes a p-dimensional Wiener process with correlation matrix equal to the identity. Because we are dealing with deterministic background dynamics, we shall avoid density representations and use, instead, the backward-Kolmogorov equation with describes the evolution of expectation values of observables Ψ:

$$\begin{equation}{\partial }_{t}{\Psi}(\cdot ,t)=\left({\mathcal{L}}_{0}^{\ast }+{\varepsilon }^{2}{\mathcal{L}}_{2}^{\ast }\right){\Psi}(\cdot ,t).\end{equation} \tag{ 22 }$$

Here, we are using the adjointness between the Fokker–Planck and the backward-Kolmogorov equation (also invoked in equation (8)), whereby the operators ${\mathcal{L}}_{0}^{\ast }$ and ${\mathcal{L}}_{2}^{\ast }$ are given by:

$$\begin{equation}{\mathcal{L}}_{0}^{\ast }{\Psi}=\mathbf{F}\cdot \nabla {\Psi};\end{equation} \tag{ 23a }$$

$$\begin{equation}{\mathcal{L}}_{2}^{\ast }{\Psi}=\frac{1}{2}\left({\Gamma}{{\Gamma}}^{\top }\right):{\nabla }^{2}{\Psi},\end{equation} \tag{ 23b }$$

for Ψ in $D({\mathcal{L}}_{2}^{\ast })$. Here we immediately observe that stochastic forcing enters at second order in ɛ when considering the operator representation of the stochastic process (21). Using the semigroup expansions of section 2, we can derive the leading order response of the steady state to the introduction of noise:

$$\begin{align}\hfill {\delta }^{(2)}\left[{\Psi}\right](t)& ={\varepsilon }^{2}\int {\rho }_{0}(\mathrm{d}\mathbf{x}){\int }_{0}^{t}{\mathcal{L}}_{2}^{\ast }{\Psi}(\mathbf{x}(t))\mathrm{d}s=\frac{{\varepsilon }^{2}}{2}\int {\rho }_{0}(\mathrm{d}\mathbf{x})\times {\int }_{0}^{t}\left({\Gamma}(\mathbf{x}){{\Gamma}}^{\top }(\mathbf{x})\right):{\nabla }^{2}{\Psi}(\mathbf{x}(s))\left.\right)\mathrm{d}s.\hfill \end{align} \tag{ 24 }$$

Note that we have evaluated the deterministic Koopman operator ${\mathrm{e}}^{t{\mathcal{L}}_{0}^{\ast }}$ to the observable Ψ to give ${\mathrm{e}}^{t{\mathcal{L}}_{0}^{\ast }}{\Psi}(\mathbf{x})={\Psi}(\mathbf{x}(t))$, where x(t) solves equation (21) for ɛ = 0.

The formula in equation (24) coincides with what was obtained in [Abr17], where instead, the author resorted to computing the impact of the perturbation directly on the tangent space in order to evaluate $\left({\Gamma}(\mathbf{x}){{\Gamma}}^{\top }(\mathbf{x})\right):{\nabla }^{2}{\Psi}(\mathbf{x}(t))$. An approach based on second-order response theory was followed by [Luc12], yet yielding a different formula compared to Abramov's [Abr17]. In such a work, the Wiener increments dW_t were treated as the time-modulation of a sequence of p perturbations with spatial patterns defined by the columns of the Γ matrix. To illustrate this idea, it is useful to rewrite equation (21) as a series of applied vector fields modulated by Wiener increments:

$$\begin{equation}\mathrm{d}\mathbf{x}=\mathbf{F}(\mathbf{x})\mathrm{d}t+\varepsilon \sum\limits _{i=1}^{p}\,{g}_{i}(t){{\Gamma}}_{:,i}(\mathbf{x}),\end{equation} \tag{ 25 }$$

where Γ_:,i is the ith column of Γ and g_i (t) are assumed to satisfy $\mathbb{E}\left[{g}_{i}(t)\right]=0$ and $\mathbb{E}\left[{g}_{k}(t){g}_{l}(s)\right]=\delta (t-s){\delta }_{k,l}$ for every i, k, l = 1, ..., p, where δ_k,l is the Kronecker delta. For a given realisation of the noise indexed by σ, one can apply equation (8) to obtain the linear response correction:

$$\begin{equation}{\delta }_{\sigma }^{(1)}\left[{\Psi}\right](t){:=}\varepsilon \sum\limits _{i=1}^{p}{\int }_{-\infty }^{\infty }{\mathcal{G}}_{i}(s){g}_{i}(t-s)\mathrm{d}s,\end{equation} \tag{ 26 }$$

where ${\mathcal{G}}_{i}(t)={\Theta}(t)\int {\rho }_{0}(\mathrm{d}\mathbf{x}){\mathcal{L}}_{1}^{\ast }{\mathrm{e}}^{t{\mathcal{L}}_{0}^{\ast }}{\Psi}(\mathbf{x})$ is the Green function associated with the vector field Γ_:,i for every i. Since, g_i corresponds to a Wiener increment, taking averages over all realisations σ makes the linear responses vanish: $\mathbb{E}\left[{\delta }_{\sigma }^{(1)}\left[{\Psi}\right]\right]=0$. Contrarily, the second order response ${\delta }_{\sigma }^{(2)}\left[{\Psi}\right]$ survives the averaging and is simplified to a single-time-variable integral of the form:

$$\begin{equation}{\tilde{\delta }}^{(2)}\left[{\Psi}\right](t){:=}\mathbb{E}\left[{\delta }_{\sigma }^{(2)}\left[{\Psi}\right]\right]=\frac{{\varepsilon }^{2}}{2}\sum\limits _{k=1}^{p}\int {\rho }_{0}(\mathrm{d}\mathbf{x}){\int }_{-\infty }^{\infty }{\Theta}(s){\mathcal{L}}_{1,k}^{\ast }{\mathcal{L}}_{1,k}^{\ast }{\mathrm{e}}^{s{\mathcal{L}}_{0}^{\ast }}{\Psi}(\mathbf{x})\mathrm{d}s,\end{equation} \tag{ 27 }$$

where ${\mathcal{L}}_{1,i}^{\ast }={{\Gamma}}_{:,i}\cdot \nabla$. To compare equation (27) with (24), it is handy to disentangle the column vector fields Γ_:,i from the correlation matrix Γ so that, after some vector-matrix manipulations, one deduces, instead, that:

$$\begin{equation}{\delta }^{(2)}\left[{\Psi}\right](t)=\frac{{\varepsilon }^{2}}{2}\sum\limits _{k=1}^{p}\int {\rho }_{0}(\mathrm{d}\mathbf{x}){\int }_{0}^{t}\left({{\Gamma}}_{:,k}(\mathbf{x}){{\Gamma}}_{:,k}^{\top }(\mathbf{x})\right):{\nabla }^{2}{\Psi}(\mathbf{x}(s))\mathrm{d}s.\end{equation} \tag{ 28 }$$

Remark 3.1. In case the applied fields Γ_;,i in equation (25) are not modulated by independent time-functions g_i (t), we would have that $\mathbb{E}\left[{g}_{k}(t){g}_{l}(s)\right]=\delta (t-s){c}_{k,l}$, for some correlation factor c_k,l . This would make cross-terms appear, so that

$$\begin{align}\hfill {\tilde{\delta }}^{(2)}\left[{\Psi}\right](t)& {:=}\mathbb{E}\left[{\delta }_{\sigma }^{(2)}\left[{\Psi}\right]\right]=\frac{{\varepsilon }^{2}}{2}\sum\limits _{k=1}^{p}\,\sum\limits _{l=1}^{p}\,{C}_{k,l}\int {\rho }_{0}(\mathrm{d}\mathbf{x})\times {\int }_{-\infty }^{\infty }{\Theta}(s){\mathcal{L}}_{1,k}^{\ast }{\mathcal{L}}_{1,l}^{\ast }{\mathrm{e}}^{s{\mathcal{L}}_{0}^{\ast }}{\Psi}(\mathbf{x})\mathrm{d}s,\hfill \end{align} \tag{ 29 }$$

and,

$$\begin{equation}{\delta }^{(2)}\left[{\Psi}\right](t)=\frac{{\varepsilon }^{2}}{2}\sum\limits _{k=1}^{p}\,\sum\limits _{l=1}^{p}\,{C}_{k,l}\int {\rho }_{0}(\mathrm{d}\mathbf{x}){\int }_{0}^{t}\left({{\Gamma}}_{:,k}(\mathbf{x}){{\Gamma}}_{:,l}^{\top }(\mathbf{x})\right):{\nabla }^{2}{\Psi}(\mathbf{x}(x))\mathrm{d}s.\end{equation} \tag{ 30 }$$

However, it is easy to recast this case into the form of equation (21) upon a modification of the covariance matrix.

We remark that the approach described in the present paper using operator expansions yields a different formula to that of [Luc12] for the seemingly same question. In order to understand this discrepancy, we need to clarify our interpretation of the noise term. The algebraic manipulations used in [Luc12] rely on using the standard rules of calculus, which is appropriate if one considers Stratonovich's interpretation of noise. And, indeed, the concept behind the use of the second order response approach relied on treating the noise terms as fast forcings, thus considering the physical limit of time-scale separation going to infinity. This, again, points in the direction of Stratonovich's interpretation [PS08, Pav14].

To test our hypothesis, we shall investigate whether the approach taken in [Luc12] can be recast more rigorously according to the formalism used in this paper. To this end and for clarity, let us consider equation (25) where a single perturbation Γ_:,1 is applied, namely, p = 1. This means we can rewrite equation (25) as:

$$\begin{equation}\mathrm{d}\mathbf{x}(t)=\mathbf{F}(\mathbf{x})\mathrm{d}t+\varepsilon \tilde{{\Gamma}}H\mathrm{d}{W}_{t},\end{equation} \tag{ 31 }$$

where x(t) is in ${\mathbb{R}}^{d}$, $\tilde{{\Gamma}}=\mathrm{diag}({{\Gamma}}_{:,1})$, W_t is a d-dimensional Wiener process and H is a d × d matrix defined as:

$$\begin{equation}H={\mathbb{1}}_{d}{\mathbf{e}}_{1}^{\top }=\left[\begin{matrix}\hfill 1\hfill \\ \hfill {\vdots}\hfill \\ \hfill 1\hfill \end{matrix}\right][1,0,\dots ,0].\end{equation} \tag{ 32 }$$

We proceed by expanding the operator product ${\mathcal{L}}_{1,1}^{\ast }{\mathcal{L}}_{1,1}^{\ast }$, which is instrumental in equation (27). Letting f denote a generic twice differentiable function and $\nabla \times$ the curl, we have:

$$\begin{equation}{\mathcal{L}}_{1,1}^{\ast }{\mathcal{L}}_{1,1}^{\ast }f={{\Gamma}}_{:,1}\cdot \mathrm{\nabla }({{\Gamma}}_{:,1}\cdot \mathrm{\nabla }f)\end{equation} \tag{ 33a }$$

$$\begin{align}\hfill & ={{\Gamma}}_{:,1}\cdot \left[\left({{\Gamma}}_{:,1}\cdot \nabla \right)\nabla f+\left(\nabla f\cdot \nabla \right){{\Gamma}}_{:,1}\right]\hfill \\ \hfill & \quad +{{\Gamma}}_{:,1}\cdot \left[{{\Gamma}}_{:,1}\times (\nabla \times \nabla f)+\nabla f\times (\nabla \times {{\Gamma}}_{:,1})\right]\hfill \end{align} \tag{ 33b }$$

$$\begin{equation}={{\Gamma}}_{:,1}\cdot \left[\left({{\Gamma}}_{:,1}\cdot \nabla \right)\nabla f+\left(\nabla f\cdot \nabla \right){{\Gamma}}_{:,1}+\nabla f\times (\nabla \times {{\Gamma}}_{:,1})\right]\end{equation} \tag{ 33c }$$

$$\begin{equation}={{\Gamma}}_{:,1}\cdot \left({{\Gamma}}_{:,1}\cdot \nabla \right)\nabla f+{{\Gamma}}_{:,1}\cdot \left(\nabla f\cdot \nabla \right){{\Gamma}}_{:,1}+{{\Gamma}}_{:,1}\cdot \left(\nabla f\times \nabla \times {{\Gamma}}_{:,1}\right)\end{equation} \tag{ 33d }$$

$$\begin{equation}={{\Gamma}}_{:,1}\cdot \left({{\Gamma}}_{:,1}\cdot \nabla \right)\nabla f+{{\Gamma}}_{:,1}\cdot \left(\nabla f\cdot \nabla \right){{\Gamma}}_{:,1}-\nabla f\cdot \left({{\Gamma}}_{:,1}\times \nabla \times {{\Gamma}}_{:,1}\right)\end{equation} \tag{ 33e }$$

$$\begin{equation}=\left({{\Gamma}}_{:,1}{{\Gamma}}_{:,1}^{\top }\right):{\nabla }^{2}f+{{\Gamma}}_{:,1}\cdot \left(\nabla f\cdot \nabla \right){{\Gamma}}_{:,1}-\nabla f\cdot \left({{\Gamma}}_{:,1}\times \nabla \times {{\Gamma}}_{:,1}\right).\end{equation} \tag{ 33f }$$

We immediately observe that if $\nabla {{\Gamma}}_{:,1}\equiv 0$, the response formula derived in equations (24) and (28) coincide. Consequently, agreement is found if the forcing is additive white noise. The second and third terms on the rhs of equation (33f) are extra terms that can be associated with the Itô-to-Stratonovich correction [Pav14]. Indeed, we see this by regarding equation (31) in the Stratonovich sense:

$$\begin{equation}\mathrm{d}\mathbf{x}(t)=\mathbf{F}(\mathbf{x})\mathrm{d}t+\varepsilon \tilde{{\Gamma}}H{\circ}\mathrm{d}{W}_{t}.\end{equation} \tag{ 34 }$$

The corresponding Itô conversion reads and expands as

$$\begin{equation}\mathrm{d}\mathbf{x}(t)=\left[\mathbf{F}(\mathbf{x})+\frac{{\varepsilon }^{2}}{2}\left[\nabla \cdot \left((\tilde{{\Gamma}}H){(\tilde{{\Gamma}}H)}^{\top }\right)-\left(\tilde{{\Gamma}}H\right)\nabla \cdot \left(\tilde{{\Gamma}}H\right)\right]\right]\mathrm{d}t+\varepsilon \tilde{{\Gamma}}H\mathrm{d}{W}_{t}\end{equation} \tag{ 35a }$$

$$\begin{equation}=\left[\mathbf{F}(\mathbf{x})+\frac{{\varepsilon }^{2}}{2}\left[\nabla \cdot \left({{\Gamma}}_{:,1}{{\Gamma}}_{:,1}^{\top }\right)-\left(\nabla \cdot {{\Gamma}}_{:,1}\right){{\Gamma}}_{:,1}\right]\right]\mathrm{d}t+\varepsilon \tilde{{\Gamma}}H\mathrm{d}{W}_{t}.\end{equation} \tag{ 35b }$$

Note that an ɛ² term has appeared in the drift component. Consequently, the backward-Kolmogorov equation associated with equation (35) will only have perturbation operators of order ɛ². Following the expansions that lead to equation (24), the perturbation operator ${\mathcal{L}}_{2}^{\ast }$ for equation (35) writes as:

$$\begin{equation}{\mathcal{L}}_{2}^{\ast }f=\frac{1}{2}\left[\nabla \cdot \left({{\Gamma}}_{:,1}{{\Gamma}}_{:,1}^{\top }\right)-\left(\nabla \cdot {{\Gamma}}_{:,1}\right){{\Gamma}}_{:,1}\right]\cdot \nabla f+\frac{1}{2}{{\Gamma}}_{:,1}{{\Gamma}}_{:,1}^{\top }:{\nabla }^{2}f\end{equation} \tag{ 36a }$$

$$\begin{equation}=\frac{1}{2}\left[\left({{\Gamma}}_{:,1}\cdot \nabla \right){{\Gamma}}_{:,1}\right]\cdot \nabla f+\frac{1}{2}{{\Gamma}}_{:,1}{{\Gamma}}_{:,1}^{\top }:{\nabla }^{2}f\end{equation} \tag{ 36b }$$

$$\begin{equation}=\frac{1}{2}\left[\frac{1}{2}\nabla \left({{\Gamma}}_{:,1}\cdot {{\Gamma}}_{:,1}\right)-{{\Gamma}}_{:,1}\times \nabla \times {{\Gamma}}_{:,1}\right]\cdot \nabla f+\frac{1}{2}{{\Gamma}}_{:,1}{{\Gamma}}_{:,1}^{\top }:{\nabla }^{2}f.\end{equation} \tag{ 36c }$$

Hence, comparing equations (36c) and (33f), we are left with showing that $\nabla \left({{\Gamma}}_{:,1}\cdot {{\Gamma}}_{:,1}\right)\cdot \nabla f=2{{\Gamma}}_{:,1}\cdot \left(\nabla f\cdot \nabla \right){{\Gamma}}_{:,1}$, which can be seen by direct evaluation:

$$\begin{equation}\frac{1}{4}\nabla \left({{\Gamma}}_{:,1}\cdot {{\Gamma}}_{:,1}\right)\cdot \nabla f=\frac{1}{4}\left[\begin{matrix}\hfill 2{{\Gamma}}_{1,1}{\partial }_{{x}_{1}}{{\Gamma}}_{1,1}+\cdots +2{{\Gamma}}_{d,1}{\partial }_{{x}_{1}}{{\Gamma}}_{d,1}\hfill \\ \hfill {\vdots}\hfill \\ \hfill 2{{\Gamma}}_{1,1}{\partial }_{{x}_{d}}{{\Gamma}}_{1,1}+\cdots +2{{\Gamma}}_{d,1}{\partial }_{{x}_{d}}{{\Gamma}}_{d,1}\hfill \end{matrix}\right]\cdot \nabla f\end{equation} \tag{ 37a }$$

$$\begin{equation}=\frac{1}{2}\left[{{\Gamma}}_{:,1}\cdot {\partial }_{{x}_{1}}f{\partial }_{{x}_{1}}{{\Gamma}}_{:,1}+\cdots +{{\Gamma}}_{:,1}\cdot {\partial }_{{x}_{d}}f{\partial }_{{x}_{d}}{{\Gamma}}_{:,1}\right]\end{equation} \tag{ 37b }$$

$$\begin{equation}=\frac{1}{2}{{\Gamma}}_{:,1}\cdot \left(\nabla f\cdot \nabla \right){{\Gamma}}_{:,1}.\end{equation} \tag{ 37c }$$

Which proves the claim. The general case is formalised in the following proposition:

Proposition 3.1. Consider the SDE (21) in Stratonovich form and let $\mathbf{H}:{\mathbb{R}}^{d}\to {\mathbb{R}}^{d}$ be the Stratonovich-to-Itô correction:

$$\begin{equation}\mathbf{H}=\frac{1}{2}\left[\nabla \cdot \left({\Gamma}{{\Gamma}}^{\top }\right)-{\Gamma}\left(\nabla \cdot {{\Gamma}}^{\top }\right)\right].\end{equation} \tag{ 38 }$$

Then,

$$\begin{equation}{\tilde{\delta }}^{(2)}\left[{\Psi}\right]-{\delta }^{(2)}\left[{\Psi}\right]={\varepsilon }^{2}\int {\rho }_{0}(\mathrm{d}\mathbf{x}){\int }_{-\infty }^{\infty }{\Theta}(s)\mathbf{H}(\mathbf{x})\cdot \nabla {\Psi}(\mathbf{x}(s))\mathrm{d}s,\end{equation} \tag{ 39 }$$

where the first and second terms of the left-hand side are defined in equations (27) and (28), respectively.

Proof. We note the following chain of equalities:

$$\begin{equation}\mathrm{d}\mathbf{x}=\mathbf{F}\mathrm{d}t+\varepsilon {\Gamma}{\circ}\mathrm{d}{W}_{t}=\left[\mathbf{F}+\mathbf{H}\right]\mathrm{d}t+\varepsilon {\Gamma}\mathrm{d}{W}_{t}\end{equation} \tag{ 40a }$$

$$\begin{equation}=\left[\mathbf{F}+{\varepsilon }^{2}\sum\limits _{k=1}^{p}\,{\mathbf{H}}_{k}\right]\mathrm{d}t+\varepsilon \sum\limits _{i=1}^{p}\,{\tilde{{\Gamma}}}_{i}{H}_{i}\mathrm{d}{W}_{t},\end{equation} \tag{ 40b }$$

where ${\tilde{{\Gamma}}}_{i}=\mathrm{diag}\left({{\Gamma}}_{:,i}\right)$, ${H}_{i}={\mathbb{1}}_{d}{\mathbf{e}}_{i}^{\top }$ and H_k is defined as:

$$\begin{equation}{\mathbf{H}}_{k}=\frac{1}{2}\left[\nabla \cdot \left({{\Gamma}}_{:,k}{{\Gamma}}_{:,k}^{\top }\right)-{\tilde{{\Gamma}}}_{k}\nabla \cdot {\tilde{{\Gamma}}}_{k}\right].\end{equation} \tag{ 41 }$$

Now, it is enough to repeat the argument started in equation (31) for each ${\tilde{{\Gamma}}}_{i}$ and then use the linearity of the leading order response to extend it for i = 1, ..., p. □

In [LW17] it was shown how to compute the change in the correlation properties of a deterministic chaotic system resulting from a static forcing applied to its dynamics. This amounts to considering the case described in equation (1) with Σ = 0 and with g(t) = 1. Indicating with ${C}_{{\Psi},{\Phi}}^{\varepsilon }$ the correlation function as two square-integrable observables Ψ and Φ in the perturbed dynamics, one gets:

$$\begin{equation}{C}_{{\Psi},{\Phi}}^{\varepsilon }(t)=\int {\rho }_{\varepsilon }(\mathrm{d}\mathbf{x}){\mathrm{e}}^{t({\mathcal{L}}_{0}^{\ast }+\varepsilon {\mathcal{L}}_{1}^{\ast })}{\Psi}(\mathbf{x}){\Phi}(\mathbf{x})\end{equation} \tag{ 42a }$$

$$\begin{equation}=\int {\rho }_{0}(\mathrm{d}\mathbf{x}){\mathrm{e}}^{t{\mathcal{L}}_{0}^{\ast }}{\Psi}(\mathbf{x}){\Phi}(\mathbf{x})\end{equation} \tag{ 42b }$$

$$\begin{equation}\quad +\varepsilon \int {\rho }_{0}(\mathrm{d}\mathbf{x}){\int }_{0}^{\infty }{\mathcal{L}}_{1}^{\ast }{\mathrm{e}}^{(t+s){\mathcal{L}}_{0}^{\ast }}{\Psi}(\mathbf{x}){\mathrm{e}}^{s{\mathcal{L}}_{0}^{\ast }}{\Phi}(\mathbf{x})\mathrm{d}s\end{equation} \tag{ 42c }$$

$$\begin{equation}\quad +\varepsilon \int {\rho }_{0}(\mathrm{d}\mathbf{x}){\int }_{0}^{t}{\mathrm{e}}^{(t-s){\mathcal{L}}_{0}^{\ast }}{\mathcal{L}}_{1}^{\ast }{\mathrm{e}}^{s{\mathcal{L}}_{0}^{\ast }}{\Psi}(\mathbf{x}){\Phi}(\mathbf{x})\mathrm{d}s+\mathcal{O}\left({\varepsilon }^{2}\right).\end{equation} \tag{ 42d }$$

where ${\mathcal{L}}_{0}^{\ast }$ and ${\mathcal{L}}_{1}^{\ast }$ are given in equations (23a) and (9b), respectively. The term shown in equation (42b) gives the unperturbed correlation. The leading order correction to correlation functions consists of two different components corresponding to the perturbation expansion of the invariant measure and the Koopman operator, respectively. The term (42c) can also be written as

$$\begin{equation}\varepsilon \int {\rho }_{1}(\mathrm{d}\mathbf{x}){\int }_{0}^{\infty }{\mathrm{e}}^{t{\mathcal{L}}_{0}^{\ast }}{\Psi}(\mathbf{x}){\Phi}(\mathbf{x}),\end{equation} \tag{ 43 }$$

which shows that it can be interpreted as the expectation value of the unperturbed form of time-delayed product of the two observables in the first order correction of the measure. Instead, the first order contribution given in equation (42d) is qualitatively different and vanishes if t = 0, because in such a case the regular response theory for observables applies.

We can generalise the previous result in two different directions. If we stick to deterministic background dynamics (so that ${\mathcal{L}}_{0}^{\ast }$ is given by equation (23a)) but consider, instead, adding stochastic forcing to the system—as in equation (21)—the change in the correlation properties of the system can be evaluated by performing the following substitutions in equations (42) and (43): ${\mathcal{L}}_{1}^{\ast }\to {\mathcal{L}}_{2}^{\ast }$, where the latter operator is defined in equation (23b); ρ₁ → ρ₂, where the second-order correction to the invariant measure is defined in equation (5); and ɛ → ɛ². If, instead, we consider a stochastic reference dynamics (so that ${\mathcal{L}}_{0}^{\ast }$ is given by equation (9a)) and want to study the impact of changing the noise law, we can insert the operator ${\mathcal{L}}_{1}^{\ast }$ given in equation (13) into equation (42). Of course equation (42) apply seamlessly in the case we alter the drift term of the SDE.

Applying the Fourier transform to ${C}_{{\Psi},{\Phi}}^{\varepsilon }(t)$ yields the response in spectral domain, again following [LW17]. We can, however, use the resolvent formalism and Laplace transforms of semigroups to express the correlation function in spectral domain [EN00]. The Laplace transform of a correlation function reads as:

$$\begin{equation}\mathfrak{L}[{C}_{{\Psi},{\Phi}}^{\varepsilon }](z)=\int {\rho }_{\varepsilon }(\mathrm{d}\mathbf{x})R(z,{\mathcal{L}}_{0}^{\ast }+\varepsilon {\mathcal{L}}_{1}^{\ast }){\Psi}(\mathbf{x}){\Phi}(\mathbf{x}),\end{equation} \tag{ 44 }$$

for $\mathfrak{R}\mathfrak{e}z > 0$ and where $R(z,{\mathcal{L}}_{0})={\left(z-{\mathcal{L}}_{0}\right)}^{-1}$ is the resolvent operator. Expanding $\mathfrak{L}[{C}_{{\Psi},{\Phi}}^{\varepsilon }](z)$ in powers of ɛ we get a linear response for correlation functions in terms of derivatives with respect to ɛ:

$$\begin{align}\hfill \frac{\mathrm{d}}{\mathrm{d}\varepsilon }\mathfrak{L}[{C}_{{\Psi},{\Phi}}^{\varepsilon }](z){\vert }_{\varepsilon =0}& =\int {\rho }_{1}(\mathrm{d}\mathbf{x})R(z,{\mathcal{L}}_{0}^{\ast }){\Psi}(\mathbf{x}){\Phi}(\mathbf{x})+\int {\rho }_{0}(\mathrm{d}\mathbf{x})R(z,{\mathcal{L}}_{0}^{\ast }){\mathcal{L}}_{1}^{\ast }R(z,{\mathcal{L}}_{0}^{\ast }){\Psi}(\mathbf{x}){\Phi}(\mathbf{x}).\hfill \end{align} \tag{ 45 }$$

This formula can be obtained by taking the Laplace transform of equation (42), although the resolvent formalism makes the calculations substantially easier. We, thus, provide an explicit link between the response in frequency domain and the spectral properties of the unperturbed generator ${\mathcal{L}}_{0}^{\ast }$ via the resolvent. Similarly to the discussion above, upon the replacement of ${\mathcal{L}}_{0}^{\ast }$ or ${\mathcal{L}}_{1}^{\ast }$, one can deduce the response formula for deterministic or stochastic perturbations.

In [Luc12] linear response of power spectral densities is related to the modulus of the susceptibility function in case of white-noise modulated forcings. Indeed, there it is argued that if $\mathfrak{F}\left[{C}_{{\Psi},{\Psi}}^{\varepsilon }\right](\omega )$ is the perturbed power spectrum of the observable Ψ then we have:

$$\begin{equation}\mathbb{E}\left[\mathfrak{F}\left[{C}_{{\Psi},{\Psi}}^{\varepsilon }\right](\omega )\right]-\mathfrak{F}\left[{C}_{{\Psi},{\Psi}}\right](\omega )\approx \mathbb{E}\left[{\left\vert {\delta }_{\sigma }^{(2)}\left[{\Psi}\right]\right\vert }^{2}\right]\approx {\varepsilon }^{2}{\left\vert \mathfrak{F}\left[\mathcal{G}\right](\omega )\right\vert }^{2}.\end{equation} \tag{ 46 }$$

In order to prove this, the author of [Luc12] invoked the Wiener–Khinchin theorem and applied it to the autocorrelation function of ${\delta }_{\sigma }^{(1)}[{\Psi}](t)$ before taking averages over σ. Equation (46) appears to have a term missing if one compares it to equation (45). This is because that the author uses $\mathfrak{F}\left[{C}_{{\Psi},{\Psi}}^{\varepsilon }\right]\approx \mathfrak{F}\left[\int {\rho }_{0}(\mathrm{d}\mathbf{x}){\Psi}(\mathbf{x}){\mathrm{e}}^{t\mathcal{L}}{\Psi}(\mathbf{x})\right]$ as an approximation. Indeed, before averaging over realisations, the leading order response of the flow is proportional to ɛ whereas the perturbed invariant measure responds proportionally to ɛ². Therefore, the rhs of equation (45) yields a single term from which equation (46) follows.

Notice that the left-hand side of the equation above is especially easy to calculate if one is let to sample the perturbed and unperturbed dynamics. Thus, it is possible to obtain estimates of the modulus of the susceptibility function in a fully empirical way. One is left with finding the actual locations of the poles, which the modulus itself cannot reveal. This would be attained by noticing that the logarithm of the susceptibility function can be written as:

$$\begin{equation}\mathrm{log}\left(\mathfrak{F}\left[\mathcal{G}\right](\omega )\right)=\mathrm{log}\left(\left\vert \mathfrak{F}\left[\mathcal{G}\right](\omega )\right\vert \right)+i\varphi (\omega ),\end{equation} \tag{ 47 }$$

for every ω and some phase function φ. Assuming, further, that $\mathrm{log}\left(\mathfrak{F}\left[\mathcal{G}\right](\omega )\right)$ is analytic, Kramers–Kronig relations are invoked to link its real and imaginary parts shown in equation (47) and, thus, that of $\mathfrak{F}\left[\mathcal{G}\right](\omega )$. This way, equation (46) provides a practical way of estimating the susceptibility function of a deterministic flow by using time series of the stochastically forced vector field.

In this note we have shown how the formalism of operator semigroups allows for a compact derivation of response formulas describing at leading order the impact of adding forcings of either deterministic or stochastic nature in both cases of systems whose background dynamics is deterministic or stochastic. Linear response gives the leading order correction to expectation values due to extra deterministic forcings acting on either stochastic or chaotic dynamical systems. The response is also linear in the intensity of the (extra) noise for background stochastic dynamics. Note that we are able to accommodate results pertaining to non-stationary noise laws.

In the case of stochastic systems, the linear response formulas have been expressed through a spectral decomposition on the basis of the eigenvectors and eigenfunctions of the unperturbed Fokker–Planck/backward-Kolmogorov generator; see equation (17). This angle gives a basis for bottom-up construction of response functions and relaxation rates for non-equilibrium systems, and clarifies the general applicability and interpretation of the FDT. In particular, it underpins the close link between the presence of (a) a vanishing spectral gap in the system of interest, (b) of a strongly amplified response to forcings, and (c) of a slow decay of correlations (sometimes referred to as critical slowing down) between observables, which acts as early warning signal indicating the proximity of a tipping point [TLD18, TLLD18, GL20, SBBBC09, LHKHL08, Kue11]. It also casts the computation of the response of an observable to perturbations in the context of the problem of finding suitable and accurate estimates for the modes of the Koopman operator [Mez05, BMM12]. This research area has benefitted considerably from the use of data-driven methods [KBBP16, LKB18]. Koopman modes are ordered according to their decay rate, thus giving a rationale for, in most cases, the greater relevance of the contributions coming from slowly decaying modes; see equations (17) and (18). We remark that such a connection between response formulas and Koopmanism applies under the conditions that the residual spectrum decays on considerably faster time-scales than the point spectrum and/or that we consider observables whose projection on the modes associated with such residual spectrum is negligible. An instance where this framework holds can be found in [WILA16, SGLCG21, section D]. Broadly, this indicates that our approach could be suitable for the cases where model reduction techniques based on the Koopman operator are applicable [BMM12, TRL+14, LKB18].

Future endeavours should aim at finding computationally amenable ways of identifying the eigenvalues and the associated timescales appearing in equation (17). We remark that this problem is attracting considerable interest in the climate community, because it would allow to draw a convincing link between climate variability and climate response [GL20]. Nonetheless, recent proposals in this direction—while indeed interesting and promising [TMPR21a, TMPR21b, BDH21]— are yet unable to deal with complex eigenvalues (see at this regard the methodology discussed in reference [103]), thus missing out an essential ingredient of the dynamics of nonequilibrium systems [SGLCG21]. Empirically-learnt Markov-matrices are a promising technique that should be explored, even taking into account the need to necessarily operate in a much reduced dimension phase space for any application where the dynamics is high-dimensional [CNK+14, CDN+20]. Note that in the case of Markov chains it is easy to estimate the actual radius of convergence of the response operators, i.e. the maximum intensity of the applied forcing [Luc16, SGL20]. Contrarily, we highlight the results of [Ken71, Ken73] which present a non-perturbative response formula based on the Mori–Zwanzig projection operator [Zwa01]; this formula is revisited in appendix A. In this paper we have on purpose avoided to discuss any aspect of convergence of the perturbative formulas presented above; this is left for future investigations.

When stochastically perturbing systems whose underlying dynamics is chaotic, the second order-response gives the leading-order correction. We have been able to reconcile two results published in the literature and showed that the formulas given in [Luc12] and in [Abr17] coincide when one interprets the results of the former study taking into account the physically relevant Stratonovich interpretation of noise, which indeed complies with the heuristic procedure followed there, where the standard rules of calculus had been used throughout the paper. This result helps to better understand how the invariant measure convergences to the physical one when the intensity of the noise becomes infinitesimal.

Finally, we have adapted the mathematical technology developed for studying the change in the expectation value of observables due to the application of external forcings combined with the perturbative expansion of the resolvent operator [EN00] to generalise the results presented in [LW17]. We have developed leading response formulas for the two-point correlations of either deterministic or stochastic dynamical systems. Also in this case, it is easy to show that the order of correction is linear in the case of deterministic perturbation to either stochastic or chaotic systems and of stochastic perturbations of stochastic systems, while it is quadratic in the case of purely stochastic perturbations applied to chaotic systems. The weak noise assumption is fundamental to derive equation (46), originally presented and discussed with the help of numerical simulations in [Luc12]. A more detailed study of this relation would be most valuable since its implementation is straightforward and could provide new operational algorithms for computing the spectrum of the Fokker–Planck operators and the subsequent susceptibility functions, which are key to understand critical transitions in complex systems [Luc18, TLD18, TLLD18].

On a more general note, the effects of stochastic noise on the forced variability of a deterministic system are yet to be fully understood. For instance, numerical studies show that under the weak-noise assumption—potentially at the linear response level, as in section 3—the low-frequency variability spectrum of a multi-modal deterministic system can be heavily amplified as results of an apparent stabilisation of the different coexisting regimes [Kwa14]. While the counterintuitive fact that noise stabilise an otherwise unstable deterministic system is well known, it is clear that the underlying dynamical structure of the deterministic system is also essential for such a behaviour [SW07, TGOS08, PKM11, CP17, MS18, KCSW19, MSdSB22]. Indeed, it has been argued that the systems displaying intermittency or activation are more likely to be stabilised under the noisy perturbations [Dor21]. The present paper does not tackle this question, although we here provide formal methods to treat systems subject to general external forcings and it is hoped that our results will help understand the interplay between deterministic and stochastic feature of the dynamics of rather general class of systems featuring non-trivial variability.

The authors would like to thank Mickaël Chekroun, Joshua Dorrington, Greg Pavliotis, and Jeroen Wouters, for their encouragement, comments, and suggestions. VL acknowledges the support received from the Horizon 2020 project TiPES (Grant No. 820970), from the Marie Curie ITN CriticalEarth (Grant No. 956170), and from the EPSRC project EP/T018178/1.

No new data were created or analysed in this study.

The limit of weak forcing strength is instrumental in deriving meaningful linear response formulas that capture statistical changes due to an applied external field. When small, yet finite-sized perturbations are introduced, power series are resorted to and their convergence depends on a radius of expansion. To obtain the full response to any order in ɛ, one seems inevitably obliged to solve the perturbed Fokker–Planck equation (2) to find the probability state at a certain time to evaluate the phase average a given observable. In the present section, we aim at deriving a scalar equation for the response function ⟨Ψ, ρ(⋅, t)⟩ in equation (7) that does not require a power expansion nor the full knowledge of perturbed probability density function. To do so, we revisit the calculations presented in [Ken71, Ken73] in the context of quantum mechanical systems, at a time when linear response and the subsequent fluctuation–dissipation theorem were just being introduced. Here, our goal is to review such results but in a stochastic framework and with references to recent advances in the treatment of operator equations. To do this, the equation of motion—regardless it being stochastic or deterministic—will suitably be projected to derive, with the aid of the Mori–Zwanzig formalism [Zwa01], a scalar equation for the response ⟨Ψ, ρ(⋅, t)⟩ where the only independent variable is ⟨Ψ, ρ(⋅, t)⟩ itself; this is, an homogeneous equation.

To achieve this goal, we shall start by considering the forced SDE (1). Since the introduction of perturbations makes the equations of motion non-autonomous, the usual exponential operator no longer serves as a notation for the solution of equation (2) unless we consider time-ordered exponentials [Gil17]. It is useful to define the full response function R_Ψ(t) as:

$$\begin{equation}{R}_{{\Psi}}(t){:=}\left\langle {\Psi},\rho (\cdot ,t)\right\rangle =\int {\Psi}(\mathbf{x})\rho (\mathbf{x},t)\mathrm{d}\mathbf{x},\end{equation} \tag{ 48 }$$

where ρ(⋅, t) is a solution of equation (2). Traditionally, the response function is defined as the difference between the perturbed and unperturbed means of the observable Ψ, but since the latter is stationary, this definition will not introduce differences with other formulations.

The idea is now to project the Fokker–Planck equation according to $\mathbb{P}f=\left(\int {\Psi}f\right){\rho }_{0}$, so that the Mori–Zwanzig formalism yields an equation in terms of projected operators instead of the full ones. Upon substitution and simplification we find that the evolution equation for the response function R_Ψ can be written in more compact terms as:

$$\begin{equation}{\partial }_{t}{R}_{{\Psi}}(t)=B(t){R}_{{\Psi}}(t)+{\int }_{0}^{t}K(t,s){R}_{{\Psi}}(s)\mathrm{d}s,\end{equation} \tag{ 49 }$$

where the newly introduced function B(t) and memory kernel K(t, s) are defined as:

$$\begin{equation}B(t)={\xi }^{-1}\int {\Psi}(\mathbf{x})\mathcal{L}(t){\rho }_{0}(\mathbf{x})\mathrm{d}\mathbf{x},\end{equation} \tag{ 50a }$$

$$\begin{equation}K(t,s)={\xi }^{-1}\int {\Psi}(\mathbf{x})\mathcal{L}(t){\mathrm{e}}^{{\int }_{s}^{t}\left(1-\mathbb{P}\right)\mathcal{L}(\tau )\mathrm{d}\tau }\left(1-\mathbb{P}\right)\mathcal{L}(s){\rho }_{0}(\mathbf{x})\mathrm{d}\mathbf{x},\end{equation} \tag{ 50b }$$

where $\xi =\left\langle {\Psi},{\rho }_{0}\right\rangle \ne 0$. Equation (49) becomes a scalar and homogeneous equation on R_Ψ, where there is not explicit reference to the time-dependent measure ρ(x, t) that would be obtained by solving the full Fokker–Planck equation (2). Collecting the leading order terms in equation (49) one obtains the usual linear response, as shown in [Ken71].

Solving the non-Markovian equation (49) would yield the response function and, in principle, should be equal to the perturbative approach in equation (7), provided that the parameter ɛ is within the radius of convergence. In the work of [Ken71], there is no proof of the convergence of equation (49), which would boil down to determining the decay of K(t) as t tends to infinity. If one assumes that g is constant and that $\mathcal{L}$ generates a quasi-compact semigroup—see the beginning of section 2—, one can view $\mathcal{L}\mathbb{P}$ as a compact perturbation of $\mathcal{L}$, in which case it would generate a quasi-compact semigroup [EN00], meaning that the decay is dominated by a sum of exponential functions, cf equation (17). In case K(t, s) = K(t − s, 0), equation (49) can be solved explicitly in spectral domain and one could relate K(t − s, 0) with the high-order responses. Finding further simplifications of the memory kernel K(t, s) would be of great interest, since it would provide a means of computing full responses R_Ψ without resorting to perturbation expansions, nor solving the perturbed Fokker–Planck equation.