Generalized second fluctuation-dissipation theorem in the nonequilibrium steady state: Theory and applications

Yuanran Zhu yzhu56@ucmerced.edu Huan Lei Changho Kim Department of Applied Mathematics, University of California, Merced, CA 95343 Department of Computational Mathematics, Science & Engineering and Department of Statistics & Probability, Michigan State University, MI 48824

Abstract

In this paper, we derive a generalized second fluctuation-dissipation theorem (FDT) for stochastic dynamical systems in the steady state. The established theory is built upon the Mori-type generalized Langevin equation for stochastic dynamical systems and only uses the properties of the Kolmogorov operator. The new second FDT expresses the memory kernel of the generalized Langevin equation as the correlation function of the fluctuation force plus an additional term. In particular, we show that for nonequilibrium states such as heat transport between two thermostats with different temperatures, the classical second FDT is valid even when the exact form of the steady state distribution is unknown. The obtained theoretical results enable us to construct a data-driven nanoscale fluctuating heat conduction model based on the second FDT. We numerically verify that the new model of heat transfer yields better predictions than the Green-Kubo formula for systems far from the equilibrium.

^†^†journal: arXiv

1 Introduction

The fluctuation-dissipation relations are one of the most important results in statistical mechanics. In 1966, Kubo, in his renowned paper [20], proposed a general relationship that connects the response of a given system to the external perturbation and the internal thermal fluctuation of the system in the absence of disturbance. This groundbreaking result, together with many generalizations for open systems, is now categorized as the fluctuation-dissipation theorems (FDTs) which provide a far reaching generalization of Einstein’s theory of Brownian motion and Nyquist’s work on thermal noise in electrical resistors. In his original paper [20], Kubo presented two results which characterize the fluctuation-dissipation relationship for arbitrary Hamiltonian systems. The aforementioned response result is called the first FDT, whereas another relationship which builds connections between the thermal fluctuation of an observable and the memory kernel for the corresponding generalized Langevin equation (GLE) is now known as the second FDT. It is not widely known that Kubo actually provided two derivation methods for these two relationships. The first method is based on a phenomenological linear GLE for a system observable $u(t)$ :

\displaystyle\frac{d}{dt}u(t)=-\int_{0}^{t}K(t-s)u(s)ds+\frac{1}{m}f(t)+\frac{1}{m}g(t),

(1)

where $K(t)$ is known as the memory kernel and $f(t)$ is the fluctuation force satisfying $\langle f(t)\rangle=0$ and $\langle u(s)f(t)\rangle=0$ for $t>s$ . Here $\langle\cdot\rangle$ is the ensemble average in the thermal equilibrium. The system is assumed to be perturbed by a periodic external force $g(t)=g_{0}\cos(\omega t)$ . For such a case, the first FDT is proved using the standard Fourier-Laplace transform and the second FDT is derived from the first FDT hence can be viewed as a corollary of it [20]. The second method differs for the derivation of these two FDTs. In particular, the first FDT is derived using a perturbation method which is known as the linear response theory. During the derivation, the phenomenological GLE (1) is not used. The second FDT, however, can be derived using Mori’s equation [33], which can be interpreted as the (1) in an operator form. In fact, since the evolution operator for the Hamiltonian system is given by $e^{t\mathcal{L}}$ , the second FDT is a natural result of the skew-symmetry of the Liouville operator $\mathcal{L}$ with respect to the Hilbert space inner product $\langle\cdot,\cdot\rangle_{\rho_{eq}}$ , where $\rho_{eq}$ is the equilibrium distribution. Detailed explanations can be found in Section 3.

Since the establishment of the FDTs in the 1960s, there has been a considerable amount of work on its verification, violation, generalization and applications. In fact, the FDT-related linear response relation has became the cornerstone of nonequilibrium statistical mechanics, which provides a powerful tool to study various transport phenomena for near-equilibrium systems [43, 35]. Here we only mention some recent studies in active matter [3], turbulence [30, 12] and heat conduction [22, 26, 7]. We note that most of these studies are about the first FDT and relatively less attention has been paid to the second FDT. For most applications, the second FDT was required for the construction of the GLE when it was used as a reduced-order model for particular observables. Its generalization and verification has not yet been widely studied except for Mae’s recent work [28]. One of the reasons why this is the case is that, unlike the first FDT which is directly related to the Green-Kubo formula for transport theory, the study on the second FDT is difficult to translate into useful, experimentally verifiable results. In fact, Kubo’s second derivation method reviewed in the previous paragraph clearly indicates the first and the second FDTs are actually quite different. Specifically, the second FDT is an intrinsic property for Hamiltonian systems in the thermodynamic equilibrium and its validity does not rely on perturbation arguments. This observation hints that one can possibly derive the second FDT or its generalized form from the operator-form GLE, i.e. the Mori-Zwanzig equation [33, 51], for generalized stochastic dynamics in the nonequilibrium steady state. In addition, noting a gradually increasing interest in open systems and far-from-equilibrium phenomena, we anticipate that a classical or generalized second FDT for the nonequilibrium system will provide insights on the system dynamics. All these thoughts motivated the current research.

The main purpose of this paper is two-fold. First, we will follow Kubo-Mori’s methodology to establish a generalized second FDT for stochastic dynamical systems driven by Gaussian white noise. Such stochastic models are widely used in coarse-graining modeling for molecules [10, 11, 24, 17, 13], nonequilibrium heat conduction models [9, 26, 22] and many other open systems [30, 29]. Secondly, the second FDT is shown to be valid in the nonequilibrium steady state, therefore can be applied to address the far-from-equilibrium transport problem.

This article is organized as follows. Section 2 provides a short summary of the theoretical results obtained in this paper. Section 3 briefly reviews the derivation of the GLE for stochastic systems using the Mori-Zwanzig equation. In Section 4, as a comparison, we first review the derivation of the first FDT for stochastic systems via perturbation analysis, then we use the GLE to derive a generalized second FDT which holds in the nonequilibrium steady state. In Section 5, this newly established relation is applied to equilibrium and nonequilibrium systems commonly used in the statistical mechanics. In particular, we show the validity of the classical second FDT for the averaged heat flux in a heat conduction model. The theoretical results are verified numerically in Section 6 via the molecular dynamics simulation of the heat conduction model. In addition, two reduced-order models are proposed and shown to faithfully characterize the thermal conduction for far-from-equilibrium systems. The main findings of the paper are summarized in Section 7.

2 The main theoretical results

In this section, we give a short summary of the main theoretical results obtained in this paper and provide simple examples to explain them. In statistical physics, open systems normally refer to dynamical systems in contact with thermal reservoirs. These reservoirs are modeled by deterministic or stochastic thermostats. If a nonequilibrium condition is imposed, for instance one may choose thermostats with different temperatures and attach them to a Hamiltonian system, then the external force exerted by the thermostats will drive the system out of the thermal equilibrium. In this paper, we are mainly concerned with the nonequilibrium systems with stochastic thermostats. Mathematically, they can be described by a general stochastic differential equation (SDE):

\displaystyle d\bm{x}(t)=\bm{F}(\bm{x}(t))+\sigma(\bm{x}(t))d\bm{\mathcal{W}}(t),\qquad\bm{x}(0)=\bm{x}_{0}\sim\rho_{0}(\bm{x}),

(2)

Stochastic analysis already told us [31] after a finite transient time, SDE (2) satisfying some suitable conditions will converge to a unique steady state. For stochastic force $\sigma(\bm{x}(t))d\bm{\mathcal{W}}(t)$ in a nonequilibrium setting, such a steady state is normally termed as the nonequilibirum steady state (NESS). The primary goal of the current paper is to investigate whether the classical second FDT still holds in the NESS of SDEs. To this end, we follow the aforementioned Mori-Kubo’s methodology and use the following GLE for the (open) stochastic system as the starting point of our analysis:

\displaystyle\frac{d}{dt}u(t)=\Omega u(t)+\int_{0}^{t}K(t-s)u(s)ds+f(t),

(3)

where $u(t)=e^{t\mathcal{K}}u(0)$ is an arbitrary observable function of the stochastic system. Using the properties of the Kolmogorov backward operator $\mathcal{K}$ (generator of the stochastic dynamics), we proved the following result regarding the validity of the second FDT:

(I) In GLE (3), the memory kernel $K(t)$ and the fluctuation force $f(t)$ satisfy a generalized second FDT:

\displaystyle K(t)=-\frac{\langle f(t),f(0)\rangle_{\rho}}{\langle u^{2}(0)\rangle_{\rho}}+\frac{\langle f(t),w(0)\rangle_{\rho}}{\langle u^{2}(0)\rangle_{\rho}},

(4)

where $\langle\cdot\rangle_{\rho}$ is the ensemble average in the NESS. In (4), the first term yields the classical second FDT. The additional term $\langle f(t),w(0)\rangle_{\rho}/\langle u^{2}(0)\rangle_{\rho}$ constitutes the generalized relation that the current paper studies. However, since $w(0)$ in (4) depends on the steady state probability density $\rho$ (see its explicit expression (5)), which is hard to obtain for most high-dimensional stochastic systems, this generalized relation cannot be used directly. To overcome this difficulty, we notice an important fact about the obtained second FDT (4) which is the second main result of this paper:

(II) The additional observable $w(0)$ can be explicitly written as

\displaystyle w(0)=\mathcal{S}u(0)=2\sigma^{2}\Delta u(0)+2\sigma^{2}\nabla(\ln\rho)\cdot\nabla,

(5)

where $\mathcal{S}$ is a second-order differential operator in the non-degenerate coordinate of the stochastic system. If in addition, observable $u(0)$ is a function of the degenerate coordinate, then $w(0)=\mathcal{S}u(0)=0$ and the generalized second FDT (4) degenerates to the classical second FDT.

We can use the Langevin dynamics as an example to explain the meaning of the degenerate coordinate and result (II). As it is well-known, the Gaussian white noise for a Langvein dynamics is only imposed in the momentum coordinate $p_{i}$ . We shall call the momentum coordinate $p_{i}$ non-degenerate and the position coordinate $q_{i}$ degenerate. Bearing this in mind, the reason why the special case discussed in (II) has $w(0)=\mathcal{S}u(0)=0$ becomes obvious. Since the initial condition of the observable $u(0)=f(q_{i}(0))$ , then we have

\displaystyle w(0)=\mathcal{S}u(0)=\mathcal{S}f(q_{i}(0))=2\sigma^{2}\partial^{2}_{p_{i}}f(q_{i}(0))+2\sigma^{2}\partial_{p_{i}}(\ln\rho)\partial_{p_{i}}f(q_{i}(0))=0.

This result can be physically interpreted as follows: For the linear GLE of an arbitrary (open) stochastic system, the validity of the classical second FDT of an observable $u(t)$ depends on whether $u(t)$ has direct interactions with the environment through the white noise. In reality, many nonequilibrium systems can be modeled by highly degenerate stochastic systems. Hence most observables in such nonequilibrium systems satisfy the classical second FDT. We will use this fact and the GLE (3) to build effective reduced-order models for low-dimensional observables. In Section 6, we show the application of such an effective model in studying the far-from equilibrium transport phenomena, which are generally hard to handle using the linear response theory.

3 Mori-type generalized Langevin equations for SDEs

The starting point of our analysis is the Mori-Zwanzig (MZ) equation for stochastic differential equations (SDEs). Such stochastic system MZ equation was derived independently by different researchers using various methods [10, 34, 49]. Here we briefly review the derivation used in [49]; more detailed explanations can be found therein. Let us consider a $d$ -dimensional SDE on a smooth manifold $M$

\displaystyle d\bm{x}(t)=\bm{F}(\bm{x}(t))dt+\bm{\sigma}(\bm{x}(t))d\bm{\mathcal{W}}(t),\qquad\bm{x}(0)=\bm{x}_{0}\sim\rho_{0}(\bm{x}),

(6)

where $\bm{F}:M\rightarrow\mathbb{R}^{d}$ and $\bm{\sigma}:M\rightarrow\mathbb{R}^{d\times m}$ are smooth functions, $\bm{\mathcal{W}}(t)$ is an $m$ -dimensional Wiener process with independent components, and $\bm{x}_{0}$ is a random initial state characterized in terms of a probability density function $\rho_{0}(\bm{x})$ . For the deterministic initial condition we have $\rho_{0}(\bm{x})=\prod_{i=1}^{d}\delta(x_{i}-x_{i}(0))$ . The solution of (6) with different initial conditions form a $d$ -dimensional stochastic flow on the manifold $M$ , which is known as Brownian flow [23]. The Markov property of (6) allows us to define a composition operator $\mathcal{M}(t,0)$ that pushes forward in time the average of the observable $\bm{u}(t)=\bm{u}(\bm{x}(t))$ over the noise, i.e.,

\displaystyle\mathbb{E}_{\bm{\mathcal{W}}(t)}[\bm{u}(\bm{x}(t))|\bm{x}_{0}]=\mathcal{M}(t,0)\bm{u}(\bm{x}_{0})=e^{t\mathcal{K}}\bm{u}(\bm{x}_{0}).

(7)

Using Itô’s interpretation of the white noise, $\mathcal{M}(t,0)$ is known as a Markovian semigroup generated by the backward Kolmogorov operator [39, 19]:

\displaystyle\mathcal{K}(\bm{x}_{0})

\displaystyle=\sum_{k=1}^{d}F_{k}(\bm{x}_{0})\frac{\partial}{\partial x_{0k}}+\frac{1}{2}\sum_{j=1}^{m}\sum_{i,k=1}^{d}\sigma_{ij}(\bm{x}_{0})\sigma_{kj}(\bm{x}_{0})\frac{\partial^{2}}{\partial x_{0i}\partial x_{0k}}.

(8)

From now on, with a slight abuse of notation, we will use $\bm{u}(\bm{x}(t))$ to represent its noise-averaged quantity (7). In fact, without any further specification, all $\bm{u}(\bm{x}(t))$ that appear in the rest of this paper represent the noise-averaged quantity (7). We also note (7) is only a conditional expectation. When the initial condition is random, (7) is still a stochastic quantity with initial randomness. We now introduce a projection operator $\mathcal{P}$ and its orthogonal operator $\mathcal{Q}=\mathcal{I}-\mathcal{P}$ . By differentiating the well-known Dyson’s identity

\displaystyle e^{t\mathcal{K}}=e^{t\mathcal{Q}\mathcal{K}}+\int_{0}^{t}e^{s\mathcal{K}}\mathcal{P}\mathcal{K}e^{(t-s)\mathcal{Q}\mathcal{K}}\mathcal{Q}\mathcal{K}ds,

it is straightforward to obtain the following MZ equation governing the evolution of the noise-averaged observable (7):

\frac{\partial}{\partial t}e^{t\mathcal{K}}\bm{u}(0)=e^{t\mathcal{K}}\mathcal{P}\mathcal{K}\bm{u}(0)+e^{t\mathcal{Q}\mathcal{K}}\mathcal{Q}\mathcal{K}\bm{u}(0)+\int_{0}^{t}e^{s\mathcal{K}}\mathcal{P}\mathcal{K}e^{(t-s)\mathcal{Q}\mathcal{K}}\mathcal{Q}\mathcal{K}\bm{u}(0)ds.

(9)

The three terms at the right hand side of (9) are called, respectively, streaming term, fluctuation (or noise) term, and memory term. Applying the projection operator $\mathcal{P}$ to (9) yields the projected MZ equation

\frac{\partial}{\partial t}\mathcal{P}e^{t\mathcal{K}}\bm{u}(0)=\mathcal{P}e^{t\mathcal{K}}\mathcal{PK}\bm{u}(0)+\int_{0}^{t}\mathcal{P}e^{s\mathcal{K}}\mathcal{PK}e^{(t-s)\mathcal{Q}\mathcal{K}}\mathcal{QK}\bm{u}(0)ds.

(10)

Note that the MZ equation (9) and its projected form (10) for stochastic systems have the same structure as the classical MZ equations for deterministic (autonomous) systems [49, 46, 48]. However, the Liouville operator $\mathcal{L}$ is now replaced by the Kolmogorov operator $\mathcal{K}$ . Let us consider the weighted Hilbert space $H=L^{2}(M,\rho)$ with inner product

\langle h,g\rangle_{\rho}=\int_{M}h(\bm{x})g(\bm{x})\rho(\bm{x})d\bm{x},\qquad h,g\in H,

(11)

where $\rho$ is a positive weight function on $M$ which is often chosen to be a certain type of probability densities. With this Hilbert space, we can introduce the Mori-type projection operator:

\displaystyle\mathcal{P}h=\sum_{i,j=1}^{N}G^{-1}_{ij}\langle u_{i}(0),h\rangle_{\rho}u_{j}(0),\qquad h\in H,

(12)

where $G_{ij}=\langle u_{i}(0),u_{j}(0)\rangle_{\rho}$ and $u_{i}(0)=u_{i}(\bm{x}(0))$ ( $i=1,...,N$ ) are $N$ linearly independent functions. It is easy to check that $\mathcal{P}$ and its orthogonal $\mathcal{Q}=\mathcal{I}-\mathcal{P}$ are both symmetric projection operators in $L^{2}(M,\rho)$ , i.e. $\mathcal{P}^{*}=\mathcal{P}=\mathcal{P}^{2}$ , $\mathcal{Q}^{*}=\mathcal{Q}=\mathcal{Q}^{2}$ , where $\mathcal{P}^{*},\mathcal{Q}^{*}$ are the $L^{2}(M,\rho)$ -adjoint operator of $\mathcal{P},\mathcal{Q}$ . With $\mathcal{P}$ defined as (12), the memory integral of the MZ equation (9) and its projected version (10) can be simplified to a convolution term. Therefore these two MZ equations can be rewritten as

	$\displaystyle\frac{d\bm{u}(t)}{dt}$	$\displaystyle=\bm{\Omega}\bm{u}(t)+\int_{0}^{t}\bm{K}(t-s)\bm{u}(s)ds+\bm{f}(t),$		(13)
	$\displaystyle\frac{d}{dt}\mathcal{P}{\bm{u}}(t)$	$\displaystyle=\bm{\Omega}\mathcal{P}{\bm{u}}(t)+\int_{0}^{t}\bm{K}(t-s)\mathcal{P}{\bm{u}}(s)\,ds,$		(14)

where $\bm{u}(t)=[u_{1}(t),\dots,u_{N}(t)]^{T}$ and


$\displaystyle G_{ij}$	$\displaystyle=\langle u_{i}(0),u_{j}(0)\rangle_{\rho}\quad\text{(Gram matrix)},$	(15a)
$\displaystyle\Omega_{ij}$	$\displaystyle=\sum_{k=1}^{N}G^{-1}_{jk}\langle u_{k}(0),\mathcal{K}u_{i}(0)\rangle_{{\rho}}\quad\text{(streaming matrix)},$	(15b)
$\displaystyle K_{ij}(t)$	$\displaystyle=\sum_{k=1}^{N}G^{-1}_{jk}\langle u_{k}(0),\mathcal{K}e^{t\mathcal{Q}\mathcal{K}}\mathcal{Q}\mathcal{K}u_{i}(0)\rangle_{\rho}\quad\text{(memory kernel)},$	(15c)
$\displaystyle f_{i}(t)$	$\displaystyle=e^{t\mathcal{Q}\mathcal{K}}\mathcal{Q}\mathcal{K}u_{i}(0)\quad\text{(fluctuation term)}.$	(15d)

Equations (13)-(14) are known as the (linear) generalized Langevin equation (GLE) in statistical physics. The projection operator method provides a systematic way to derive such closed equations of motion from the first principle. Depending on the choice of the Hilbert space weight function $\rho$ , the linear GLE (13) and its projected form (14) yield evolution equations for different dynamical quantities. When considering SDE (6) in the context of statistical physics, the most common setting of $\rho$ is $\rho=\rho_{0}=\rho_{S}$ , where $\rho_{S}$ is the steady state distribution of the stochastic system satisfying the Fokker-Planck equation $\partial_{t}\rho_{S}=\mathcal{K}^{*}\rho_{S}=0$ (see details in Section 4). For such a case, GLE (13) yields the full dynamics of the noise-averaged quantity $\mathbb{E}_{\bm{\mathcal{W}}(t)}[\bm{u}(\bm{x}(t))|\bm{x}_{0}]$ and the projected GLE (14) yields the evolution equation of the steady state time-autocorrelation function $C_{ij}(t)$ , which is defined as [36, 50]

\displaystyle C_{ij}(t):=\mathbb{E}_{\bm{x}_{0}}[\mathbb{E}_{\bm{\mathcal{W}}(t)}[u_{i}(t)u_{j}(0)|\bm{x}_{0}]]=\langle\mathcal{M}(t,0)u_{i}(0),u_{j}(0)\rangle_{\rho_{0}}=\langle\mathcal{M}(t,0)u_{i}(0),u_{j}(0)\rangle_{\rho_{S}}.

For deterministic Hamiltonian systems with the unitary evolution operator $e^{t\mathcal{L}}$ , the Mori-Zwanizg equation is as (9) but with the Kolmogorov operator $\mathcal{K}$ replaced by the Liouville operator $\mathcal{L}$ . As we mentioned in the introduction, the second FDT the the thermal equilibirum holds naturally as a result of the skew-symmetry of the Liouville operator $\mathcal{L}$ with respect to $\langle\cdot,\cdot\rangle_{\rho_{eq}}$ , where $\rho_{eq}$ is the equilibrium distribution [43, 46]. In fact we have

$\displaystyle K_{ij}(t)=\sum_{k=1}^{N}G^{-1}_{jk}\langle u_{k}(0),\mathcal{L}e^{t\mathcal{Q}\mathcal{L}}\mathcal{Q}\mathcal{L}u_{i}(0)\rangle_{\rho_{eq}}$	$\displaystyle=-\sum_{k=1}^{N}G^{-1}_{jk}\langle\mathcal{L}u_{k}(0),e^{t\mathcal{Q}\mathcal{L}}\mathcal{Q}\mathcal{L}u_{i}(0)\rangle_{\rho_{eq}}$	(16)
	$\displaystyle=-\sum_{k=1}^{N}G^{-1}_{jk}\langle\mathcal{Q}\mathcal{L}u_{k}(0),e^{t\mathcal{Q}\mathcal{L}}\mathcal{Q}\mathcal{L}u_{i}(0)\rangle_{\rho_{eq}}$
	$\displaystyle=-\sum_{k=1}^{N}G^{-1}_{jk}\langle f_{k}(0),f_{i}(t)\rangle_{\rho_{eq}},$

where we used the idempotence of the projection operator, i.e. $\mathcal{Q}^{2}=\mathcal{Q}$ and the fact that Mori-type projection operator $\mathcal{Q}$ is symmetric in $L^{2}(M,\rho_{eq})$ . For stochastic systems studied in this paper, since the Kolmogorov operator $\mathcal{K}$ is not skew-symmetric with respect to the inner product $\langle\cdot,\cdot\rangle_{\rho}$ , the classical second FDT has to be generalized accordingly.

4 Generalized fluctuation-dissipation theorem

In this section, we first review the derivation of the generalized first FDT for stochastic systems and emphasize its difference from the second FDT. As a mathematical preparation, we need to assume that there exists a unique probability measure $d\mu=\rho d\bm{x}$ which is invariant under the Markovian semigroup generated by SDE (6). This implies the existence and uniqueness of a time-independent probability density $\rho$ that satisfies the Kolmogorov forward (Fokker-Planck) equation:

\displaystyle\partial_{t}\rho=\mathcal{K}^{*}\rho=0,

(17)

where $\mathcal{K}^{*}$ is the $L^{2}(M)$ -adjoint operator of $\mathcal{K}$ . Throughout the paper, we further assume that $\rho$ is a smooth function which decays to $0$ when approaching boundary $\partial M$ . For frequently used statistical physics models, these two assumptions are generally hard to prove since it involves the analysis of the degenerate elliptic operator $\mathcal{K}$ . Some recent studies on hypoellipticity [6, 8] have shown the possibility to obtain affirmative answers using the rather complicated Hörmander analysis. These theoretical results, together with numerical studies such as [25, 26], suggest that the uniqueness, smoothness and the decaying properties of $\rho$ at $\partial M$ are rather technical assumptions.

4.1 Derivation of the first FDT

The first fluctuation–dissipation theorem describes the linear response of a given dynamical system to external perturbations. This relationship was established by Kubo [20] for classical and quantum Hamiltonian systems and then has been extended to open systems by different researchers [1, 14, 3, 12]. As we briefly mentioned in the introduction, the derivation of the first FDT relies on the perturbation theory. To demonstrate this point, we consider a general perturbation to the stochastic system (6):

\displaystyle d\bm{x}(t)=\bm{F}(\bm{x}(t))dt+\bm{\sigma}(\bm{x}(t))d\bm{\mathcal{W}}(t)+\delta\bm{G}(t)\cdot\tilde{\bm{F}}(\bm{x}(t))dt+\sqrt{\delta}\bm{H}(t)\cdot\tilde{\bm{\sigma}}(\bm{x}(t))d\bm{\mathcal{W}}(t),\qquad\bm{x}(0)\sim\rho,

(18)

where $0<\delta\ll 1$ , $\delta\bm{G}(t)\cdot\tilde{\bm{F}}(\bm{x}(t))dt$ and $\sqrt{\delta}\bm{H}(t)\cdot\tilde{\bm{\sigma}}(\bm{x}(t))d\bm{\mathcal{W}}(t)$ characterize respectively the dynamics of the pertubative frictional and fluctuating forces added to the stochastic system. We denote $\rho$ as the steady state distribution for the unperturbed stochastic system (6) and $\rho_{\delta}$ as the one for the perturbed stochastic system (18). For the sake of simplicity, we assume that $\delta\bm{G}(t)$ , $\sqrt{\delta}\bm{H}(t)$ , the diffusion matrix $\bm{\sigma}(\bm{x}(t))$ and the perturbed diffusion matrix $\tilde{\bm{\sigma}}(\bm{x}(t))$ are all $d\times d$ diagonal matrices. As a consequence, the diffusion part (second-order derivatives) of the Kolmogorov backward operator corresponding to the perturbed system is of the diagonal form and we have

\displaystyle\mathcal{K}=\sum_{i=1}^{d}F_{i}(\bm{x})\partial_{x_{i}}+\sigma_{i}^{2}(\bm{x})\partial_{x_{i}}^{2}+\delta G_{i}(t)\tilde{F_{i}}(\bm{x})\partial_{x_{i}}+\delta H^{2}_{i}(t)\tilde{\sigma}_{i}^{2}(\bm{x})\partial_{x_{i}}^{2},

where $\sigma_{i}(\bm{x}),\tilde{\sigma}_{i}(\bm{x})$ , $\delta G_{i}(t)$ and $\delta H_{i}^{2}(t)$ are the $i$ -th diagonal elements of the corresponding matrices. Using the standard perturbation theory, we obtain the following generalized first FDT:

Theorem 1.

(Generalized-1st-FDT) Assuming that the steady state distribution $\rho=\rho(\bm{x})$ is a smooth function of $\bm{x}$ and decays to 0 when approaching the boundary $\partial M$ . For the perturbed SDE (18), the following generalized first FDT holds for state space observable $u=u(\bm{x})$ :

	$\displaystyle\langle u(t)\rangle_{\rho_{\delta}}-\langle u(t)\rangle_{\rho}=-\sum_{i=1}^{d}\int_{0}^{t}$	$\displaystyle\left\langle\frac{1}{\rho}[\partial_{x_{i}}\tilde{F}_{i}\rho]u(t-s)\right\rangle_{\rho}\delta G_{i}(s)ds$		(19)
		$\displaystyle+\sum_{i=1}^{d}\int_{0}^{t}\left\langle\left(\partial_{x_{i}}^{2}\tilde{\sigma}_{i}+\frac{2}{\rho}\partial_{x_{i}}\tilde{\sigma}_{i}\partial_{x_{i}}\rho+\frac{1}{\rho}\tilde{\sigma}_{i}\partial_{x_{i}}^{2}\rho\right)u(t-s)\right\rangle_{\rho}\delta H_{i}^{2}(s)ds+O(\delta^{2}).$		(19)

If the perturbative forces $\delta\bm{G}(t)$ and $\sqrt{\delta}\bm{H}(t)$ are homogeneous, i.e. $\delta G_{i}(t)=\delta G(t)$ and $\sqrt{\delta}H_{i}(t)=\sqrt{\delta}H(t)$ , then the generalized first FDT (19) has a vector-form representation:

	$\displaystyle\langle u(t)\rangle_{\rho_{\delta}}-\langle u(t)\rangle_{\rho}=-\int_{0}^{t}$	$\displaystyle\left\langle\frac{1}{\rho}\nabla\cdot[\tilde{\bm{F}}\rho]u(t-s)\right\rangle_{\rho}\delta G(s)ds$		(20)
	$\displaystyle+$	$\displaystyle\int_{0}^{t}\left\langle\left(\nabla\cdot\overrightarrow{L_{\tilde{\bm{\sigma}}}}+\frac{2}{\rho}\overrightarrow{L_{\tilde{\bm{\sigma}}}}\cdot\nabla\rho+\frac{1}{\rho}\mathcal{D}_{\tilde{\bm{\sigma}}}\rho\right)u(t-s)\right\rangle_{\rho}\delta H^{2}(s)ds+O(\delta^{2}),$		(20)

where the vector function $\overrightarrow{L_{\tilde{\bm{\sigma}}}}=[\partial_{x_{i}}\tilde{\sigma_{i}}]_{i=1}^{d}$ and $\mathcal{D}_{\tilde{\bm{\sigma}}}$ is a second-order differential operator defined as $\mathcal{D}_{\tilde{\bm{\sigma}}}=\sum_{i=1}^{d}\tilde{\sigma_{i}}\partial_{x_{i}}^{2}$ .

The proof of Theorem 1 follows a perturbation analysis similar to the one used in [20] or more recent papers such as [12, 3]. We also outline the procedure in A. In the absence of the stochastic perturbation, i.e. $\sqrt{\delta}\bm{H}(t)=0$ , for homogeneous perturbation $\delta G_{i}(t)=\delta G(t)$ , we get the commonly used linear-response relation [12, 3]:

\displaystyle\langle u(t)\rangle_{\rho_{\delta}}=-\int_{0}^{t}\left\langle\frac{1}{\rho}\nabla\cdot[\tilde{\bm{F}}\rho]u(t-s)\right\rangle_{\rho}\delta G(s)ds.

(21)

We note that to apply formulas (19)-(21), it is required to know the steady state distribution $\rho$ . For nonequilibrium systems, this is generally hard to obtain due to high dimensionality of the Fokker-Planck equation (17). Various approaches such as the Gaussian approximation [12] and the information theory method [30] were proposed to address this issue. As we will see in the following section, the generalized second FDT also contains additional terms that involve $\rho$ . However, for specifically chosen observables $u(\bm{x})$ , calculation of the density $\rho$ can be avoided.

4.2 Derivation of the second FDT

The second flutucation-dissipation theorem gives the proportionality between the noise amplitude and the friction kernel in the GLE. For open systems, depending on the form of the GLE and the observable of interest, the explicit expression of the second FDT will be different. A well-organized review in this regard is given by Maes in [28]. In this section, we provide a novel way to generalize the second FDT using the first-principle GLE derived in Section 3. The derivation only uses the mathematical properties of the Kolmogorov operator $\mathcal{K}$ . Without loss of generality, we consider the SDE (6) with the Kolmogorov (backward) operator (8) of the form:

\displaystyle\mathcal{K}=\mathcal{L}+\mathcal{D}=\sum_{i=1}^{N}F_{i}(\bm{x})\partial_{x_{i}}+\sum_{i,j=1}^{N}\sigma_{i}\sigma_{j}\partial^{2}_{x_{i}x_{j}},

(22)

where $\sigma_{i},\sigma_{j}$ are constants. Note that here $\mathcal{L}$ represents a general advection operator instead of the skew-symmetric Liouville operator. For such a stochastic system in the steady state $\rho$ , by introducing the Mori-type projection operator (12) in weighted Hilbert space $L^{2}(M,\rho)$ , we can prove the following generalized second FDT for the GLE (13)-(14). This is the main theoretical result of this paper.

Theorem 2.

(Generalized-2nd-FDT) For SDE with infinitesimal generator (22), if we assume that the steady state distribution $\rho$ is a smooth function of $\bm{x}$ and decays to $0$ when approaching the boundary $\partial M$ , then the following generalized second FDT holds for any state space observable $u=u(\bm{x})$ in GLE (13)-(14):

\displaystyle K(t)=-\frac{\langle f(0),f(t)\rangle_{\rho}}{\langle u^{2}(0)\rangle_{\rho}}+\frac{\langle w(0),f(t)\rangle_{\rho}}{\langle u^{2}(0)\rangle_{\rho}},

(23)

where $w(0)=\mathcal{S}u(0)$ and $\mathcal{S}$ is a second-order differential operator defined as

\displaystyle\mathcal{S}=2\sum_{i,j=1}^{N}\sigma_{i}\sigma_{j}\partial^{2}_{x_{i}x_{j}}+\frac{1}{\rho}\sum_{i,j=1}^{N}\sigma_{i}\sigma_{j}\partial_{x_{i}}\rho\partial_{x_{j}}+\sigma_{i}\sigma_{j}\partial_{x_{j}}\rho\partial_{x_{i}}.

(24)

For a special case where $\mathcal{K}$ is a degenerate elliptic operator and the diffusion term is of the diagonal form $\mathcal{D}=\sigma^{2}\Delta$ , where $\Delta$ is a $n$ -dimensional Laplacian with $n\leq N$ , then $\mathcal{S}$ is an $n$ -dimensional operator and admits a simple form:

\displaystyle\mathcal{S}=2\sigma^{2}\Delta+2\sigma^{2}\nabla(\ln\rho)\cdot\nabla.

(25)

Proof.

For any operator $\mathcal{O}$ defined in weighted Hilbert space $L^{2}(M,\rho)$ , where $\rho$ is the steady state distribution, we denote $\mathcal{O}^{*}_{\rho}$ as the adjoint operator of $\mathcal{O}$ in $L^{2}(M,\rho)$ with respect to inner product $\langle\cdot,\cdot\rangle_{\rho}$ . Correspondingly, we denote $\mathcal{O}^{*}$ as the adjoint operator of the flat Hilbert space $\mathcal{O}$ in $L^{2}(M)$ with respect to the weight-less inner product $\langle\cdot,\cdot\rangle$ . First we aim to find the $L^{2}(M,\rho)$ adjoint of the Kolmogorov operator $\mathcal{K}$ (22). Using the integration by parts formula and the fact that $\rho(\bm{x})=0$ at $\partial M$ , it is easy to obtain

\displaystyle\mathcal{L}+\mathcal{L}^{*}_{\rho}

\displaystyle=-\nabla\cdot\bm{F}(\bm{x})-\frac{1}{\rho}\mathcal{L}\rho.

(26)

Similarly, for the second-order differential operator $\mathcal{D}$ in (22), using the integration by parts formula twice leads to

\displaystyle\mathcal{D}^{*}_{\rho}=\mathcal{D}+\frac{1}{\rho}\sum_{i,j=1}^{N}\sigma_{i}\sigma_{j}\partial_{x_{i}}\rho\partial_{x_{j}}+\sigma_{i}\sigma_{j}\partial_{x_{j}}\rho\partial_{x_{i}}+\frac{1}{\rho}\mathcal{D}\rho.

(27)

Naturally we have

\displaystyle\mathcal{D}+\mathcal{D}^{*}_{\rho}=2\mathcal{D}+\frac{1}{\rho}\sum_{i,j=1}^{N}\sigma_{i}\sigma_{j}\partial_{x_{i}}\rho\partial_{x_{j}}+\sigma_{i}\sigma_{j}\partial_{x_{j}}\rho\partial_{x_{i}}+\frac{1}{\rho}\mathcal{D}\rho=\mathcal{S}+\frac{1}{\rho}\mathcal{D}\rho,

(28)

where operator $\mathcal{S}$ is defined as (24). For the special case where $\mathcal{D}=\sigma^{2}\Delta$ , it is easy to see that $\mathcal{S}$ has the simple form (25). Summing up formulas (26) and (28), we can get that

\displaystyle\mathcal{K}_{\rho}^{*}+\mathcal{K}=-\nabla\cdot\bm{F}(\bm{x})-\frac{1}{\rho}\mathcal{L}\rho+\frac{1}{\rho}\mathcal{D}\rho+\mathcal{S}.

(29)

Using a similar procedure for operator $\mathcal{L}$ and $\mathcal{D}$ defined in the flat Hilbert space $L^{2}(M)$ , we have

\displaystyle\mathcal{L}+\mathcal{L}^{*}=-\nabla\cdot\bm{F}(\bm{x}),\qquad\mathcal{D}^{*}=\mathcal{D}.

(30)

On the other hand, since $\rho$ satisfies the steady state Fokker-Planck equation $\partial_{t}\rho=\mathcal{K}^{*}\rho=0$ , we obtain an operator identity $\partial_{t}\rho=\mathcal{K}^{*}\rho=\mathcal{L}^{*}\rho+\mathcal{D}^{*}\rho=0$ . Combining this with formula (30) and noting that $-\nabla\cdot\bm{F}(\bm{x})$ is a multiplication operator, we obtain

\left.\begin{aligned} \mathcal{L}\rho+\mathcal{L}^{*}\rho&=-\rho\nabla\cdot\bm{F}(\bm{x})\\ \mathcal{D}^{*}\rho&=\mathcal{D}\rho\end{aligned}\right\}\Rightarrow-\frac{1}{\rho}\mathcal{L}\rho+\frac{1}{\rho}\mathcal{D}\rho-\nabla\cdot\bm{F}(\bm{x})=0.

(31)

Substituting this relation into (29) we get $\mathcal{K}^{*}_{\rho}=-\mathcal{K}+\mathcal{S}$ . Now we note that symmetric projection operator $\mathcal{Q}$ satisfies $\mathcal{Q}^{*}=\mathcal{Q}=\mathcal{Q}^{2}$ . For the formally defined GLE memory kernel (15c), we obtain the generalized second FDT:

	$\displaystyle\langle u^{2}(0)\rangle_{\rho}K(t)=\langle u(0),\mathcal{K}e^{t\mathcal{Q}\mathcal{K}}u(0)\rangle_{\rho}$	$\displaystyle=\langle\mathcal{K}^{*}_{\rho}u(0),e^{t\mathcal{Q}\mathcal{K}}\mathcal{Q}\mathcal{K}u(0)\rangle_{\rho}$
		$\displaystyle=-\langle\mathcal{Q}\mathcal{K}u(0),e^{t\mathcal{Q}\mathcal{K}}\mathcal{Q}\mathcal{K}u(0)\rangle_{\rho}+\langle w(0),e^{t\mathcal{Q}\mathcal{K}}\mathcal{Q}\mathcal{K}u(0)\rangle_{\rho}$
		$\displaystyle=-\langle f(0),f(t)\rangle_{\rho}+\langle w(0),f(t)\rangle_{\rho}.\$

∎

Theorem 2 can be readily generalized to the $N$ -dimensional GLE (13)-(14) and for the general Kolmogorov operator (8). When the dissipative term in the Kolmogorov operator (22) is given by $\mathcal{D}=\sum_{i,j=1}^{N}\sigma_{i}(\bm{x})\sigma_{j}(\bm{x})\partial^{2}_{x_{i},x_{j}}$ , then we obtain $w(0)=\mathcal{S}u(0)$ with

\displaystyle\mathcal{S}=\sum_{i,j=1}^{N}2\sigma_{i}\sigma_{j}\partial^{2}_{x_{i}x_{j}}+\frac{1}{\rho}\sigma_{i}\sigma_{j}\partial_{x_{i}}\rho\partial_{x_{j}}+\frac{1}{\rho}\sigma_{i}\sigma_{j}\partial_{x_{j}}\rho\partial_{x_{i}}+\partial_{x_{i}}[\sigma_{i}\sigma_{j}]\partial_{x_{j}}+\partial_{x_{j}}[\sigma_{i}\sigma_{j}]\partial_{x_{i}},

(32)

where the shorthand notation $\sigma_{i}(\bm{x})=\sigma_{i}$ and $\sigma_{j}(\bm{x})=\sigma_{j}$ are used. The derivation of (32) follows directly from the proof of Theorem 2 and is given in B. On the other hand, the result of Theorem 2 holds for an arbitrary Mori-type projection operator $\mathcal{P}$ . Hence for the $N$ -dimensional GLEs (13)-(14), we have

\displaystyle K_{ij}(t)=\sum_{k=1}^{N}G^{-1}_{jk}(-\langle f_{k}(0),f_{i}(t)\rangle_{\rho}+\langle\mathcal{S}u_{k}(0),f_{i}(t)\rangle_{\rho}).

(33)

When one tries to apply the first and the second FDTs such as (21) and (23), it is normally necessary to know the exact form of the steady state distribution $\rho$ , which, however, is generally hard to deduce for nonequilibrium systems such as turbulence [30] and heat conduction models [26]. This difficulty can be bypassed because of the following corollary which is already announced in Section 2:

Corollary 2.1.

If the phase space observable $u=u(\bm{p}(t),\bm{q}(t))$ is only a function of the degenerate coordinates of the stochastic system, then $w(0)=\mathcal{S}u(0)=0$ and the GLE (13)-(14) for $u(t)$ satisfies the classical second FDT.

The proof is obvious since $\mathcal{S}$ is an operator in the non-degenerate coordinate. A typical example is the Langevin dynamics for a molecular dynamical system. Since the Gaussian white noise is only imposed in the momentum coordinate $p_{i}$ , the position coordinate $q_{i}$ is therefore degenerate. Hence we have $w(0)=\mathcal{S}u(0)=\mathcal{S}f(q_{i}(0))=0$ . For more general, nonequilibrium systems in the steady state, this result still holds which is somewhat surprising because the nonequilibrium steady state (NESS) measure is generally unknown! As we will see in the following section, many statistical physics models are generated by highly degenerate elliptic operator $\mathcal{K}$ , i.e. the dissipative forces act on a small subset of the phase space coordinates. For such models, if the observable $u(\bm{x})$ of interest is a function of the degenerate coordinate, we can simply avoid the calculation of $\rho$ in evaluating $w(0)$ and use the classical second FDT to build reduced-order models for the observable. An application of this fact in heat conduction problems is presented in Section 6. We also have the following result:

Corollary 2.2.

For SDE with infinitesimal generator (8), if the linear GLEs (13)-(14) for observable $\bm{u}(t)$ satisfying the classical second FDT with $\bm{w}(0)=0$ , then $\Omega_{ij}=0$ in (13)-(14).

Proof.

It is sufficient to prove the one-dimensional case. According to the definition of $\Omega_{ij}$ (15b), for one-dimensional GLE (13)-(14), we have

\displaystyle\Omega=\frac{\langle\mathcal{K}u(0),u(0)\rangle_{\rho}}{\langle u^{2}(0)\rangle_{\rho}}=\frac{\langle u(0),\mathcal{K}_{\rho}^{*}u(0)\rangle_{\rho}}{\langle u^{2}(0)\rangle_{\rho}}=-\frac{\langle u(0),\mathcal{K}u(0)\rangle_{\rho}}{\langle u^{2}(0)\rangle_{\rho}}+\frac{\langle u(0),w(0)\rangle_{\rho}}{\langle u^{2}(0)\rangle_{\rho}}=-\frac{\langle u(0),\mathcal{K}u(0)\rangle_{\rho}}{\langle u^{2}(0)\rangle_{\rho}}=0.

∎

Remark

When compared with previous methods, e.g. [28], in deriving the (generalized) second FDT for stochastic systems, our approach does not rely on weak perturbation assumptions therefore is generally applicable to arbitrary stochastic systems driven by white noise. Moreover, the derivation only uses the properties of the Kolmogorov operator without applying any ad hoc approximations. On the other hand, although the presented example is based on the Mori-type projection operator which leads to linear GLEs, the above theory also applies to nonlinear GLEs such as the Zwanzig’s equation [52], essentially because the Zwanzig-type projection operator is an infinite-rank operator which is similar to the Mori-type projection operator. The proof is rather technical therefore will be deferred to C.

5 The generalized second FDT for specific systems

5.1 The generalized second FDT for equilibrium systems

The application of the generalized FDTs to equilibrium systems leads to explicit expressions of formula (19), (20) and (23) since the equilibrium distribution $\rho$ is given by the Gibbs-Boltzmann form $\rho=e^{-\beta\mathcal{H}}/Z$ for the canonical ensemble. In this section, we will derive such expression for some frequently used statistical mechanics models. Before we move onto analyzing stochastic systems, it is worth noticing that for equilibrium systems generated by deterministic forces such as the Nosé-Hoover thermostats, the classical second FDT holds as a result of the skew-symmetry of the Liouville operator $\mathcal{L}$ .

Langevin dynamics

The Langevin dynamics for a $d$ -dimensional system of $N$ interacting particles is given by the following SDE in $\mathbb{R}^{2d\times N}$ :

\displaystyle\begin{dcases}d\bm{q}_{i}=\frac{1}{m_{i}}\bm{p}_{i}dt\\ d\bm{p}_{i}=\sum_{i\neq j}^{N}\bm{F}_{i,j}^{C}dt-\frac{\gamma}{m_{i}}\bm{p}_{i}dt+\sigma d\bm{\mathcal{W}}_{i}(t)\end{dcases},

(34)

where $m_{i}$ is the mass of each particle, $\sum_{i\neq j}^{N}\bm{F}_{i,j}^{C}$ is the total conservative force acting on particle $i$ , and $\bm{\mathcal{W}}_{i}(t)$ is a $d$ -dimensional Wiener process which satisfies $d\bm{\mathcal{W}}_{j}(t)=\bm{\xi}_{i}(t)dt$ with

\displaystyle\langle\xi_{ij}(t)\rangle=0,\qquad\langle\xi_{ij}(t)\xi_{i^{\prime}j^{\prime}}(s)\rangle=2\gamma k_{B}T\delta_{ii^{\prime}}\delta_{jj^{\prime}}\delta(t-s).

The parameter $\gamma$ is the friction coefficient which is related to $\sigma$ through the fluctuation-dissipation relation $\sigma=(2\gamma/\beta)^{1/2}$ , where $\beta=1/k_{B}T$ , $k_{B}$ is the Boltzmann constant and $T$ the temperature of the equilibrium system. The stochastic dynamical system (34) is widely used in the mesoscopic modelling of liquids and gases. The Kolmogorov backward operator (8) associated with the SDE (34) is given by

\displaystyle\mathcal{K}=\sum_{i=1}^{N}\frac{\bm{p}_{i}}{m_{i}}\cdot\partial_{\bm{q}_{i}}+\sum_{i,j\neq i}^{N}\bm{F}_{i,j}^{C}\cdot\partial_{\bm{p}_{i}}-\sum_{i}^{N}\frac{\gamma\bm{p}_{i}}{m_{i}}\cdot\partial_{\bm{p}_{i}}+\sum_{i}^{N}\frac{\gamma}{\beta}\partial_{\bm{p}_{i}}\cdot\partial_{\bm{p}_{i}},

(35)

where “ $\cdot$ ” denotes the standard dot product. If the interaction potential $V(\bm{q})$ is strictly positive at infinity then the Langevin equation (34) admits a unique invariant Gibbs measure given by

\rho_{eq}(\bm{p},\bm{q})=\frac{1}{Z}e^{-\beta\mathcal{H}(\bm{p},\bm{q})},

(36)

where

\mathcal{H}(\bm{p},\bm{q})=\sum_{i=1}^{N}\frac{\|\bm{p}_{i}\|_{2}^{2}}{2m_{i}}+V(\bm{q})

(37)

is the Hamiltonian and $Z$ is the partition function. The formal expression of the Gibbs-Boltzmann distribution enables us to get the explicit expression of the additional term $w(0)$ in the generalized second FDT (23). As an example, the Langevin dynamics (34) is often used to study the self-diffusion of Brownian particles, for which a relevant physical observable is the tagged particle velocity $\bm{v}_{j}$ . By choosing $u(0)=p_{jx}/m_{j}$ and using (23) and (35), we obtain

\displaystyle w(0)=\sum_{i=1}^{N}\frac{\gamma}{\beta}\partial_{\bm{p}_{i}}\cdot\partial_{\bm{p}_{j}}\frac{p_{jx}}{m_{j}}+\sum_{i=1}^{N}\frac{\sigma^{2}}{\rho_{eq}}\partial_{\bm{p}_{i}}\rho_{eq}\cdot\partial_{\bm{p}_{i}}\frac{p_{jx}}{m_{j}}=-\frac{2\gamma}{m_{j}^{2}}p_{jx}.

(38)

This implies that the GLE (13) for observable $p_{jx}$ can be rewritten as

\displaystyle\frac{d}{dt}p_{jx}(t)=\Omega p_{jx}(t)-\int_{0}^{t}\left\langle f(0)+\frac{2\gamma}{m_{j}^{2}}p_{jx}(0),f(t-s)\right\rangle_{\rho_{eq}}p_{jx}(s)ds+f(t).

Remark

Here we obtained a somewhat counterintuitive conclusion stating that the classical second FDT for a Brownian particle does not hold in the statistical equilibrium since we have an additional term (38) in the GLE memory kernel. This is because instead of using Hamiltonian dynamics, we used stochastic dynamics, i.e. Eqn (34), to simulate the equilibrium. On the other hand, the form of the second FDT strongly depends on the GLE under investigation. For other GLEs such as the nonlinear ones derived by Zwanzig-type projection operators [24, 27, 17], it is possible that the classical second FDT still holds for a Brownian particle generated by the Langevin dynamics.

Dissipative particle dynamics

For a $d$ -dimensional interacting particle system of $N$ particles, the SDE that governs the particle position $\bm{q}_{i}$ and momentum $\bm{p}_{i}$ in dissipative particle dynamics (DPD) is given by [16, 11]

\displaystyle\begin{dcases}d\bm{q}_{i}&=\frac{\bm{p}_{i}}{m_{i}}dt\\ d\bm{p}_{i}&=\sum_{i\neq j}^{N}\bm{F}_{ij}^{C}(\bm{q}_{ij})dt-\sum_{i\neq j}^{N}\gamma\omega(q_{ij})(\bm{e}_{ij}\cdot\bm{v}_{ij})\bm{e}_{ij}dt+\sum_{i\neq j}^{N}\sigma\omega^{1/2}(q_{ij})\bm{e}_{ij}d\mathcal{W}_{ij}(t)\end{dcases}

(39)

where $m_{i}$ is the mass of $i$ -th particle, $\bm{q}_{ij}=\bm{q}_{i}-\bm{q}_{j}$ , $q_{ij}=\|\bm{q}_{i}-\bm{q}_{j}\|$ , $\bm{e}_{ij}=\bm{q}_{ij}/q_{ij}$ , $\bm{v}_{ij}=\bm{v}_{i}-\bm{v}_{j}$ , $\bm{v}_{i}=\bm{p}_{i}/m_{i}$ and $\bm{F}^{C}_{ij}(\bm{q}_{ij})$ is the conservative force exerted on particle $i$ by particle $j$ . The dimensionless weight function $\omega(q_{ij})$ provides the range of interactions of the dissipative and random forces. The friction coefficient and the noise intensity are linked with each other through the fluctuation-dissipation relation $\sigma=(2\gamma/\beta)^{1/2}$ , where $\beta=1/k_{B}T$ . For the DPD model, the frictional forces are applied in a pair-wise form, such that the sum of thermostating forces acting on a particle pair equals zero. Hence for $d\mathcal{W}_{ij}(t)=\bm{\xi}_{i}(t)dt$ , we have $\xi_{ij}(t)=\xi_{ji}(t)$ and

\displaystyle\langle\xi_{ij}(t)\xi_{i^{\prime}j^{\prime}}(t)\rangle=(\delta_{ii^{\prime}}\delta_{jj^{\prime}}+\delta_{ij^{\prime}}\delta_{ji^{\prime}})\delta(t-s).

The Kolmogorov backward operator associated with the DPD model (39) is given by [10]

\displaystyle\mathcal{K}=\sum_{i=1}^{N}\frac{\bm{p}_{i}}{m_{i}}\cdot\partial_{\bm{q}_{i}}+\sum_{i,j\neq i}^{N}\bm{F}_{i,j}^{C}\cdot\partial_{\bm{p}_{i}}-\sum_{i,j\neq i}^{N}\gamma\omega(q_{ij})(\bm{e}_{ij}\cdot\bm{v}_{ij})\bm{e}_{ij}\cdot\partial_{\bm{p}_{i}}+\frac{1}{\beta}\sum_{i,j\neq i}^{N}\gamma\omega(q_{ij})(\partial_{\bm{p}_{i}}\cdot\partial_{\bm{p}_{i}}-\partial_{\bm{p}_{i}}\cdot\partial_{\bm{p}_{j}}).

Similarly, for the $x$ -directional velocity of the $h$ -th particle in the DPD model, the additional term $w(0)$ in the generalized second FDT can be calculated using (32) as:

	$\displaystyle w(0)$	$\displaystyle=\mathcal{S}\frac{p_{hx}}{m_{h}}$
		$\displaystyle=\frac{1}{\beta}\sum_{i,j}^{N}\frac{\sigma_{i}\sigma_{j}\delta_{hj}\delta_{ij}}{m_{h}\rho_{eq}}\partial_{p_{ix}}\rho_{eq}+\frac{\sigma_{i}\sigma_{j}\delta_{hi}\delta_{ij}}{m_{h}\rho_{eq}}\partial_{p_{jx}}\rho_{eq}-\frac{1}{\beta}\sum_{i,i\neq j}^{N}\frac{\sigma_{i}\sigma_{j}\delta_{hj}}{m_{h}\rho_{eq}}\partial_{p_{ix}}\rho_{eq}+\frac{\sigma_{i}\sigma_{j}\delta_{hi}}{m_{h}\rho_{eq}}\partial_{p_{jx}}\rho_{eq}$
		$\displaystyle=-\sum_{i,j}^{N}\frac{\sigma_{i}\sigma_{j}\delta_{hj}\delta_{ij}}{m_{h}m_{i}}p_{ix}+\frac{\sigma_{i}\sigma_{j}\delta_{hj}\delta_{ij}}{m_{h}m_{j}}p_{jx}+\sum_{i,i\neq j}^{N}\frac{\sigma_{i}\sigma_{j}\delta_{hj}}{m_{h}m_{i}}p_{ix}+\frac{\sigma_{i}\sigma_{j}\delta_{hi}}{m_{h}m_{j}}p_{jx},$

where $\sigma_{i}\sigma_{j}=\gamma\omega(q_{ij})$ and we used the fact that $\partial_{\bm{p}_{i}}[\sigma_{i}\sigma_{j}]=0$ . The resulting GLEs (13)-(14) can be obtained accordingly. For the Langevin dynamics and the DPD model, the corresponding Kolmogorov operators are both degenerate elliptic operators since white noise is only imposed in the momentum space. Hence if we choose the tagged particle position $\bm{q}_{j}$ as the quantity of interest, then the classical second FDT holds for GLEs (13)-(14) as claimed in Section 4.

5.2 The generalized second FDT for nonequilibrium systems

Refer to caption — Figure 1: Schematic of the $1$ -D and $2$ -D heat transport model (41). The onsite potential and boundary Langevin forces are shown in the top plot for the $1$ -D model. Note that only the boundary oscillators, marked in red, interact with the heat bath.

As an example of nonequilibrium systems, we consider an $d$ -dimensional heat conduction model in [6, 26]. To this end, we consider a lattice $\mathcal{G}$ of interacting oscillators. For each oscillator $i\in\mathcal{G}$ , the position and momentum are given respectively by $\bm{q}_{i}\in\mathbb{R}^{d}$ and $\bm{p}_{i}\in\mathbb{R}^{d}$ . The phase space is therefore given by $\Omega=\mathbb{R}^{2|\mathcal{G}|d}$ , where $|\mathcal{G}|$ gives the cardinality of the set $\mathcal{G}$ which corresponds to the total number of particles. These oscillators are interacting with the substrate and each other through the Hamiltonian

\displaystyle\mathcal{H}(\bm{p},\bm{q})=\sum_{i\in\mathcal{G}}\left(\frac{\bm{p}_{i}^{2}}{2m_{i}}+U_{i}(\bm{q}_{i})\right)+\sum_{e\in\mathcal{E}}V_{e}(\delta\bm{q}_{e}),

(40)

where $U_{i}$ and $V_{e}$ are the pining potential and interactive potential, respectivly, For $e=(i,i^{\prime})\in\mathcal{E}=\mathcal{G}\times\mathcal{G}$ we have $\delta\bm{q}_{e}=\bm{q}_{i^{\prime}}-\bm{q}_{i}\in\mathbb{R}^{d}$ . The model $(\mathcal{G},\mathcal{E})$ can be viewed as an undirected graph with no on-site loop, i.e. self-interactions of the kind $V_{0}(\delta\bm{q}_{e}=0)$ are not allowed. The graph is undirected in the sense that the interactive potential $V_{e}(\bm{q}_{i^{\prime}}-\bm{q}_{i})=V_{e}(\bm{q}_{i}-\bm{q}_{i^{\prime}})$ appears only once in the Hamiltonian (40). Without the loss of generality, we assume the uniform mass condition $m_{i}=m=1$ . We now choose a subset of the boundary oscillators $\mathcal{B}\subset\partial\mathcal{G}$ to impose thermal baths. For each $b\in\mathcal{B}$ we assume that a thermostat of the temperature $T_{b}>0$ is given, along with a coupling constant $\gamma_{b}>0$ . With this setting, we obtain an $n$ -dimensional heat conduction model given by the following system of stochastic differential equations:

\begin{dcases}\begin{aligned} d\bm{q}_{i}&=\bm{p}_{i}dt\\ d\bm{p}_{i}&=-\partial_{\bm{q}_{i}}\mathcal{H}(\bm{p},\bm{q})dt-\gamma_{i}\bm{p}_{i}dt+\sqrt{2k_{B}T_{i}\gamma_{i}}d\bm{\mathcal{W}}_{i}(t)\end{aligned}\end{dcases}\qquad\qquad i\in\mathcal{G}

(41)

where $\gamma_{i}=0$ for $i\in\mathcal{G}\setminus\mathcal{B}$ and $\gamma_{i}=\gamma_{b}$ , $T_{i}=T_{b}$ for $b\in\mathcal{B}$ . Since the Langevin forces only act on the boundary oscillators, for a “bulk” of oscillators that are away from the boundaries, the dynamics are kept deterministic and Hamiltonian. For modeling purposes, different boundary conditions can be specified. Typical choices are periodic, fixed or free boundaries. No matter which condition is used, the form of the Kolmogorov backward operator $\mathcal{K}$ corresponding to SDE (41) can be written as:

\displaystyle\mathcal{K}=\sum_{b\in\mathcal{B}}(\gamma_{b}k_{B}T_{b}\partial_{\bm{p}_{b}}\cdot\partial_{\bm{p}_{b}}-\gamma_{b}\bm{p}_{b}\cdot\partial_{\bm{p}_{b}})+\sum_{i\in\mathcal{G}}(\bm{p}_{i}\cdot\partial_{\bm{q}_{i}}-\partial_{\bm{q}_{i}}\mathcal{H}(\bm{p},\bm{q})\cdot\partial_{\bm{p}_{i}}).

(42)

When all added thermostats have the same temperature $T_{b}=1/\beta,b\in\mathcal{B}$ , then the system admits a unique invariant measure

\displaystyle d\mu_{\beta}=\rho_{eq}d\bm{p}d\bm{q}=\frac{1}{Z}e^{-\beta\mathcal{H}(\bm{p},\bm{q})}d\bm{p}d\bm{q},

(43)

which is known as the Gibbs measure for thermal equilibrium. When the boundary temperatures are different, the Gibbs measure is no longer invariant and (41) describes the dynamics of heat flowing from the higher-temperature thermostats to the lower-temperature thermostats. Under certain assumptions to the potential energy (see [6]), it can be proved that the system approaches to a unique, steady state exponentially fast. In the literature, such a state is often called the nonequilibrium steady state (NESS). We further denote the steady state probability density as $\rho_{\text{NESS}}$ . By introducing the Mori-type projection operator (12) in the Hilbert space $L^{2}(\mathbb{R}^{2|\mathcal{G}|d},\rho_{\text{NESS}})$ , we can derive the GLE (13)-(14). According to Theorem 2, we can obtain the following proposition:

Proposition 1.

For an n-dimensional heat conduction model given by the SDE (41), if the potential energy $U_{i}(\bm{q}_{i})$ and $V_{e}(\delta\bm{q}_{e})$ satisfy certain conditions which ensure the smoothness and uniqueness of $\rho_{NESS}$ , e.g. the one outlined in [6], then the following generalized second FDT holds for state space observable $u=u(\bm{p},\bm{q})$ in GLEs (13)-(14):

\displaystyle K(t)=-\frac{\langle f(0),f(t)\rangle_{\rho_{\text{NESS}}}}{\langle u^{2}(0)\rangle_{\rho_{\text{NESS}}}}+\frac{\langle w(0),f(t)\rangle_{\rho_{\text{NESS}}}}{\langle u^{2}(0)\rangle_{\rho_{\text{NESS}}}},

where the additional term $w(0)$ is given by

\displaystyle w(0)=2\sum_{b\in\mathcal{B}}\gamma_{b}k_{B}T_{b}\partial_{\bm{p}_{b}}\cdot\partial_{\bm{p}_{b}}u(0)+\frac{2}{\rho_{\text{NESS}}}\sum_{b\in\mathcal{B}}\gamma_{b}k_{B}T_{b}\partial_{\bm{p}_{b}}\rho_{NESS}\cdot\partial_{\bm{p}_{b}}u(0).

In particular, if observable $u$ is a function of bulk coordinates $\mathcal{G}\setminus\mathcal{B}$ , i.e. $u=u(\bm{p},\bm{q})=u([\bm{p}_{i}]_{i=1}^{|\mathcal{G}\setminus\mathcal{B}|},[\bm{q}_{i}]_{i=1}^{|\mathcal{G}\setminus\mathcal{B}|})$ , then the classical second FDT holds in the nonequilibrium steady state:

\displaystyle K(t)=-\frac{\langle f(0),f(t)\rangle_{\rho_{\text{NESS}}}}{\langle u^{2}(0)\rangle_{\rho_{\text{NESS}}}}.

(44)

Proof.

The Kolmogrov operator (42) can be decomposed as

\displaystyle\mathcal{K}=\underbrace{\sum_{b\in\mathcal{B}}\gamma_{b}k_{B}T_{b}\partial_{\bm{p}_{b}}\cdot\partial_{\bm{p}_{b}}}_{\mathcal{D}}-\underbrace{\sum_{b\in\mathcal{B}}\gamma_{b}\bm{p}_{b}\cdot\partial_{\bm{p}_{b}}+\sum_{i\in\mathcal{G}}(\bm{p}_{i}\cdot\partial_{\bm{q}_{i}}-\partial_{\bm{q}_{i}}\cdot\mathcal{H}(\bm{p},\bm{q})\partial_{\bm{p}_{i}})}_{\mathcal{L}}=\mathcal{L}+\mathcal{D}.

It is easy to get the desired result using Theorem 2 and Corollary 2.1. ∎

For the heat conduction model (41), a physically meaningful observable $u=u(\bm{q},\bm{p})$ is the heat flux. If we consider an oscillator chain (one-dimensional case) with symmetric on-site and neighbourhood interaction potential energy. i.e. $U_{i}=U(q_{i})$ and $V_{e}(\delta q_{e})=V(q_{i+1}-q_{i})$ . The local and total heat flux of the system can be defined as [25, 44, 26]:

\displaystyle J_{i}=p_{i}V^{\prime}(q_{i+1}-q_{i}),\quad i\in\mathcal{G}\setminus\mathcal{B},\qquad J_{tot,N}=\sum_{j}p_{i}V^{\prime}(q_{j+1}-q_{j}),\quad j\in\mathcal{G}\setminus\mathcal{B}

(45)

where $J_{i}$ is the local heat flux and $J_{tot,N}$ is the total one with $N=|\mathcal{G}\setminus\mathcal{B}|$ . Proposition 1 ensures the validity of the classical second FDT for observables of the bulk coordinates. This implies that for the local and total heat flux defined as (45), the classical second FDT holds in the NESS even though the explicit form of the steady state probability density $\rho_{NESS}$ is unknown. Note that our definition of the total heat flux excludes the heat flux at the chain boundary. Another frequently used definition of $J_{tot,N}(t)$ contains such contributions and can be decomposed as

\displaystyle J_{tot,N}(t)=J_{\mathcal{G}\setminus\mathcal{B}}(t)+J_{\mathcal{B}}(t),

(46)

where $J_{\mathcal{G}\setminus\mathcal{B}}(t)$ is the bulk contribution as (45) and $J_{\mathcal{B}}(t)$ is the boundary contribution. When applying the generalized second FDT to (46), we have a non-zero additional term $w(0)$ which breaks the classical second FDT. However, with some weak assumptions, we can show (see D) that the dynamics of the averaged heat flux $J_{av}(t)=J_{tot,N}(t)/N$ can be approximated by the bulk contribution $J_{\mathcal{G}\setminus\mathcal{B}}(t)/N$ in the thermodynamic limit as $N\rightarrow\infty$ . Hence we conclude that the GLE (13)-(14) model for the averaged heat flux satisfies the classical second FDT in the thermodynamic limit.

Lastly, we want to comment on the mathematical difficulty to get similar results on the second FDT for nonequilibrium systems generated by deterministic forces. Consider a similar heat conduction chain model driven by Nosé-Hoover thermostats. By assuming that $m_{i}=m=1$ and there are only two thermostats with temperature $T_{L}$ and $T_{R}$ , the dynamics is described by the following equations of motion [25]:

\displaystyle\begin{dcases}\frac{dq_{i}}{dt}&=p_{i}\\ \frac{dp_{i}}{dt}&=-\partial_{q_{i}}\mathcal{H}(\bm{p},\bm{q})-\begin{dcases}\gamma_{L}p_{i},\qquad\text{if $i\in\mathcal{B}_{T_{L}}$}\\ \gamma_{R}p_{i},\qquad\text{if $i\in\mathcal{B}_{T_{R}}$}\end{dcases}\end{dcases},

(47)

where $\mathcal{B}_{L,R}$ are the set of boundary oscillators which interact with thermostats at the temperature $T_{L,R}$ . The cardinality of $\mathcal{B}_{L,R}$ are denoted as $|\mathcal{B}_{L,R}|$ . The dynamics of the auxiliary variables $\gamma_{L,R}$ are given by:

\displaystyle\frac{d\gamma_{L,R}}{dt}=\frac{1}{\theta_{L,R}}\left(\frac{1}{k_{B}T_{L,R}|\mathcal{B}_{T_{L,R}}|}\sum_{n\in\mathcal{B}_{T_{L,R}}}p_{n}^{2}-1\right),

(48)

where $\theta_{L,R}$ are the thermostat response times. It is easy to check that the velocity field for the combined system (47)-(48) has divergence

\displaystyle\nabla\cdot\bm{F}(\bm{p},\bm{q},\bm{\gamma})=-\gamma_{L}(t)-\gamma_{R}(t)

which changes back and forth between positive and negative values depending on the kinetic temperature of the boundary oscillators [25]. As a consequence, the whole system oscillates between energy dissipating state and increasing state. Mathematically speaking, it is very hard to define a proper probability measure to quantify such a nonequilibrium steady state. Even the SRB measure [2, 40], which works for dissipative systems, is not applicable to this case. In practice, one may assume the ergodicity of the deterministic system so that the ensemble average $\langle\cdot\rangle_{\rho_{NESS}}$ can be replaced by the time average. Then due to the similarity of the dynamics, it is reasonable to expect the heat conduction model generated by deterministic thermostats shares some properties with the stochastic model (41), including the classical and the generalized second FDT. This conjecture, however, needs to be verified.

6 Application to reduced-order modeling

In this section, we apply the generalized second FDT to the reduced-order modeling of the heat conduction problem. We first propose suitable methods to approximate the memory kernel $K(t)$ and the fluctuating force $f(t)$ for a GLE satisfying the classical second FDT. The resulting stochastic integro-differential equation serves as the reduced-order model for Gaussian observables of the nonequilibrium system. This model enables us to numerically verify the generalized second FDT. Secondly, we propose a polynomial chaos expansion method to approximate the dynamics of non-Gaussian observables which satisfy the generalized second FDT. Applying these two methods to the averaged heat flux $J_{av}(t)$ leads to dynamical models which characterize the steady state heat transfer in nonequilibirum systems. We note that a similar approach was used in [5] for Hamiltonian systems and the resulting stochastic model is often referred to as the fluctuating heat conduction model.

6.1 Methodology

Without loss of generality, we consider a one-dimensional GLE for scalar observable $u(t)$

\displaystyle\frac{d}{dt}u(t)=\Omega u(t)+\int_{0}^{t}K(t-s)u(s)ds+f_{u}(t).

(49)

In (49), the streaming constant $\Omega$ is easy to obtain using the definition (15b). However, evaluating the memory kernel $K(t)$ and the fluctuation force $f_{u}(t)$ from the first principle is rather challenging since it involves the approximation of the high dimensional orthogonal flow $e^{t\mathcal{Q}\mathcal{K}}$ . To avoid such technical difficulties, here we adopt a data-driven method introduced in [47] to approximate the memory kernel $K(t)$ . As for the fluctuation force $f_{u}(t)$ , it can be approximated by suitable series expansions of a stochastic process. The whole procedure can be described as follows. First of all, we recall that the projected GLE yields (see Section 3) the evolution of the steady state correlation function $C(t)$ of $u(t)$ :

\displaystyle\frac{d}{dt}C_{u}(t)=\Omega C_{u}(t)+\int_{0}^{t}K(t-s)C_{u}(s)ds,\qquad\text{where}\qquad C_{u}(t)=\frac{\langle u(t),u(0)\rangle_{\rho}}{\langle u^{2}(0)\rangle_{\rho}}.

(50)

Hence with equation (50), the memory kernel $K(t)$ can be represented formally using the inverse Laplace transform if we know $C_{u}(t)$ . If we further assume that the dynamics of the observable $u(t)$ is a stationary Gaussian process, since GLE (49) is a linear equation for $u(t)$ , this implies the fluctuation force $f_{u}(t)$ is also a stationary Gaussian process. Then we can use the truncated Karhunen-Loéve (KL) expansion series to approximate $f_{u}(t)$ , namely

\displaystyle f_{u}(t)\simeq\langle f_{u}(t)\rangle_{\rho}+\sum_{k=1}^{K}\xi_{k}\sqrt{\lambda_{k}}e_{k}(t)=\bar{f}_{u}+\sum_{k=1}^{K}\xi_{k}\sqrt{\lambda_{k}}e_{k}(t).

(51)

In the steady state, the mean value of the stochastic process satisfies $\langle f_{u}(t)\rangle_{\rho}=\langle f_{u}(0)\rangle_{\rho}=\bar{f}_{u}$ , which can be obtained by taking the ensemble average of the GLE (49) and then evaluating it at $t=0$ :

\displaystyle\bar{f}_{u}=\langle f(0)\rangle_{\rho}

\displaystyle=\langle u(0)\rangle_{\rho}-\Omega\langle u(0)\rangle_{\rho}=\langle u(0)\rangle_{\rho}-\frac{\langle\dot{u}(0),u(0)\rangle_{\rho}}{\langle u^{2}(0)\rangle_{\rho}}\langle u(0)\rangle_{\rho}.

(52)

The KL expansion random coefficients $\{\xi_{k}\}_{k=1}^{K}$ are necessarily independent Gaussian random variables satisfying $\langle\xi_{i}\xi_{j}\rangle=\delta_{ij}$ , and $\{\lambda_{k},e_{k}\}_{k=1}^{K}$ are, respectively, eigenvalues and eigenfunctions of the homogeneous Fredholm integral equation of the second kind:

\displaystyle\int_{0}^{T}\langle f_{u}(t),f_{u}(0)\rangle_{\rho}e_{k}(s)ds=\lambda_{k}e_{k}(t),\qquad t\in[0,T],

(53)

where $T$ is a certain numerical integration time. If the classical second FDT holds for GLE (49), by substituting $\langle f_{u}(t),f_{u}(0)\rangle_{\rho}=-K(t)\langle u^{2}(0)\rangle_{\rho}$ into eqn (53) and solving for $\{\lambda_{k},e_{k}(t)\}_{k=1}^{K}$ , we can get the exact KL series representation (51) for the fluctuation force $f_{u}(t)$ . With all these terms available, we propose the following data-driven modeling diagram for an arbitrary Gaussian observable $u(t)$ :

\displaystyle\langle u(t),u(0)\rangle_{\rho},\langle u(0)\rangle_{\rho}\Rightarrow\begin{cases}\displaystyle K(t)=\mathfrak{L}^{-1}\left[s-\Omega-\frac{C_{u}(0)}{\tilde{C}_{u}(s)}\right]\\ \displaystyle f_{u}(t)\simeq\bar{f}_{u}+\sum_{k=1}^{K}\xi_{k}\sqrt{\lambda_{k}}e_{k}(t)\end{cases}\Rightarrow\text{Eqn}\ \eqref{eqn:uni_GLE_model},

(54)

where $\mathfrak{L}[C_{u}(t)]=\tilde{C}_{u}(s)$ is the Laplace transform of $C_{u}(t)$ and $\mathfrak{L}^{-1}[\cdot](t)$ is the inverse transform. In flowchart (54), the leftmost three terms are inputs of the diagram which can be obtained by solving numerically the SDE (6) and then averaging samples collected from the steady state simulation. We note that solving (50) for $K(t)$ is a well-known inverse problem which is ill-conditioned. In this paper, we use the series expansion method [47] and LASSO regression to approximate $K(t)$ . By combining all these approximations and executing (54), we can get the first reduced-order model for $u(t)$ :

\displaystyle\frac{d}{dt}u(t)=\Omega u(t)+\sum_{i=1}^{I}\int_{0}^{t}k_{i}b_{i}(t-s)u(s)ds+f(t,\bm{\xi}),

(55)

where $K(t)=\sum_{i=1}^{I}k_{i}b_{i}(t)$ is the series expansion approximation for $K(t)$ and $f(t,\bm{\xi})$ is the truncated KL expansion approximateing the fluctuation force $f_{u}(t)$ . We want to emphasize that the key relation that rationalizes the whole algorithm is the classical second FDT.

The above reduced-order modeling diagram cannot be applied to non-Gaussian cases nor the case where the generalized second FDT holds. The main modeling difficulties stem from the fact that the steady state distribution of the fluctuation force $f(t)$ and the additional term $w(0)$ in the memory kernel $K(t)$ are generally unknown and cannot be easily constructed from MD simulation. However, we can use the polynomial chaos expansion to directly simulate $u(t)$ . To this end, we propose the following modeling diagram for non-Gaussian $u(t)$ :

\displaystyle\langle u(t),u(0)\rangle_{\rho},\rho_{u},\langle u(0)\rangle_{\rho}\Rightarrow\begin{cases}\displaystyle K(t)=\mathfrak{L}^{-1}\left[s-\Omega-\frac{C_{u}(0)}{\tilde{C}_{u}(s)}\right]\\ \displaystyle\text{Eqn}\ \eqref{eqn:PGLE_heat}\end{cases}\Rightarrow u(t)=\sum_{i=1}^{M}u_{i}H_{i}(\gamma(t,\xi)).

(56)

In (56), $u(t)=\sum_{i=1}^{M}u_{i}H_{i}(\gamma(t,\bm{\xi}))$ is the polynomial chaos expansion for a stationary non-Gaussian process, which can be constructed from the time-autocorrelation and the steady state probability density $\rho_{u}$ . Specifically, the expansion coefficient $u_{i}$ and a Gaussian process $\gamma(t,\bm{\xi})$ is calculated via a modified Sakamoto-Ghanem algorithm [41]. In E, we explain the procedure in detail. By directly simulating the non-Gaussian processes in the state space, we avoid the computation of the infinite Kramers-Moyal coefficients [39] or the effective Fokker-Planck diffusion coefficient [5].

For these two reduced-order modeling methods, by approximating the full GLE or using the polynomial chaos expansion method, we can construct a surrogate model for $u(t)$ . Moreover, it also enables us to use a short-time MD simulation data to predict the long-time dynamics of $u(t)$ . This part will be verified later via numerical simulations in the following section.

6.2 Numerical result for a one-dimensional heat conduction model

We now study numerically the one-dimensional heat conduction model (41). In particular, we will use reduced order models introduced in Section 6.1 to build effective models for different phase space observables, from which we can verify the validity of the generalized second FDT and demonstrate the effectiveness of these reduced order models. Moreover, we will study in detail the model for the averaged heat flux $J_{av}(t)$ and discuss its usefulness in characterizing the heat transport intensity for systems in and out of the statistical equilibrium.

To this end, we set the on-site potential energy in (41) to be 0 and the neighbourhood interaction potential energy to be the Lennard-Jones (LJ) potential energy, i.e.:

\displaystyle V_{e}(\delta q_{e})=V(q_{i+1}-q_{i})=4\epsilon\left[\left(\frac{\sigma}{q_{i+1}-q_{i}}\right)^{12}-2\left(\frac{\sigma}{q_{i+1}-q_{i}}\right)^{6}\right].

The whole chain is linked with two thermostats with temperatures $T_{L}$ and $T_{R}$ which will be specified later for different cases. Free boundary conditions are imposed and the modeling parameters are set as follows: $N=|\mathcal{G}|=256$ , $\epsilon=0.2$ , $\sigma=1$ , $\gamma_{L}=\gamma_{R}=1$ . To solve (41) numerically, we use the Euler-Maruyama scheme with step size $dt=10^{-5}$ . In Figure 1, we show the sample trajectories of selected observables of this stochastic system.

6.2.1 Verification of the generalized second FDT

The reduced-order modeling method we introduced in Section 6.1 enables us to numerically verify the generalized second FDT by hypothesis testing. We will first choose an observable for which the classical second FDT holds ( $\Omega=0$ and $w(0)=0$ ) and show that the stochastic model (55) gives correct statistics for observable $u(t)$ . Then we will repeat the procedure for an observable for which the generalized second FDT holds ( $\Omega\neq 0$ and $w(0)\neq 0$ ). Since the algorithm works only when the classical second FDT holds, stochastic model (55) should give wrong statistics for observable $u(t)$ . Throughout this subsection, the thermostat temperatures are set to be $T_{L}=1$ and $T_{R}=5$ which makes the system approaching an NESS as $t\rightarrow+\infty$ .

Observable $p_{M}$

We choose the momentum of the oscillator in midst as the observable. According to Proposition 1, the momentum $p_{M}$ satisfies the classical second FDT (44) in the NESS. Figure 3 shows that $p_{M}$ is a Gaussian variable and the marginal probability density function (PDF) satisfies $p_{M}\sim\mathcal{N}(0.2955,0.42^{2})$ approximately. Therefore the standard KL expansion can be used to represent the the fluctuation term, which leads to the reduced-order model for $p_{M}(t)$ :

\displaystyle\frac{d}{dt}p_{M}(t)=\sum_{i=1}^{I}\int_{0}^{t}k_{i}b_{i}(t-s)p_{M}(s)ds+\sum_{i=1}^{M}\sqrt{\lambda_{i}}\xi_{i}e_{i}(t).

(57)

In this paper, the memory kernel expansion basis $b_{i}(t)$ are set to be the Laguerre polynomials, and a LASSO regression method is used to approximate the expansion coefficient $k_{i}$ [47]. In (57), $I=20$ and $M=500$ (the same hereinafter). Since a stationary Gaussian process $p_{M}(t)$ is fully characterized by its marginal distribution and the time autocorrelation function, we solve (57) numerically and compare these two statistics with the exact results obtained through MD simulation. In Figure 3, we can see that the solution of the reduced-order model (57) reproduces the correct statistics of the observable $p_{M}(t)$ .

Observable $p_{L}$

We choose the momentum of the leftmost oscillator as the observable. According to Proposition 1, the momentum $p_{L}$ satisfies the generalized second FDT (44) ( $w(0)\neq 0$ ) in the NESS. Figure 3 shows that $p_{L}$ is also a Gaussian variable with the marginal PDF satisfying $p_{L}\sim\mathcal{N}(-0.4169,1)$ approximately. If we assume that the classical second FDT holds, then the reduced-order model for $p_{L}(t)$ is given by

\displaystyle\frac{d}{dt}p_{L}(t)=\Omega p_{L}(t)+\sum_{i=1}^{I}\int_{0}^{t}k_{i}b_{i}(t-s)p_{L}(s)ds+\sum_{i=1}^{M}\sqrt{\lambda_{i}}\xi_{i}e_{i}(t),

(58)

Using the data we obtained in the NESS we obtain $\Omega\approx-0.8622$ . When comparing the marginal distribution and the time autocorrelation function with the exact MD simulation results, we find that the reduced-order model reproduces the correct autocorrelation function for observable $p_{L}(t)$ , which is reasonable since the evolution equation for $C_{p_{L}}(t)$ is given by (50) and only depends on the memory kernel $K(t)$ . However, since the classical second FDT does not hold ( $w(0)\neq 0$ ) for the GLE of $p_{L}(t)$ , the reduced-order model with the pre-assumed classical second FDT must reproduce the wrong NESS marginal distribution $\rho_{p_{L}}$ .

Remark

Combining these numerical simulation results for observable $p_{M}(t)$ and $p_{L}(t)$ , we verify the existence of the additional term $w(0)$ in the memory kernel $K(t)$ that violates the classical second FDT. In this paper, we would not determine the specific form of $w(0)$ and leave it as an independent research topic.

6.2.2 Stochastic modeling of the averaged heat flux

In this subsection, we use the reduced-order technique to build dynamical models for the averaged heat flux $J_{av}(t)$ and show that (55) is an effective, generic model for heat transport close to and far-from the statistical equilibrium. In addition, we show that (55) with short MD simulation data can predict the correct long-term dynamics of the averaged heat flux $J_{av}(t)$ .

We firstly briefly review the classical Kubo’s linear response theory for heat transport. For a stochastic model such as (41), the system is initially set to be in the statistical equilibrium (43) and then perturbed by a small temperature difference acted on the boundary. This can be achieved by alternating the temperature of a thermostat as $T_{eq}+\Delta T$ . If $\Delta T$ is sufficiently small, then the thermal conductivity $\kappa$ of the lattice system can be calculated using the first FDT, which is known as the Green-Kubo formula [21, 25, 26, 22]:

\displaystyle\kappa=\frac{N}{k_{B}T_{eq}^{2}}\int_{0}^{\infty}\langle J_{av}(t),J_{av}(0)\rangle_{eq}dt,

(59)

where the ensemble average $\langle\cdot\rangle_{eq}$ is taken with respect to the equilibrium measure (43) with temperature $T_{eq}$ and $N$ is the total particle number. Hence the Green-Kubo formula (59) links the equilibrium time autocorrelation function of the flux with the transport coefficient near the statistical equilibrium. We note that many low-dimensional heat conduction models exhibit a violation of Fourier’s law with $\kappa=\kappa(N)$ depending on $N$ . This, however, will not invalidate the Green-Kubo formula since the $N$ -dependence for such special systems is embedded in the time correlator $\langle J_{av}(t),J_{av}(0)\rangle_{eq}$ . A more difficult case for some low-dimensional systems is the anomalous long-time tail of the heat flux autocorrelation function since it scales as $t^{-d}$ , $0<d<1$ . This anomaly is normally associated with a chain-length dependent conductivity scaling $\kappa(N)\propto N^{\alpha}$ with $\alpha>1$ , which will lead to divergent Green-Kubo integral (59) and an infinite conductivity $\kappa$ . This phenomenon is rather system-dependent and strongly related to the boundary conditions. The LJ system under our investigation does not exhibit such an anomaly, which agrees with the finding in [42] with fixed boundary conditions.

Equilibrium case

We first focus on the equilibrium case and show that the stochastic model we introduced in Section 6.1 yields a correct prediction of the equilibrium time autocorrelation function $\langle J_{av}(t),J_{av}(0)\rangle_{eq}$ . Note that the equilibrium case is a special case of the nonequilibrium model with $T_{L}=T_{R}=T_{eq}$ , and we already proved that the GLE for the averaged heat flux $J_{av}(t)$ in the NESS satisfies the classical second FDT. If the equilibrium distribution for $J_{av}(t)$ is Gaussian, then the full stochastic model can be constructed via (54), which generates a simulated sample trajectory of $J_{av}(t)$ as the solution of (55) with $\Omega=0$ . Figure 4 shows that equilibrium marginal distribution for $J_{av}(t)$ at low temperature $T_{eq}=1$ is approximately Gaussian with $J_{av}(t)\sim\mathcal{N}(0,0.21^{2})$ . Hence the corresponding fluctuation force $f_{J}(t)$ in (55) can be approximated by a truncated KL series, which leads to the reduced-order model:

\displaystyle\frac{d}{dt}J_{av}(t)=\sum_{i=1}^{I}\int_{0}^{t}k_{i}b_{i}(t-s)J_{av}(s)ds+\sum_{i=1}^{M}\sqrt{\lambda_{i}}\xi_{i}e_{i}(t).

(60)

The heat flux has a longer correlation time scale when comparing to observables such as $p_{L},p_{M}$ . When solving (60) numerically, we use a short MD trajectory data for $t\in[0,10]$ to construct the memory kernel $K(t)$ . From Figure 4, we can see that the stochastic model (60) predicts the long time tail of the correlation function.

The situation is different when increasing the system temperature to $T_{eq}=5$ . In particular, the estimated PDF of $J_{av}(t)$ is clearly non-Gaussian and has a long tail. For modeling such heat flux, we adopt the second model (56) where $J_{av}(t)$ is approximated by a polynomial chaos expansion

\displaystyle J_{av}(t)=\sum_{i=1}^{M}J_{i}H_{i}(\gamma(t,\bm{\xi})).

(61)

The simulation result is shown in the second row of Figure 4. We find that with short-term data, the generated stochastic process (61) has a target distribution which agrees with the MD simulation result. Moreover, the extrapolated correlation function of (61) predicates the correct long-time tail of the exact result.

Remark

The reduced-order modeling we introduced for $J_{av}(t)$ is only semi-analytical. Namely, while the second FDT induced GLE is closed and derived from the first principle, the calculation of the memory kernel is extrapolated using a short-term data-driven method. This is less satisfying from a theoretical point of view. Here we want to mention the technical difficulty of developing a pure analytical method to get similar results for the system under our investigation. We note that the time autocorrelation function of the heat flux has multiple time scales. Specifically, there are a short-time scale, $\delta$ -function alike fast decaying of the correlation close to $t=0$ and a long-time scale, slow decaying of the correlation as $t\rightarrow+\infty$ . Generally speaking, it is hard to find an analytical method which generates such two-time-scale dynamics of the correlation simultaneously. Some recently developed analytical methods such as the nonlinear fluctuating hydrodynamics [44, 32] used mode-coupling theory to approximate the memory kernel and successfully predicted the long-time tail of the correlation function. But since the short-time scale dynamics cannot be captured within this framework, which is important for accurately evaluating the Green-Kubo integral (59), it is hard to get the correct transport coefficient based on purely analytical calculations.

Close-to-Equilibrium case

We now consider the close-to equilibrium case and use the classical Green-Kubo theory to verify the validity of stochastic heat conduction model. To this end, we choose $T_{L}=1$ and $T_{R}=1.1$ in the MD simulation of the SDE (41). By collecting the data from the nonequilibrium steady state, we find $J_{av}(t)$ is also approximately Gaussian with $J_{av}\sim\mathcal{N}(0.0170,0.24^{2})$ . Hence the reduced-order model for the close-to-equilibrium averaged heat flux $J_{av}(t)$ is given by

\displaystyle\frac{d}{dt}J_{av}(t)=\sum_{i=1}^{I}\int_{0}^{t}k_{i}b_{i}(t-s)J_{av}(s)ds+\sum_{i=1}^{M}\sqrt{\lambda_{i}}\xi_{i}e_{i}(t).

(62)

In Figure 5, we compare the MD simulation result of the normalized time-autocorrelation function:

\displaystyle C_{J}(t)=\langle J_{av}(t)-\langle J_{av}\rangle_{\Delta T},J_{av}(0)-\langle J_{av}\rangle_{\Delta T}\rangle_{\Delta T},

where $\langle\cdot\rangle_{\Delta T}$ is the ensemble average with respect to the NESS probability density with $T_{R}-T_{L}=\Delta T$ . Similar to the previous case, we used the NESS simulation data for $t\in[0,10]$ to construct the memory kernel. The long-term prediction of the correlation function is shown to be accurate. In addition, we calculate the averaged long-term energy accumulation:

\displaystyle\Delta E(t)=\left\langle\int_{0}^{t}J_{av}(s)ds\right\rangle_{\Delta T},

where $J_{av}(t)$ is the MD sample path for system in the near-equilibrium steady state. For stochastic model (62), this quantity can be calculated by averaging the solution with respect to the probability measure introduced by the i.i.d Gaussian random variables $\{\xi_{i}\}_{i=1}^{M}$ , i.e.

\displaystyle\Delta E(t)=\left\langle\int_{0}^{t}J_{av}(s)ds\right\rangle_{\bm{\xi}=\{\xi_{i}\}_{i=1}^{M}}.

Here $J_{av}(t)$ is the numerical solution of (62). On the other hand, the Green-Kubo formula can also be used to calculate the energy accumulation as:

\displaystyle\Delta E(t)=\kappa\frac{dT}{dx}t\approx\kappa\frac{\Delta T}{aN}t=\frac{t}{k_{B}T^{2}_{eq}a}\int_{0}^{\infty}\langle J_{av}(s),J_{av}(0)\rangle_{eq}ds,

(63)

where $a$ is the equilibrium position spacing of the chain. For our example $a=\sigma=1$ . $\langle J_{av}(s),J_{av}(0)\rangle_{eq}$ can be obtained from the equilibrium MD simulation result or the solution of (60). Note that the second approximation in (63) is valid since normally a linear temperature profile is formed for the heat conduction model [25]. The simulation result verifies that both the Green-Kubo formula and the reduced-order model (62) predicate rather accurately the long-term energy accumulation. In fact, stochastic model (62) gives a different definition of thermal conductivity $\kappa$ based on the second FDT:

\displaystyle\kappa(T,\Delta T,N):=\lim_{t\rightarrow+\infty}a\frac{1}{\Delta T}\frac{\langle\int_{0}^{t}J_{tot,N}(s)ds\rangle_{\Delta T}}{t}\approx\lim_{t\rightarrow+\infty}a\frac{1}{\Delta T}\frac{\langle\int_{0}^{t}J_{tot,N}(s)ds\rangle_{\bm{\xi}}}{t}.

(64)

For our example, $\kappa(T=1,\Delta T=0.1,N=256)\approx\kappa=0.0178$ , where $\kappa$ is the Green-Kubo transport coefficient (59). We also consider a high temperature case. By choosing $T_{L}=5$ and $T_{R}=5.25$ , we found the steady state distribution for $J_{av}(t)$ is similar to the result we obtained for the high-temperature equilibrium case. Because of the non-Gaussian feature of $J_{av}(t)$ , we need to use the aforementioned polynomial chaos expansion method to generate the stochastic model. The simulation result is shown in Figure 5, and we see that the constructed fluctuating heat conduction model faithfully captures and predicates the static and dynamical properties of $J_{av}(t)$ . When comparing the accumulated energy $\Delta E(t)$ , the calculated result is more accurate than the predication of the Green-Kubo formula, with an estimated conductivity (given by (64)) $\kappa(T=5,\Delta T=0.25,N=256)\approx 0.0271$ . This might stems from the fact that for temperature difference $\Delta T=0.25$ , the system is already out of the linear response regime (see explanations below).

Far-from-equilibrium case

When a large temperature gradient $\Delta T$ is imposed to the system through the boundary thermostats, the system is outside of the linear response regime. For such far-from-equilibrium cases, the validity of the Green-Kubo formula (59) is questionable since it is based on perturbation theory. However, we can still use methods introduced in Section 6.1 to build reduced-order models for $J_{av}(t)$ and use it to quantify the intensity of heat transfer. Form the MD simulation result presented in Figure 6, we find that the heat flux probability density in the far-from equilibrium steady state is non-Gaussian, strongly asymmetric and has long tails. Similar results were reported for Hamiltonian systems [32], where the PDF was shown to approximate the Baik-Rains distribution. Here we did not pre-assume the form of the PDF and directly used the numerically evaluated density function to build a model. Specifically, since the heat flux is a non-Gaussian process, we adopt the second methodology (56) and repeat what we have done for the high temperature near-equilibrium case to calculate $J_{av}(t)$ using $J_{av}(t)=\sum_{i=1}^{M}J_{i}H_{i}(\gamma(t,\bm{\xi}))$ . Figure 6 clearly shows that the simulated $J_{av}(t)$ predicts the energy accumulation more accurately than the Green-Kubo formula (with $T_{eq}=1,5$ ), as expected. The estimated conductivity for the two cases we considered are $\kappa(T=1,\Delta T=4,N=256)\approx 0.28$ and $\kappa(T=5,\Delta T=4,N=256)\approx 0.29$ .

Remark

We may compare the classical Green-Kubo transport theory with the second FDT-induced transport theory from the information theory point of view. The Green-Kubo theory uses the equilibrium information (correlation function) of the flux to predict the near-equilibrium transport. In our method, we used the short-time information of flux to predict the long-time transport. These two methods have their own merits and drawbacks. Specifically, the Green-Kubo theory is uniformly valid for equilibrium systems with various perturbations, such as thermostats with different temperature gradients $\Delta T$ . However, it only applies to near-equilibrium systems. The second FDT-induced transport theory has larger range of applicability and can be used to predict far-from-equilibrium transport. But up to this point, we can only reply on data-driven methods to calculate the GLE, which means that the heat transport for systems with different temperature gradients has to be handled on a one-to-one basis.

7 Summary

In this paper, we discussed the generalization of the second fluctuation-dissipation theorem (FDT) for stochastic dynamical systems driven by white noise and its application to the heat conduction model. Following Kubo-Mori’s methodology, the derivation of this new FDT only uses the property of the Komogorov operator $\mathcal{K}$ , and hence generally holds for observables in both equilibrium and nonequilibrium steady states (NESS). We also note a surprising fact that for observables in the degenerate coordinate of the Kolmogorov operator, the classical second FDT holds in the NESS even when the steady state probability measure is generally unknown. The generalized second FDT was applied to various statistical physics models in and out of equilibrium such as the Langevin dynamics and the DPD model. In particular, we focused on a low-dimensional heat conduction model and proved the validity of the classical second FDT for the averaged heat flux. We also introduced two reduced-order modeling methods based on the second FDT. The first one works for Gaussian observables which satisfy the classical second FDT, whereas the other method is based on the polynomial chaos expansion approach which is generally applicable to observables satisfying the generalized second FDT. With these two methods, we were be able to numerically verify our main theoretical result, Theorem 2. Moreover, we derived a fluctuating heat conduction model and introduced a new definition of thermal conductivity which holds for both near-equilibrium and far-from-equilibrium heat transport. For a nanoscale open system out of the linear response regime, this second-FDT-induced conductivity yields more accurate prediction of heat transfer than the classical Green-Kubo theory. We conclude by emphasizing that the presented work can be generalized to other transport processes. Further applications of the proposed methodology can be expected.

Note added in proof. We note that similar conclusions on the validity of the classical second FDT for far-from-equilibrium systems was discovered independently by Jung and Schmid [18] soon after we submitted this manuscript to arXiv. In their work, more numerical evidence from the DPD simulation of an colloid particle immersed in fluids is provided which supports our claims, in particular, the result (II) announced in Section 2. Moreover, a mathematical argument based on the deterministic Volterra equation is shown to provide a different proof of the classical second FDT in the nonequilibrium. In fact, their derivation shows that the necessary condition of Corollary 2.2 is also sufficient which means $\Omega=0$ if and only if $w(0)=0$ . We believe our methodology provides a solid mathematical validation of Jung and Schmid’s result using the stochastic system Mori-Zwanzig equation, and it is constructive to compare the theoretical and numerical results in these two manuscripts.

Appendix A Proof of the generalized first FDT

For the perturbed stochastic system (18) evolving from the steady state $\rho(0)=\rho$ , the evolution equation of the probability density is given by the Fokker-Planck equation:

\displaystyle\frac{\partial}{\partial t}\rho(t)=\mathcal{K}^{*}(t)\rho(t),

(65)

where $\mathcal{K}^{*}(t)$ is the Kolmogorov forward (Fokker-Planck) operator of the form

\displaystyle\mathcal{K}^{*}(t)=\mathcal{K}^{*}+\mathcal{K}^{*}_{ext}(t)=\mathcal{K}^{*}+\left[\sum_{i}\delta G_{i}(t)\tilde{F}_{i}(\bm{x})\partial_{x_{i}}+\sum_{i}\delta H^{2}_{i}(t)\tilde{\sigma}_{i}^{2}(\bm{x})\partial_{x_{i}}^{2}\right]^{*}.

Since the added forces are of the order $O(\delta)$ , which is assumed to be a small quantity, the solution of equation (65) can be written as a perturbative expansion of $\rho(t)$ , i.e.

\displaystyle\rho(t)=\rho+\delta\rho_{1}(t)+\delta^{2}\rho_{2}(t)+\cdots,

(66)

where the steady state solution satisfies $\mathcal{K}^{*}\rho=0$ . By substituting the perturbation series (66) into equation (65) and retaining only the first order term $O(\delta)$ , we get the formal expression of $\delta\rho_{1}(t)$ :

\displaystyle\delta\rho_{1}(t)=\int_{0}^{t}dse^{(t-s)\mathcal{K}^{*}}\mathcal{K}^{*}_{ext}(s)\rho.

Hence in the phase space, for observable $u=u(\bm{x})$ , the change of ensemble average induced by the perturbation can be expressed as:

	$\displaystyle\langle u(t)\rangle_{\rho_{\delta}}-\langle u(t)\rangle_{\rho}$	$\displaystyle=\left\langle\int_{0}^{t}e^{(t-s)\mathcal{K}^{}}\mathcal{K}^{}_{ext}(s)ds\rho u(0)\right\rangle+O(\delta^{2})$
		$\displaystyle=\left\langle\rho\int_{0}^{t}\mathcal{K}_{ext}(s)e^{(t-s)\mathcal{K}}dsu(0)\right\rangle+O(\delta^{2})=\int_{0}^{t}\langle\mathcal{K}_{ext}(s)u(t-s)\rangle_{\rho}ds+O(\delta^{2}),$

where $\rho_{\delta}=\rho_{\delta}(t)$ is the formal solution of the Fokker-Planck equation (65). Using the integration by parts formula to simplify $\langle\mathcal{K}_{ext}(s)u(t-s)\rangle_{\rho}$ , we obtain the generalized linear-response formulas (19)-(20).

Appendix B Proof of the generalized second FDT (32)

Following the proof of Theorem 2, we have

\displaystyle\mathcal{L}+\mathcal{L}_{\rho}^{*}=-\nabla\cdot\bm{F}(\bm{x})-\frac{1}{\rho}\mathcal{L}\rho.

Using integration by parts twice and a shorthand notation $\mathcal{D}=\sum_{i,j=1}^{N}\sigma_{i}(\bm{x})\sigma_{j}(\bm{x})\partial_{x_{i}x_{j}}^{2}=\sum_{i,j=1}^{N}\sigma_{i}\sigma_{j}\partial_{x_{i}x_{j}}^{2}$ , we have

	$\displaystyle\mathcal{D}_{\rho}^{*}$	$\displaystyle=\sum_{i,j=1}^{N}\partial_{x_{i}}\partial_{x_{j}}[\sigma_{i}\sigma_{j}(\cdot)]+\frac{1}{\rho}\partial_{x_{i}}\rho\partial_{x_{j}}[\sigma_{i}\sigma_{j}(\cdot)]+\frac{1}{\rho}\partial_{x_{j}}\rho\partial_{x_{i}}[\sigma_{i}\sigma_{j}(\cdot)]+\frac{1}{\rho}\sigma_{i}\sigma_{j}\partial_{x_{i}x_{j}}^{2}\rho$
		$\displaystyle=\mathcal{D}+\frac{1}{\rho}\mathcal{D}\rho+\frac{1}{\rho}\partial_{x_{i}}\rho\partial_{x_{j}}[\sigma_{i}\sigma_{j}(\cdot)]+\frac{1}{\rho}\partial_{x_{j}}\rho\partial_{x_{i}}[\sigma_{i}\sigma_{j}(\cdot)]+\partial_{x_{i}}[\sigma_{i}\sigma_{j}]\partial_{x_{i}}+\partial_{x_{i}}[\sigma_{i}\sigma_{j}]\partial_{x_{i}}+\partial_{x_{i}x_{j}}^{2}[\sigma_{i}\sigma_{j}]$
		$\displaystyle=\mathcal{D}+\frac{1}{\rho}\mathcal{D}\rho+\mathcal{A},$

where $\mathcal{A}$ is a second-order differential operator. Then it is straightforward to obtain

	$\displaystyle\mathcal{D}_{\rho}^{*}+\mathcal{D}$	$\displaystyle=2\mathcal{D}+\frac{1}{\rho}\mathcal{D}\rho+\mathcal{A},$
	$\displaystyle\mathcal{K}^{*}_{\rho}+\mathcal{K}$	$\displaystyle=-\nabla\cdot\bm{F}(\bm{x})-\frac{1}{\rho}\mathcal{L}\rho+\frac{1}{\rho}\mathcal{D}\rho+2\mathcal{D}+\mathcal{A}.$

Since $\mathcal{K}^{*}\rho=\mathcal{L}^{*}\rho+\mathcal{D}^{*}\rho=0$ , $\mathcal{L}+\mathcal{L}^{*}=-\nabla\cdot\bm{F}(\bm{x})$ and $\mathcal{D}^{*}=\mathcal{D}+\mathcal{B}$ , where operator $\mathcal{B}$ is defined as

\displaystyle\mathcal{B}=\sum_{i,j=1}^{N}\partial_{x_{i}}[\sigma_{i}\sigma_{j}]\partial_{x_{j}}+\partial_{x_{j}}[\sigma_{i}\sigma_{j}]\partial_{x_{i}}+\partial^{2}_{x_{i}x_{j}}[\sigma_{i}\sigma_{j}],

we now obtain

\left.\begin{aligned} \mathcal{L}\rho+\mathcal{L}^{*}\rho&=-\rho\nabla\cdot\bm{F}(\bm{x})\\ \mathcal{D}^{*}\rho&=\mathcal{D}\rho+\mathcal{B}\rho\end{aligned}\right\}\Rightarrow-\frac{1}{\rho}\mathcal{L}\rho+\frac{1}{\rho}\mathcal{D}\rho-\nabla\cdot\bm{F}(\bm{x})=-\frac{1}{\rho}\mathcal{B}\rho.

(67)

By substituting the above equality into $\mathcal{K}_{\rho}^{*}+\mathcal{K}$ , we can get the desired result:

$\displaystyle\mathcal{K}^{*}_{\rho}+\mathcal{K}$	$\displaystyle=-\frac{1}{\rho}\mathcal{B}\rho+2\mathcal{D}+\mathcal{A}$	(68)
	$\displaystyle=2\mathcal{D}+\sum_{i,j=1}^{N}\frac{1}{\rho}\sigma_{i}\sigma_{j}\partial_{x_{i}}\rho\partial_{x_{j}}+\frac{1}{\rho}\sigma_{i}\sigma_{j}\partial_{x_{j}}\rho\partial_{x_{i}}+\partial_{x_{i}}[\sigma_{i}\sigma_{j}]\partial_{x_{j}}+\partial_{x_{j}}[\sigma_{i}\sigma_{j}]\partial_{x_{i}}$
	$\displaystyle=\mathcal{S},$

where $\mathcal{S}$ is defined as (32).

Appendix C The generalized second FDT for nonlinear GLEs

The nonlinear GLEs are normally derived using the Zwanzig-type projection operator [51, 4, 46]:

\displaystyle(\mathcal{P}f)(\hat{\bm{x}})=\frac{\int f(\bm{x})\rho(\bm{x})d\tilde{\bm{x}}}{\int\rho(\bm{x})d\tilde{\bm{x}}},\qquad d\bm{x}=\underbrace{dx_{1}dx_{2},\cdots dx_{m}}_{d\hat{\bm{x}}}\underbrace{dx_{m+1}dx_{m+2},\cdots dx_{N}}_{d\tilde{\bm{x}}}=d\hat{\bm{x}}d\tilde{\bm{x}}.

(69)

Zwanzig’s projection (69) is a conditional expectation operator satisfying $\mathcal{P}^{2}=\mathcal{P}$ . It is also a symmetric operator in $L^{2}(M,\rho)$ , i.e. $\langle\cdot,\mathcal{P}\cdot\rangle_{\rho}=\langle\mathcal{P}\cdot,\cdot\rangle_{\rho}$ . Using (69) to derive the Mori-Zwanzig equation for stochastic system (6), we obtain a nonlinear GLE:

\displaystyle\frac{\partial}{\partial t}\underbrace{e^{t\mathcal{K}}\bm{u}(0)}_{\bm{u}(t)}=\underbrace{e^{t\mathcal{K}}\mathcal{P}\mathcal{K}\bm{u}(0)}_{\bm{F}(\bm{u}(t))}+\underbrace{\int_{0}^{t}e^{s\mathcal{K}}\mathcal{P}\mathcal{K}e^{(t-s)\mathcal{Q}\mathcal{K}}\mathcal{Q}\mathcal{K}\bm{u}(0)ds}_{\int_{0}^{t}\bm{R}(t-s,s)ds}+\underbrace{e^{t\mathcal{Q}\mathcal{K}}\mathcal{Q}\mathcal{K}\bm{u}(0)}_{\bm{f}(t)},

(70)

where $\bm{F}(\bm{u}(t))$ is a nonlinear vector function of $\bm{u}(t)$ . The memory integral $\int_{0}^{t}\bm{R}(t-s,s)ds$ is not directly given by a convolution form such as the one in (13). To simplify this equation, we note that in the the conditional projection (69) is an infinite-rank operator [4, 46], therefore can be formally written as [15]:

\displaystyle\mathcal{P}=\sum_{j=1}^{\infty}\langle\cdot,\psi_{j}(0)\rangle_{\rho}\phi_{j}(0),

where $\{\psi_{j}(0)\}_{j=1}^{\infty}$ and $\{\phi_{j}(0)\}_{j=1}^{\infty}$ are set of phase space functions which are linearly independent in $L^{2}(M,\rho)$ . As a result, the Zwanzig-equation memory kernel is given by:

	$\displaystyle\int_{0}^{t}e^{s\mathcal{K}}\mathcal{P}\mathcal{K}e^{(t-s)\mathcal{Q}\mathcal{K}}\mathcal{Q}\mathcal{K}u_{i}(0)ds$	$\displaystyle=\int_{0}^{t}e^{s\mathcal{K}}\sum_{j=1}^{\infty}\langle\psi_{j}(0),\mathcal{K}e^{(t-s)\mathcal{Q}\mathcal{K}}\mathcal{Q}\mathcal{K}u_{i}(0)\rangle_{\rho}\phi_{j}(0)ds$
		$\displaystyle=\int_{0}^{t}\underbrace{\sum_{j=1}^{\infty}\langle\psi_{j}(0),\mathcal{K}e^{(t-s)\mathcal{Q}\mathcal{K}}\mathcal{Q}\mathcal{K}u_{i}(0)\rangle_{\rho}}_{K_{ij}(t-s)}\phi_{j}(s)ds.$

Hence using the proof of Theorem 2, the memory kernel can be rewritten as:

$\displaystyle K_{ij}(t)=\sum_{j=1}^{\infty}\langle\psi_{j}(0),\mathcal{K}e^{t\mathcal{Q}\mathcal{K}}\mathcal{Q}\mathcal{K}u_{i}(0)\rangle_{\rho}$	$\displaystyle=\sum_{j=1}^{\infty}-\langle\mathcal{K}\psi_{j}(0),\mathcal{K}e^{t\mathcal{Q}\mathcal{K}}\mathcal{Q}\mathcal{K}u_{i}(0)\rangle_{\rho}+\langle\mathcal{S}\psi_{j}(0),\mathcal{K}e^{t\mathcal{Q}\mathcal{K}}\mathcal{Q}\mathcal{K}u_{i}(0)\rangle_{\rho}$
	$\displaystyle=\sum_{j=1}^{\infty}-\langle\mathcal{Q}\mathcal{K}\psi_{j}(0),e^{t\mathcal{Q}\mathcal{K}}\mathcal{Q}\mathcal{K}u_{i}(0)\rangle_{\rho}+\langle\mathcal{S}\psi_{j}(0),e^{t\mathcal{Q}\mathcal{K}}\mathcal{Q}\mathcal{K}u_{i}(0)\rangle_{\rho}$
	$\displaystyle=\sum_{j=1}^{\infty}-\langle g_{j}(0),f_{i}(t)\rangle_{\rho}+\langle h_{j}(0),f_{i}(t)\rangle_{\rho},$	(71)

where $g_{j}(0)=\mathcal{Q}\mathcal{K}\psi_{j}(0)$ and $h_{j}(0)=\mathcal{S}\psi_{j}(0)$ . After all these steps, we obtain a nonlinear GLE:

\displaystyle\frac{\partial}{\partial t}u_{i}(t)=F_{i}(\bm{u}(t))+\sum_{j=1}^{\infty}\int_{0}^{t}K_{ij}(t-s)\phi_{j}(s)ds+f_{i}(t),

(72)

where the memory kernel and the fluctuation force satisfy the generalized second FDT (71). The above derivation is of course quite formal. However such mathematical rigorous discussion clearly indicates for stochastic systems the classical second FDT has to be generalized by including more observables such as $g_{j}(0)$ and $h_{j}(0)$ in the time correlations with the noise $f_{i}(t)$ . At the current stage, it seems difficult or even unnecessary to get the explicit expressions for terms such as $\psi_{j}(0)$ , $\phi_{j}(0)$ and $g_{j}(0)$ , $h_{j}(0)$ . However, nonlinear GLE (72) and the generalized second FDT (71) may be used as the ansatz for proper reduced-order modellings for stochastic systems, such as the ones used in [24].

Appendix D Validity of the classical second FDT for the averaged heat flux

It is sufficient to prove the result for $d=1$ system and the generalization to multi-dimensional cases is direct. The local heat flux at the boundaries have to be handled independently. Since the Langevin forces act on the momentum coordinates of the boundary oscillators, the heat reservoir only exchanges energy with the oscillators through the kinetic part [25]. From this observation we may define the boundary heat flux as the non-Hamiltonian contribution of the time derivative of the kinetic energy i.e.,

\displaystyle J_{b}

\displaystyle=\gamma_{b}k_{B}T_{b}\partial_{p_{b}}^{2}\frac{p_{b}^{2}}{2}-\gamma_{b}p_{b}\partial_{p_{b}}\frac{p_{b}^{2}}{2}=\gamma_{b}k_{B}T_{b}-\gamma_{b}p_{b}^{2},\qquad b\in\mathcal{B}.

In the steady state, a stable kinetic temperature profile is often obtained through numerical experiments and the boundary oscillators admit a kinetic temperature $T_{b}^{k}\neq T_{b}$ [25]. We further assume the marginal distribution for the momentum of the boundary oscillators $p_{b}$ , $b\in\mathcal{B}$ is Gaussian, i.e.

\displaystyle\rho_{p_{b}}=\int\rho\ d\bm{q}d\bm{p}\setminus\{p_{b}\}\propto e^{-\beta_{b}\frac{p_{b}^{2}}{2}},

(73)

where $\bm{p}\setminus\{p_{b}\}$ represents all momenta $\bm{p}$ but $p_{b}$ , $\beta_{b}=1/k_{B}T_{b}^{k}$ and $\rho=\rho_{NESS}$ is the probability density corresponding to the NESS. To be noticed that this assumption is verified numerically in Section 6.2.1 (see Figure 3). In the stationary regime, we get the boundary heat flux ensemble average

	$\displaystyle\langle J_{b}(t)\rangle_{\rho}$	$\displaystyle=\gamma_{b}k_{B}(T_{b}-T_{b}^{k}),$
	$\displaystyle\langle J_{b}^{2}(t)\rangle_{\rho}$	$\displaystyle=\gamma_{b}^{2}k_{B}^{2}(3(T_{b}^{k})^{2}+T^{2}_{b}-2T_{b}^{k}T_{b}).$

Hence for the total heat flux we have

\displaystyle\langle J_{tot,N}(t)\rangle_{\rho}=\langle J_{\mathcal{G}\setminus\mathcal{B}}(t)\rangle_{\rho}+\langle J_{\mathcal{B}}(t)\rangle_{\rho}=\langle J_{\mathcal{G}\setminus\mathcal{B}}(t)\rangle_{\rho}+\sum_{b\in\mathcal{B}}\gamma_{b}k_{B}(T_{b}-T_{b}^{k}),

where $J_{\mathcal{G}\setminus\mathcal{B}}$ and $J_{\mathcal{B}}$ denote respectively the bulk and boundary contributions to the total heat flux. As it is shown in [25], in the stationary regime the average bulk heat flux $\langle J_{\mathcal{G}\setminus\mathcal{B}}(t)\rangle_{\rho}$ is an extensive quantity which scales as the order of volume for the lattice system i.e. $O(L^{d})$ , where $L$ is the length of the lattice system. For fixed coupling constant $\lambda_{b}$ and thermal bath temperature $T_{b}$ , the average boundary heat flux $\langle J_{\mathcal{B}}(t)\rangle_{\rho}$ scales as the order of surface area for the lattice system. i.e. $O(L^{d-1})$ . Naturally since $N\propto L^{d}$ , we have

\displaystyle\lim_{N\rightarrow\infty}\frac{1}{N}\langle J_{tot,N}(t)\rangle_{\rho}=\frac{1}{N}\langle J_{\mathcal{G}\setminus\mathcal{B}}(t)\rangle_{\rho}.

(74)

The above limit can be seen as the first-order moment estimation of the total heat flux in the stationary regime. We also need to verify whether the above estimation hold in the second order. To this end, we use the triangle inequality to get

\displaystyle\|J_{\mathcal{G}\setminus\mathcal{B}}(t)\|_{L^{2}(\rho)}-\|J_{\mathcal{B}}(t)\|_{L^{2}(\rho)}\leq\|J_{tot,N}(t)\|_{L^{2}(\rho)}\leq\|J_{\mathcal{G}\setminus\mathcal{B}}(t)\|_{L^{2}(\rho)}+\|J_{\mathcal{B}}(t)\|_{L^{2}(\rho)}.

(75)

Then by combining another triangle inequality

\displaystyle\|J_{\mathcal{B}}(t)\|_{L^{2}(\rho)}\leq\sum_{b}\|J_{b}(t)\|_{L^{2}(\rho)}=\sum_{b}[\gamma_{b}^{2}k_{B}^{2}(3(T_{b}^{k})^{2}+T^{2}_{b}-2T_{b}^{k}T_{b})]^{1/2},

(76)

we have

\displaystyle\|J_{\mathcal{G}\setminus\mathcal{B}}(t)\|_{L^{2}(\rho)}-\sum_{b}\|J_{b}(t)\|_{L^{2}(\rho)}\leq\|J_{tot,N}(t)\|_{L^{2}(\rho)}\leq\|J_{\mathcal{G}\setminus\mathcal{B}}(t)\|_{L^{2}(\rho)}+\sum_{b}\|J_{b}(t)\|_{L^{2}(\rho)}.

Using the embedding inequality associated with $L^{2}(\rho)\rightarrow L^{1}(\rho)$ and Cauchy-Schwartz inequality, we get that

\displaystyle C|\langle J_{\mathcal{G}\setminus\mathcal{B}}(t)\rangle_{\rho}|\leq C\|J_{\mathcal{G}\setminus\mathcal{B}}(t)\|_{L^{1}(\rho)}\leq\|J_{\mathcal{G}\setminus\mathcal{B}}(t)\|_{L^{2}(\rho)},

(77)

where $C$ is the embedding constant. Inequality (76) also implies $\|J_{\mathcal{B}}(t)\|_{L^{2}(\rho)}$ scales at most of the order $O(L^{d-1})$ , while (77) indicates $\|J_{\mathcal{G}\setminus\mathcal{B}}(t)\|_{L^{2}(\rho)}$ scales at least of the order $O(L^{d})$ . Hence dividing by $N$ and then taking the limit $N\rightarrow\infty$ in both sides of inequality (75), we have

\displaystyle\lim_{N\rightarrow\infty}\frac{1}{N}\|J_{tot,N}(t)\|_{L^{2}(\rho)}=\frac{1}{N}\|J_{\mathcal{G}\setminus\mathcal{B}}(t)\|_{L^{2}(\rho)}.

(78)

Using the similar technique and the stationary condition $\|J_{\mathcal{G}\setminus\mathcal{B}}(t)\|_{L^{2}(\rho)}=\|J_{\mathcal{G}\setminus\mathcal{B}}(0)\|_{L^{2}(\rho)}$ , $\|J_{\mathcal{B}}(t)\|_{L^{2}(\rho)}=\|J_{\mathcal{B}}(0)\|_{L^{2}(\rho)}$ , it is easy to obtain the following estimate

	$\displaystyle\langle J_{\mathcal{G}\setminus\mathcal{B}}(t),J_{\mathcal{G}\setminus\mathcal{B}}(0)\rangle_{\rho}-2$	$\displaystyle\\|J_{\mathcal{G}\setminus\mathcal{B}}(t)\\|_{L^{2}(\rho)}\\|J_{\mathcal{B}}(0)\\|_{L^{2}(\rho)}-\\|J_{\mathcal{B}}(0)\\|^{2}_{L^{2}(\rho)}\leq\langle J_{tot,N}(t),J_{tot,N}(0)\rangle_{\rho}$
		$\displaystyle\leq\langle J_{\mathcal{G}\setminus\mathcal{B}}(t),J_{\mathcal{G}\setminus\mathcal{B}}(0)\rangle_{\rho}+2\\|J_{\mathcal{G}\setminus\mathcal{B}}(t)\\|_{L^{2}(\rho)}\\|J_{\mathcal{B}}(0)\\|_{L^{2}(\rho)}+\\|J_{\mathcal{B}}(0)\\|^{2}_{L^{2}(\rho)}.$		(79)

For the same reason, dividing by $N^{2}$ and then taking the limit $N\rightarrow\infty$ in both sides of inequality (D), we obtain

\displaystyle\lim_{N\rightarrow\infty}\frac{1}{N^{2}}\langle J_{tot,N}(t),J_{tot,N}(0)\rangle_{\rho}=\frac{1}{N^{2}}\langle J_{\mathcal{G}\setminus\mathcal{B}}(t),J_{\mathcal{G}\setminus\mathcal{B}}(0)\rangle_{\rho}.

(80)

Estimates (74) and (80) imply that the averaged total heat flux $J_{av,N}(t)$ , as a stationary second-order process, can be approximated by its bulk contribution $J_{\mathcal{G}\setminus\mathcal{B},N}(t)/N$ in the $L^{2}$ sense.

Appendix E Polynomial chaos expansion for stationary non-Gaussian processes

For strongly non-Gaussian stochastic processes such as the far-from equilibrium heat flux $J_{av}(t)$ , the KL expansion is no longer suitable to represent such processes since the random coefficients $\xi_{k}$ are not i.i.d Gaussian and cannot be easily determined. Some methods were proposed to address the approximation problem of a non-Gaussian process. For instance, Chu and Li [5] suggested to use a Gaussian multiplicative noise to approximate the random force in the extended stochastic system. Zhu and Venturi [49] used Phoon’s algorithm [38, 37] to generate a sample-based, iterated KL expansion to represent the non-Gaussian process. In this paper, we adopt a modified Sakamoto-Ghanem [41] method to approximate $J_{av}(t)$ . The numerical merits of the new algorithm are highlighted in two aspects. First, it is generally appliable to non-Gaussian stochastic process $u(t)$ with arbitrary steady state distribution $\rho_{u}$ and correlation function $C(t_{1},t_{2})=\langle u(t_{1}),u(t_{2})\rangle$ . Secondly, it works much more efficiently when comparing to similar approaches such as Phoon’s algorithm [38, 37], which enables us to rapidly generate sample trajectories from the GLE. Here we only briefly explain the main steps of Sakamoto-Ghanem algorithm and our modification of it. More details and extension of such a method can be found in [41].

Suppose $u(t)$ is a second order, stationary non-Gaussian stochastic process with an arbitrary steady state distribution $\rho_{u}$ and the stationary correlation function $C_{u}(t+s,s)=C_{u}(t,0)=\langle u(t),u(0)\rangle$ . Then the following polynomial chaos expansion approximates $u(t)$ in the $L^{2}$ sense as $M\rightarrow+\infty$ and $K\rightarrow+\infty$ :

\displaystyle u(t)=\sum_{i=1}^{M}u_{i}H_{i}(\gamma(t,\bm{\xi})),\qquad\text{where}\quad\gamma(t,\bm{\xi})=\sum_{j=1}^{K}\sqrt{\lambda_{j}}\xi_{j}e_{j}(t).

(81)

Here $H_{i}$ is the $i$ -th order probabilist’s Hermite polynomials. Being represented as a truncated KL expansion which can be determined by solving the Fredholm equation (53), $\gamma(t,\bm{\xi})$ has steady correlation function $C_{\gamma}(t)=\langle\gamma(t,\bm{\xi}),\gamma(0,\bm{\xi})\rangle$ satisfying the following algebraic equation:

\displaystyle C_{u}(t)=\sum_{j=1}^{M}\frac{u_{j}^{2}j!}{\sum_{k=1}^{M}u_{k}^{2}k!}C^{j}_{\gamma}(t),

(82)

where $u_{j}$ are the coefficients of the Hermite-chaos expansion of the random variable $u$ with probability density $\rho_{u}$ . Specifically, $u_{i}$ can be obtained (see more details in [45], Section 6) using Gaussian quadrature by evaluating the integral:

\displaystyle u_{i}=\frac{1}{\langle H_{i}^{2}\rangle}\int_{0}^{1}U^{-1}(x)H_{i}(G^{-1}(x))dx,

(83)

where $U^{-1}(x)$ and $G^{-1}(x)$ are, respectively, the inverse of the cumulative distribution function (CDF) for an arbitrary random variable $u\sim\rho_{u}$ and a Gaussian random variable $g\sim\mathcal{N}(0,1)$ . In the original version of the Sakamoto-Ghanem algorithm [41], algebraic equation (82) is solved exactly for $C_{u}(t)$ and $C_{\gamma}(t)$ in a discrete lattice $\{0,\Delta t,2\Delta t,\cdots,T\}$ by assuming that $C_{u}(t)>0$ for $t\in[0,T]$ . Here we generalize the algorithm by solving (82) approximately in the same lattice for the part $C_{u}(t)<0$ . i.e. For the lattice point $i\Delta t$ such that $C_{u}(i\Delta t)<0$ , we find an approximated solution $C_{\gamma}(i\Delta t)$ in $[-1,1]$ for the following algebraic equation:

\displaystyle C_{u}(i\Delta t)=\sum_{j=1}^{M}\frac{u_{j}^{2}j!}{\sum_{k=1}^{M}u_{k}^{2}k!}C^{j}_{\gamma}(i\Delta t).

(84)

To apply the modified Sakamoto-Ghanem algorithm in the stochastic modeling of the far-from equilibrium heat flux $J_{av}(t)$ , we only need to input the NESS probability density $\rho_{J}$ and correlation function $C_{J}(t)$ which can be obtained numerically via the MD simulation. We finally obtain the polynomial chaos expansion of $J_{av}(t)$

\displaystyle J_{av}(t)=\sum_{i=1}^{M}J_{i}H_{i}(\gamma(t,\bm{\xi})).

(85)

References

[1] G. S. Agarwal. Fluctuation-dissipation theorems for systems in non-thermal equilibrium and applications. Zeitschrift für Physik A Hadrons and nuclei, 252(1):25–38, 1972.
[2] F. Bonetto, J. L. Lebowitz, and L. Rey-Bellet. Fourier’s law: a challenge to theorists. In Mathematical physics 2000, pages 128–150. World Scientific, 2000.
[3] S. Dal Cengio, D. Levis, and I. Pagonabarraga. Linear response theory and Green-Kubo relations for active matter. Phys. Rev. Lett., 123(23):238003, 2019.
[4] A. J. Chorin, O. H. Hald, and R. Kupferman. Optimal prediction with memory. Physica D, 166(3-4):239–257, 2002.
[5] W. Chu and X. Li. The Mori–Zwanzig formalism for the derivation of a fluctuating heat conduction model from molecular dynamics. Commun Math Sci, 17(2), 2019.
[6] N. Cuneo, J-P Eckmann, M. Hairer, and L Rey-Bellet. Non-equilibrium steady states for networks of oscillators. Electron. J. Probab., 23, 2018.
[7] A. Dhar. Heat transport in low-dimensional systems. Advances in Physics, 57(5):457–537, 2008.
[8] J. P. Eckmann and M. Hairer. Non-equilibrium statistical mechanics of strongly anharmonic chains of oscillators. Commun. Math. Phys., 212(1):105–164, 2000.
[9] J. P. Eckmann, C. A. Pillet, and L. Rey-Bellet. Non-equilibrium statistical mechanics of anharmonic chains coupled to two heat baths at different temperatures. Commun. Math. Phys., 201(3):657–697, 1999.
[10] P. Español. Hydrodynamics from dissipative particle dynamics. Phys. Rev. E, 52(2):1734, 1995.
[11] P. Español and P. Warren. Statistical mechanics of dissipative particle dynamics. EPL, 30(4):191, 1995.
[12] A. Gritsun, G. Branstator, and A. Majda. Climate response of linear and quadratic functionals using the fluctuation–dissipation theorem. J. Atmos. Sci, 65(9):2824–2841, 2008.
[13] F. Grogan, H. Lei, X. Li, and N. A. Baker. Data-driven molecular modeling with the generalized Langevin equation. J. Comput. Phys., 418:109633–109641, 2020.
[14] M. Hairer and A. J. Majda. A simple framework to justify linear response theory. Nonlinearity, 23(4):909, 2010.
[15] P. D. Hislop and I. S. Sigal. Introduction to spectral theory: With applications to Schrödinger operators, volume 113. Springer Science & Business Media, 2012.
[16] P. J. Hoogerbrugge and J. M. V. A. Koelman. Simulating microscopic hydrodynamic phenomena with dissipative particle dynamics. EPL, 19(3):155–160, 1992.
[17] T. Hudson and X. H. Li. Coarse-graining of overdamped Langevin dynamics via the Mori–Zwanzig formalism. Multiscale Modeling & Simulation, 18(2):1113–1135, 2020.
[18] G. Jung and F. Schmid. Fluctuation-Dissipation eelations far from equilibrium: A case study. arXiv preprint arXiv:2106.00818, 2021.
[19] P. E. Kloeden and E. Platen. Numerical solution of stochastic differential equations, volume 23. Springer Science & Business Media, 2013.
[20] R. Kubo. The fluctuation-dissipation theorem. Reports on progress in physics, 29(1):255, 1966.
[21] R. Kubo, M. Toda, and N. Hashitsume. Statistical physics II: nonequilibrium statistical mechanics, volume 31. Springer Science & Business Media, 2012.
[22] A. Kundu, A. Dhar, and O. Narayan. The Green–Kubo formula for heat conduction in open systems. J. Stat. Mech. Theory Exp., 2009(03):L03001, 2009.
[23] H. Kunita. Stochastic flows and stochastic differential equations. Cambridge university press, 1997.
[24] H. Lei, N.A. Baker, and X. Li. Data-driven parameterization of the generalized Langevin equation. Proc. Natl. Acad. Sci., 113(50):14183–14188, 2016.
[25] S. Lepri, R. Livi, and A. Politi. Thermal conduction in classical low-dimensional lattices. Physics reports, 377(1):1–80, 2003.
[26] S. Lepri, R. Roberto, and A. Politi. On the anomalous thermal conductivity of one-dimensional lattices. EPL (Europhysics Letters), 43(3):271, 1998.
[27] Z. Li, X. Bian, X. Li, and G. E. Karniadakis. Incorporation of memory effects in coarse-grained modeling via the Mori-Zwanzig formalism. J. Chem. Phys, 143:243128, 2015.
[28] C. Maes. On the second fluctuation–dissipation theorem for nonequilibrium baths. J. Stat. Phys, 154(3):705–722, 2014.
[29] A. Majda and B. Gershgorin. Quantifying uncertainty in climate change science through empirical information theory. Proc. Natl. Acad. Sci., 107(34):14958–14963, 2010.
[30] A. Majda, A. V. Rafail, and M. J. Grote. Information theory and stochastics for multiscale nonlinear systems, volume 25. American Mathematical Soc., 2005.
[31] J. C. Mattingly, A. M. Stuart, and D. J. Higham. Ergodicity for SDEs and approximations: locally Lipschitz vector fields and degenerate noise. Stochastic processes and their applications, 101(2):185–232, 2002.
[32] C. B. Mendl and H. Spohn. Current fluctuations for anharmonic chains in thermal equilibrium. J. Stat. Mech. Theory Exp., 2015(3):P03007, 2015.
[33] H. Mori. Transport, collective motion, and Brownian motion. Prog. Theor. Phys., 33(3):423–455, 1965.
[34] T. Morita, H. Mori, and K Mashiyama. Contraction of state variables in non-equilibrium open systems. II. Prog. Theor. Phys, 64(2):500–521, 1980.
[35] G. P. Morriss and D. J. Evans. Statistical Mechanics of Nonequilbrium Liquids. ANU Press, 2013.
[36] G. A. Pavliotis. Stochastic processes and applications: diffusion processes, the Fokker-Planck and Langevin equations, volume 60. Springer, 2014.
[37] K.K. Phoon, H.W. Huang, and S.T. Quek. Simulation of strongly non-Gaussian processes using Karhunen–Loeve expansion. Prob. Eng. Mech., 20(2):188–198, 2005.
[38] K.K. Phoon, S.P. Huang, and S.T. Quek. Simulation of second-order processes using Karhunen–Loeve expansion. Computers & structures, 80(12):1049–1060, 2002.
[39] H. Risken. The Fokker-Planck equation: methods of solution and applications. Springer-Verlag, second edition, 1989. Mathematics in science and engineering, vol. 60.
[40] D. Ruelle. Smooth dynamics and new theoretical ideas in nonequilibrium statistical mechanics. J. Stat. Phys., 95(1):393–468, 1999.
[41] S. Sakamoto and R. Ghanem. Polynomial chaos decomposition for the simulation of non-Gaussian nonstationary stochastic processes. Journal of engineering mechanics, 128(2):190–201, 2002.
[42] A. V. Savin and Y. A. Kosevich. Thermal conductivity of molecular chains with asymmetric potentials of pair interactions. Phys. Rev. E, 89(3):032102, 2014.
[43] I. Snook. The Langevin and generalised Langevin approach to the dynamics of atomic, polymeric and colloidal systems. Elsevier, 2006.
[44] H. Spohn. Fluctuating hydrodynamics approach to equilibrium time correlations for anharmonic chains. In Thermal Transport in Low Dimensions, pages 107–158. Springer, 2016.
[45] D. Xiu and G. Karniadakis. The Wiener–Askey polynomial chaos for stochastic differential equations. SIAM journal on scientific computing, 24(2):619–644, 2002.
[46] Y. Zhu, J. M. Dominy, and D. Venturi. On the estimation of the Mori–Zwanzig memory integral. J. Math. Phys, 59(10):103501, 2018.
[47] Y. Zhu and H. Lei. Effective Mori-Zwanzig equation for the reduced-order modeling of stochastic systems. arXiv preprint arXiv:2102.01377, 2021.
[48] Y. Zhu and D. Venturi. Faber approximation of the Mori–Zwanzig equation. J. Comp. Phys., (372):694–718, 2018.
[49] Y. Zhu and D. Venturi. Generalized Langevin equations for systems with local interactions. J. Stat. Phys., 2020.
[50] Y. Zhu and D. Venturi. Hypoellipticity and the Mori-Zwanzig formulation of stochastic differential equations. arXiv preprint arXiv:2001.04565, 2020.
[51] R. Zwanzig. Memory effects in irreversible thermodynamics. Phys. Rev, 124(4):983, 1961.
[52] R. Zwanzig. Nonlinear generalized Langevin equations. J. Stat. Phys., 9(3):215–220, 1973.

Generalized second fluctuation-dissipation theorem in the nonequilibrium steady state: Theory and applications

Abstract

1 Introduction

2 The main theoretical results

3 Mori-type generalized Langevin equations for SDEs

4 Generalized fluctuation-dissipation theorem

4.1 Derivation of the first FDT

Theorem 1.

4.2 Derivation of the second FDT

Theorem 2.

Proof.

Corollary 2.1.

Corollary 2.2.

Proof.

Remark

5 The generalized second FDT for specific systems

5.1 The generalized second FDT for equilibrium systems

Langevin dynamics

Remark

Dissipative particle dynamics

5.2 The generalized second FDT for nonequilibrium systems

Proposition 1.

Proof.

6 Application to reduced-order modeling

6.1 Methodology

6.2 Numerical result for a one-dimensional heat conduction model

6.2.1 Verification of the generalized second FDT

Observable pMp_{M}

Observable pLp_{L}

Remark

6.2.2 Stochastic modeling of the averaged heat flux

Equilibrium case

Remark

Close-to-Equilibrium case

Far-from-equilibrium case

Remark

7 Summary

Appendix A Proof of the generalized first FDT

Appendix B Proof of the generalized second FDT (32)

Appendix C The generalized second FDT for nonlinear GLEs

Appendix D Validity of the classical second FDT for the averaged heat flux

Appendix E Polynomial chaos expansion for stationary non-Gaussian processes

References

Observable $p_{M}$

Observable $p_{L}$