On the "generalized Generalized Langevin Equation"

TL;DR

This paper extends the generalized Langevin equation to non-stationary, non-equilibrium systems using time-dependent projection operators, enabling the derivation of equations of motion from molecular dynamics data.

Contribution

It introduces a formalism for deriving non-stationary Langevin equations with time-dependent kernels, linking them to simulation data and initial conditions.

Findings

Derived a non-stationary Langevin equation with a time-dependent memory kernel.

Established a fluctuation-dissipation-like relation for non-equilibrium conditions.

Validated the approach with Brownian motion simulations.

Abstract

In molecular dynamics simulations and single molecule experiments, observables are usually measured along dynamic trajectories and then averaged over an ensemble ("bundle") of trajectories. Under stationary conditions, the time-evolution of such averages is described by the generalized Langevin equation. In contrast, if the dynamics is not stationary, it is not a priori clear which form the equation of motion for an averaged observable has. We employ the formalism of time-dependent projection operator techniques to derive the equation of motion for a non-equilibrium trajectory-averaged observable as well as for its non-stationary auto-correlation function. The equation is similar in structure to the generalized Langevin equation, but exhibits a time-dependent memory kernel as well as a fluctuating force that implicitly depends on the initial conditions of the process. We also derive a relation between this memory kernel and the autocorrelation function of the fluctuating force that has a structure similar to a fluctuation-dissipation relation. In addition, we show how the choice of the projection operator allows to relate the Taylor expansion of the memory kernel to data that is accessible in MD simulations and experiments, thus allowing to construct the equation of motion. As a numerical example, the procedure is applied to Brownian motion initialized in non-equilibrium conditions, and is shown to be consistent with direct measurements from simulations.

Equations

$$\frac{d{A}_t}{dt}=\omega A_t+{\int }_0^{t}K(t-\tau )A_{\tau }d\tau +\eta (t)$$

$$\frac{dC(t)}{dt}=\omega C(t)+{\int }_0^{t}K(t-\tau )C(\tau )d\tau$$

$$\frac{d{A}_t}{dt}=i\mathcal{L}{A}_t$$

$$\frac{d{A}_t}{dt}=e^{i\mathcal{L}t}i\mathcal{L}{A}_0$$

$$P_{t^{\prime }}{P}_t=P_t$$

$$P_t{\dot{P}}_t$$

$${\dot{P}}_t{P}_t$$

$$P_t{\dot{P}}_t(1-P_t)={\dot{P}}_t$$

$${\dot{Z}}_t=e^{i\mathcal{L}t}i\mathcal{L}(1-P_t)-e^{i\mathcal{L}t}{\dot{P}}_t$$

$${\dot{Z}}_t=Z_{t}i\mathcal{L}(1-P_t)+e^{i\mathcal{L}t}{P}_t(i\mathcal{L}-{\dot{P}}_t)(1-P_t)$$

$$Z_t$$

$$+{\int }_s^{t}d\tau e^{i\mathcal{L}\tau }\left[i\mathcal{L}-{\dot{P}}_{\tau }\right]\left[1-P_{\tau }\right]G_{\tau ,t}$$

$$G_{\tau ,t}={\exp }_{-}[{\int }_{\tau }^{t}d{t}^{\prime }i\mathcal{L}(1-P_{t^{\prime }})]\equiv 1+\sum _{n=1}^{\infty }{\int }_{\tau }^{t}d{t}_1{\int }_{\tau }^{t_1}d{t}_2\cdots {\int }_{\tau }^{t_{n-1}}d{t}_{n}i\mathcal{L}(1-P_{t_n})\cdots i\mathcal{L}(1-P_{t_2})i\mathcal{L}(1-P_{t_1})$$

$$\frac{d{A}_t}{dt}=e^{i\mathcal{L}t}{P}_{t}i\mathcal{L}{A}_0+Z_{t}i\mathcal{L}{A}_0$$

$$\frac{d{A}_t}{dt}=$$

$$+{\int }_s^{t}d\tau e^{i\mathcal{L}\tau }{P}_{\tau }\left[i\mathcal{L}-{\dot{P}}_{\tau }\right]\left[1-P_{\tau }\right]G_{\tau ,t}i\mathcal{L}{A}_0$$

$$+e^{i\mathcal{L}s}\left[1-P_s\right]G_{s,t}i\mathcal{L}{A}_0$$

$$P_t^{\mathrm{b}}F:=\frac{⟨A_t^{\ast }{F}_t⟩}{⟨|A_t{|}^2⟩}{A}_0$$

$$⟨X_t⟩=\int d{\mathbf{\Omega }}_0\rho ({\mathbf{\Omega }}_0)e^{i\mathcal{L}t}X({\mathbf{\Omega }}_0)$$

$$⟨X_t{Y}_{t^{\prime }}⟩=\int d{\mathbf{\Omega }}_0\rho ({\mathbf{\Omega }}_0)\left[e^{i\mathcal{L}{t}^{\prime }}Y({\mathbf{\Omega }}_0)\right]\left[e^{i\mathcal{L}t}X({\mathbf{\Omega }}_0)\right]$$

$${\omega }_n(t)=\frac{⟨A_t^{\ast }{A}_t^{(n)}⟩}{⟨|A_t{|}^2⟩}$$

$$K(t,\tau )=$$

$${\dot{P}}_{\tau }^{b}F=\frac{⟨A_{\tau }^{\ast }i\mathcal{L}{e}^{i\mathcal{L}\tau }[1-P_{\tau }^b]F⟩+⟨e^{i\mathcal{L}\tau }[1-P_{\tau }^b]Fi\mathcal{L}{A}_{\tau }^{\ast }⟩}{⟨|A_{\tau }{|}^2⟩}{A}_0$$

$$K(t,\tau )=$$

$$K(t,\tau )$$

$$+\sum _{n=1}^{\infty }{\int }_{\tau }^{t}d{t}_1...{\int }_{\tau }^{t_{n-1}}d{t}_n{k}_n(\tau ,t_1,...,t_n)$$

$$k_n(\tau ,t_1,...,t_n)=$$

$$i\mathcal{L}\left[1-P_{t_n}^{\mathrm{b}}\right]...i\mathcal{L}\left[1-P_{t_1}^{\mathrm{b}}\right]i\mathcal{L}{A}_0⟩{⟨|A_0{|}^2⟩}^{-1}$$

$$\frac{d{A}_t}{dt}={\omega }_1(t)A_t+{\int }_s^{t}d\tau K(t,\tau )A_{\tau }+{\eta }_s(s,t)$$

$${\eta }_s(s,t)=e^{i\mathcal{L}s}\left[1-P_s^{\mathrm{b}}\right]G_{s,t}i\mathcal{L}{A}_0$$

$$⟨{\eta }^{MZ}(t)⟩=0$$

$$⟨{\eta }^{MZ}(t){\eta }^{MZ}(t^{\prime })⟩=-K^{MZ}(t-t^{\prime }){⟨|A{|}^2⟩}_{\mathrm{eq}}$$

$$⟨A_s^{\ast }{e}^{i\mathcal{L}s}\left[1-P_s^{\mathrm{b}}\right]F⟩=⟨A_s^{\ast }{F}_s⟩-⟨A_s^{\ast }\frac{⟨A_s^{\ast }{F}_s⟩}{⟨|A_s{|}^2⟩}{A}_s⟩=0$$

Full text

.1 Introduction

Consider a system of $N$ classical particles that evolve according to Hamilton’s equations of motion, and a phase space variable $A_{t}\equiv A(\mathbf{q}^{N}(t),\mathbf{p}^{N}(t))$, where $\mathbf{q}^{N}(t)$ and $\mathbf{p}^{N}(t)$ are the positions and momenta of the particles. In several sub-fields of statistical physics quantities such as $\langle A_{t}A_{t^{\prime}}^{*}\rangle_{\rm Trajectories}$ are studied, i.e. averages of auto-correlations in $A$ taken over bundles of trajectories. (If the reader wonders why we introduce the auto-correlation rather than the simpler expression $\langle A_{t}\rangle_{\rm Trajectories}$, which is of equal practical importance, we suggest to skip ahead and compare eqns. (24) and (29).  Note that the latter expression is easier to analyze, because the fluctuating force averages out.)

One context, in which this type of quantity is relevant, is the field of molecular dynamics simulation of microscopic systems with the aim of constructing coarse-grained models: for instance, the system could be a polymer melt and the aim could be to develop a rheological model padding:11 ; or the system could be a biomolecule, which a researcher might simulate with classical atomistic force-fields and monitor the collective motion of specific groups of atoms in order to deduce a simplified model of a biological mechanism berendsen:07 ; kmiecik:16 . A different and equally important context is the experimental study of non-equilibrium work-relations in single molecule experiments harris:04 ; hoffmann:12 , where $A_{t}$ would be e.g. the extension of a piece of DNA or of a protein under an applied force.

In any of these contexts, it is useful to have general information on the properties of the equation of motion that governs $\langle A_{t}A_{t^{\prime}}^{*}\rangle_{\rm Trajectories}$. Therefore, in this article we discuss the form of the equation of motion in the general (i.e. the non-stationary) case.

The problem of integrating out a large number of degrees of freedom of a thermodynamic system can be treated by several approaches. Most widely used is the framework of stochastic equations of motion, which has been introduced by Langevin langevin:1908 with his description of Brownian motion. In the 1960’s a formalism was developed by Zwanzig zwanzig:1961 and Mori mori:1965 to account for memory effects in non-trivial systems. This formalism is based on the definition of projection operators that are aimed at integrating out the fast dynamics of a process in order to study slow variables.

Since then, several approaches were introduced to extend the formalism to non-stationary dynamics kawasaki:1973 ; willis:1974 ; furukawa:1979 ; grabert:1982 ; linden:1998 ; fuchs:2009 . We will, in the following base our arguments on Grabert’s approach grabert:1982 and introduce a new projection operator, that is particularly suited to study variables of the type $\langle A_{t}A^{*}_{t^{\prime}}\rangle$. (To shorten the notation, we have dropped the subscript “Trajectories”. Unless stated otherwise, in the following all averages are meant as trajectory averages.) Note that the relations discussed here do not require time-scale separation, i.e. the averaged observable can be any variable of interest, regardless of whether it is ’slow’ or ’fast’, ’relevant’ or not.

.2 The stationary case

Before we address the non-stationary case, let us briefly recall the structure of the equations in equilibrium or other stationary situations. The evolution of the the variable $A_{t}$ is governed by the equation

[TABLE]

where $\omega$ is a drift coefficient, $K(t)$ is a so-called memory kernel, and $\eta(t)$ is a so-called fluctuating force hansen:1990 . (We denote the time-dependence of operators by a subscript in contrast to the time-dependence of functions, which we denote in brackets). This equation has the form of a generalized Langevin equation (GLE). The corresponding auto-correlation function $C(t)=\left<A_{t}A_{0}^{*}\right>=\left<A_{t+t^{\prime}}A_{t^{\prime}}^{*}\right>$ evolves according to

[TABLE]

The dependence of $K$ on $t-\tau$ in the integral is convenient because the convolution theorem can be used to Laplace transform equation (2) mokshin:2005 . Thus if the autocorrelation function $C(t)$ and the drift coefficient $\alpha$ can be obtained in an experiment or simulation, the stationary memory kernel can be constructed. Note that under non-stationary conditions this property does not hold, i.e. one can in general not obtain the memory kernel by a simple Laplace transform of the observed dynamics of a coarse-grained variable.

.3 Time-dependent projection operators: a brief reminder

To introduce the non-stationary case, we briefly recall Grabert’s approach grabert:1982 . Consider the time-evolution of the dynamical variable $A_{t}$

[TABLE]

where $i\mathcal{L}$ is the propagator, e.g. the Liouvillian operator in the case of Hamiltonian dynamics. (Note that the following arguments are not restricted to Hamiltonian dynamics. They hold for any dynamics that can be described by an equation of the form of eqn. (3). Thus our conclusions also apply, in particular, to simulations with a Nosé-Hoover thermostat.) Equation (3) can be formally integrated and then again differentiated with respect to time to yield

[TABLE]

In the equilibirum Mori-Zwanzig formalism, one usually defines a stationary projection operator which is used to split the dynamics into a parallel (slow) part and on orthogonal (fast) part. When extending the formalism to non-stationary processes, Grabert introduced the time dependence directly in the projection operator. Assume that one can define a time-dependent operator $P_{t}$ that acts on phase space variables and that fullfills

[TABLE]

for all times $t$ and $t^{\prime}$. Note that if $t^{\prime}=t$, one obtains $P_{t}^{2}=P_{t}$, i.e. this operator is a projector. In fact, it indicates that $P_{t}$ projects onto a fixed subspace for all $t$, but the orientation of the projection might change with $t$. This implies that the projection of a vector with respect to a certain $t$ is part of the fixed subspace and thus remains constant once projected with respect to another time $t^{\prime}$. We can then take the derivatives with respect to either $t$ or $t^{\prime}$, and then take the limit $t^{\prime}\rightarrow t$ to find

[TABLE]

Thus, we can write

[TABLE]

Now we define the operator $Z_{t}=e^{i\mathcal{L}t}(1-P_{t})$. Its time-derivative is

[TABLE]

Inserting $1=\mathcal{P}_{t}+(1-\mathcal{P}_{t})$ in the first term and using eqn. (8) yields

[TABLE]

This differential equation for $Z_{t}$ can be solved using time-ordered exponentials,

[TABLE]

which is valid for any reference time $s\leq t$. In this expression $G_{\tau,t}$ is the negatively time-ordered exponential operator, i.e. the unique solution of the differential equation $dY_{\tau,t}/dt=Y_{\tau,t}i\mathcal{L}(1-P_{t})$ with initial condition $Y_{t,t}=1$. It can be written as

[TABLE]

We now split the propagator $e^{i\mathcal{L}t}$ into $e^{i\mathcal{L}t}P_{t}+e^{i\mathcal{L}t}(1-P_{t})$, which allows us to rewrite eqn. (4)

[TABLE]

By making use of eqn. (.3), we obtain the following equation of motion for $A_{t}$

[TABLE]

This equation is valid for any projection operator as long as it satisfies the identity (5). Here we close the reminder of Grabert’s work and come to new aspects.

.4 A time-dependent projection operator for bundles of trajectories

We specify a particular projector $P^{b}_{t}$ by defining its action on any function of phase space $F$. As we intend to apply the technique to bundles of trajectories – be they generated by a set of experiments or of molecular dynamics simulations – we introduce a definition that is natural in this context (this is where our work differs from previous work on time-dependent projection operators). In an MD simulation, one typically initializes a bundle of trajectories at a given distribution of points in phase space, $\rho(\mathbf{\Omega}_{0})$, then one numerically propagates them and computes the variable of interest $A$ on each trajectory at certain times $t$. Finally one takes the average of $A_{t}A^{*}_{0}$ over all simulated trajectories. (A set of experiments is carried out and analyzed in exactly the same way, although $\rho(\mathbf{\Omega}_{0})$ can usually not be prescribed.)

We thus define a time-dependent projection operator $P^{\rm b}_{t}$ by its action on a dynamical variable $F$.

[TABLE]

where the brackets mean an average over all possible trajectories (indicated by the superscript b for “bundle”) starting from a well-defined distribution of initial configurations, i.e.

[TABLE]

This average is also used to define the correlation function between two dynamical variables $X$ and $Y$ by

[TABLE]

Here, the exponential operators are understood to act only on arguments inside the enclosing square brackets. (In later expressions, operators are understood as acting on all arguments that appear to their right inside the angular brackets that denote subsequent averaging.)

Compare $P^{\rm b}_{t}$ to the projection operator which is used in the stationary Mori-Zwanzig formalism, $P^{\rm MZ}F=\left\langle AF\right\rangle_{\rm eq}\left\langle A^{2}\right\rangle_{\rm eq}^{-1}A_{0}$, where $\left\langle\cdots\right\rangle_{\rm eq}$ stands for the equilibrium ensemble average. In the case of equilibrium processes, our projector collapses with $P^{\rm MZ}$. Note that the property (5) is satisfied. $P^{\rm b}_{t}$ can be seen as an operator that projects onto the fixed vector $A$ but whose orientation of projection evolves with time (see fig. (1)).

We now apply definition (15) to eqn. (.3). To simplify the resulting equation, we introduce the following notation :

[TABLE]

where $A^{(n)}_{t}$ stands for the $n$-th time derivative of $A_{t}$, i.e. $A^{(n)}_{t}\equiv(i\mathcal{L})^{n}A_{t}$. Using this notation, the first term of eqn. (.3) becomes $\omega_{1}(t)A_{t}$. The second term requires more attention. In eqn. (.3) the integrand has on its left side the operator $e^{i\mathcal{L}\tau}P_{\tau}$, which implies that it is proportional $A_{\tau}$. Thus, it becomes of the form $K(t,\tau)A_{\tau}$, where

[TABLE]

This expression further simplifies once one evaluates $\dot{P}_{\tau}^{b}$ explicity. There holds (see appendix)

[TABLE]

Inserting this into eq. (19), one obtains

[TABLE]

Interestingly, the derivative of the projector vanishes in favor of shifting the Liouville operator in $[i\mathcal{L}-\dot{P}_{\tau}^{b}]$ to act “to the left” (an operation that superficially resembles taking the adjoint with respect to the weighted scalar product).

The definition of $G_{\tau,t}$ allows to write

[TABLE]

The objects $k_{n}$ are functions of the instantaneous correlations of $A$ and its first $n+2$ time-derivatives. They are defined as

[TABLE]

Combining these terms, we obtain an equation of motion for $A_{t}$ :

[TABLE]

with

[TABLE]

The structure of this equation resembles the Generalized Langevin Equation (GLE), however, as the system is not in a steady state, the friction kernel $K(t,\tau)$ does not necessarily depend on only $t-\tau$.

Note that in the stationary Mori-Zwanzig case, eq. (24) is covariant under an arbitrary time translation $s\mapsto s+t_{0}$, $t\mapsto t+t_{0}$, and the first term in eq. (25) guarantees this covariance for $\eta_{s}(s,t)$. One can then set $s=0$ to recover the term that is usually identified as a noise kubo:1966 with

[TABLE]

and

[TABLE]

In the present case, the situation is more complex. From the definition of $P^{\rm b}_{t}$ follows that $\eta_{s}(s,t)$ is perpendicular to $A_{s}$, i.e. $\left\langle A^{*}_{s}\eta_{s}(s,t)\right\rangle=0$, but not to $A_{T}$ for $T\neq s$. The orthogonality can be seen by calculating $\left\langle A^{*}_{s}e^{i\mathcal{L}s}\left[1-P^{\rm b}_{s}\right]F\right\rangle$ with $F=G_{s,t}i\mathcal{L}A_{0}$, i.e.

[TABLE]

This allows us to write an equation of motion for the two-time auto-correlation function $C(s,t)=\left\langle A^{*}_{s}A_{t}\right\rangle$ by multiplying eq. (24) by $A^{*}_{s}$ and taking the trajectory average defined in (16)

[TABLE]

This equation, just as in the stationary case of eqn. (2), does not contain the fluctuating force $\eta_{s}(s,t)$ anymore, but its form respects the non-stationarity of the problem. If one shifts the time origin by an amount $t_{0}$, i.e. one considers the equation for $C(s+t_{0},t+t_{0})$, its solution will be a priori different because $K(t,\tau)\neq K(t+t_{0},\tau+t_{0})$. Moreover, note that it is not possible to simply Laplace transform eqn. (29), because the convolution theroem does not apply anymore.

Eqn. (24) and eqn. (29) are the central results of our work. Remarkably they differ from their stationary counterparts only in the explicit dependence of the drift, the memory kernel and the fluctuating force on one additional time. Apart from this the structure is the same as in equilibrium.

.5 Consequences

.5.1 An FDT-like equation

We now derive a relation between the auto-correlation of the fluctuating force and the memory kernel. The ”$\eta$-term” defined in eqn. (25) and the memory kernel as written in eqn. (21) are related via

[TABLE]

Since $\left\langle A^{*}_{\tau}\eta_{\tau}(\tau,t)\right\rangle=0$ (see eqn.(28)), one can write $\left\langle\eta_{\tau}(\tau,t)i\mathcal{L}A_{\tau}^{*}\right\rangle=\left\langle\eta_{\tau}(\tau,t)e^{i\mathcal{L}\tau}(1-P^{b}_{\tau})i\mathcal{L}A_{0}^{*}\right\rangle$. Then, by noticing that $G(\tau,\tau)=1$, one has $\eta_{\tau}(\tau,\tau)=e^{i\mathcal{L}\tau}(1-P^{b}_{\tau})i\mathcal{L}A_{0}$, which finally yields

[TABLE]

This equation is analogous to the fluctuation-dissipation theorem, but it holds for non-stationary processes. (Other fluctuation-dissipation-like theorems have been derived for non-stationary processes, but in a Fokker-Planck picture calabrese:2005 ; verley:2011 ). Note, however, that the analogy is just in terms of mathematical structure and not in terms of interpretation.

We can go one step further and try to find a similar relation for the correlation function $\left\langle\eta^{*}_{\tau}(\tau,t)\eta_{\tau}(\tau,t^{\prime})\right\rangle$. To do this, we write

[TABLE]

where

[TABLE]

Then, to compute $\bar{\eta_{\tau}}^{(n)}$, we first show from the definition of $\eta_{\tau}(\tau,t)$, using $\partial_{\tau}G_{\tau,t}=-i\mathcal{L}(1-P_{\tau}^{b})G_{\tau,t}$, the identity

[TABLE]

This relation will become useful later, but we first need to prove the following relation :

[TABLE]

First, the case $n=0$ is true and consists in eqn.(31). Let us assume that eqn.(35) is true for a certain $n$. Thus, we obtain

[TABLE]

Since $\bar{\eta_{\tau}}^{(n)}=(1-P_{\tau})\bar{\eta_{\tau}}^{(n)}$, and $\partial_{\tau}\eta_{\tau}(\tau,t)$ is proportional to $A_{\tau}$ (see eqn.(34)), we have $\left\langle\bar{\eta_{\tau}}^{(n)*}\partial_{\tau}\eta_{\tau}(\tau,t)\right\rangle=0$. We now use

[TABLE]

that we can rewrite as

[TABLE]

where we used eqn.(34). By inserting this last equation into (.5.1), and using again $\left\langle A_{\tau}^{*}\eta_{\tau}(\tau,t)\right\rangle=0$, we prove that eqn.(35) is also true at order $n+1$. thus, one can finally write

[TABLE]

The meaning of the ”noise”-term $\eta_{s}(s,t)$ is still not fully clear at this point. In particular, $\eta_{s}(s,t)$ is perpendicular to $A$ only with respect to $P^{b}_{s}$, i.e. at the initial time $s$. It might thus be interesting to look at its correlation with $A_{t}$. Let us take the derivative of order $n$ with respect to $t$ of equation (24) for $A$, multiply it by $A^{*}_{t}$, and take the trajectory average of the result. Then many terms cancel each other, such that one finally finds the following identity

[TABLE]

which is true for all $n$. The time derivative of this relation for the case $n=0$ yields

[TABLE]

In the limit of slow processes, we can neglect the last two terms to find

[TABLE]

which is always negative.

.5.2 Time-evolution of the memory kernel

The main difference between the equations of motion that we derived here and the well-known GLE for the stationary case lies in the explicit dependence of the Kernel on two times. In order to discuss this dependence further, we define the objects $\tilde{K}_{n}(t,\tau)$ by

[TABLE]

such that $K(t,\tau)=\tilde{K}_{1}(t,\tau)$. Because of the property $\partial_{t}G_{\tau,t}=G_{\tau,t}i\mathcal{L}(1-P_{t})$, we have

[TABLE]

By applying successive time-derivatives to $K(t,\tau)$ we obtain

[TABLE]

Eqn. (45) is a linear differential equation of infinite order for which the coefficients are well-controlled functions, meaning that one knows a priori their global properties in most of the situations. For instance, in most physical many-particle processes, $\omega_{n}(t)$ are bounded and infinitely differentiable functions, which ensures the solution of eqn. (45) to be a smooth differentiable function. We do thus not expect to see a discontinuous evolution of the memory kernel in such systems.

.5.3 A Taylor expansion of the memory kernel

As shown in the previous paragraphs, the non-stationary memory kernel $K(t,\tau)$ depends explicitly on the times $t$ and $\tau$, while in the stationary case it depends depend only on $t-\tau$. The kernel that we have derived recovers the $t-\tau$ behavior in the stationary limit. To show this, we come back to the definition of $K(t,\tau)$ from eqn. (.3), i.e. $K(t,\tau)A_{\tau}=e^{i\mathcal{L}\tau}P_{\tau}^{b}\left[i\mathcal{L}-\dot{P}_{\tau}^{b}\right]\left[1-P_{\tau}^{b}\right]G_{\tau,t}i\mathcal{L}A_{0}$. In the stationary limit, $P_{\tau}^{b}$ is a constant operator (i.e. it does not depend on $\tau$), such that $P_{\tau}^{b}=P^{b}$, $\dot{P}_{\tau}^{b}=0$ and

[TABLE]

Thus, we obtain

[TABLE]

with $k_{n}^{\infty}=\left\langle A_{0}^{*}P^{b}\left[i\mathcal{L}(1-P^{b})\right]^{n+1}i\mathcal{L}A_{0}\right\rangle\left\langle|A_{0}|^{2}\right\rangle^{-1}$. This limit proves that the well-known dependence of $K$ in $t-\tau$ is recovered as long as the projector $P_{\tau}^{b}$ is constant.

As pointed out above, in the stationary case the dependence of $K$ on $t-\tau$ only, allows to relate the stationary auto-correlation function $C(s,t)=C(0,t-s)$ and $K(t)$ in the Laplace (or Fourier) space, by making use of the convolution theorem in eqn. (29). In the non-stationary case, this is no longer possible, thus we need to find another way to evaluate the memory kernel. To do this, and since the integration in eqn. (29) runs over $\tau$, we perform a one-dimensional Taylor expansion of $K(t,\tau)$ at fixed $t$ and in the direction of $\tau$, around the point $\tau=t$, i.e.

[TABLE]

where $\kappa_{n}(t)\equiv\partial_{\tau}^{n}K|_{t=\tau}$. These coefficients can be directly computed from the formal definitions (19) and (.4) of $K(t,\tau)$ and $k_{n}(\tau,t_{1},\cdots,t_{n})$, in which the projection operators are applied only to objects of the form $i\mathcal{L}^{n}A_{0}$, with $n\in\mathbb{N}$. Therefore, $\kappa_{n}(t)$ can be expressed only in terms of the functions $\omega_{p}(t)$ defined in eqn. (18), with $p\leq n+2$. Let us describe here an example of computation, e.g. for $\kappa_{1}(t)$. First, we show from eqn. (19) that $\kappa_{1}(t)=\partial_{\tau}k_{0}(t)-k_{1}(t,t)$. Then, from eqn. (.4) we get $k_{0}(\tau)=\omega_{2}(\tau)-\omega_{1}^{2}(\tau)-\dot{\omega_{1}}(\tau)$, which yields $\partial_{\tau}k_{0}(t)=\dot{\omega_{2}}(t)-2\dot{\omega_{1}}\omega_{1}(t)-\ddot{\omega_{1}}(t)$, and $k_{1}(\tau,t_{1})=\omega_{3}(\tau)-\omega_{2}(\tau)[\omega_{1}(t_{1})+\omega_{1}(\tau)]+\omega_{1}(t_{1})[\dot{\omega}_{1}(\tau)-\omega_{1}(\tau)^{2}]-\dot{\omega_{2}}(\tau)$, which gives $k_{1}(t,t)=\omega_{3}(t)-2\omega_{2}(t)\omega_{1}(t)+\omega_{1}(t)\dot{\omega}_{1}(t)-\omega_{1}(t)^{3}-\dot{\omega_{2}}(t)$. We finally obtain $\kappa_{1}(t)=2\dot{\omega_{2}}(t)-3\dot{\omega_{1}}(t)\omega_{1}(t)-\ddot{\omega_{1}}(t)-\omega_{3}(t)+2\omega_{2}(t)\omega_{1}(t)-\omega_{1}(t)^{3}$. Such a procedure can be applied for any order, with increasing complexity. The first orders are reported in the appendix, in which we also show that the number of terms involved in $\kappa_{n}(t)$ grows exponentially with $n$.

We have thus derived a relation between the Taylor coefficients of the memory kernel and the dynamics of the coarse-grained variable as it can be obtained in a MD simulation. In addition, we have shown that despite the complexity of eqn. 19 the functional form of the kernel can be constructed without the need to compute an infinite number of nested integrals.

From this formalism, one can naturally define a timescale $\mathcal{T}$ associated to the time extent of the memory kernel, which may change as the process evolves. In fact,

[TABLE]

provides useful information about the timescale on which the memory kernel is relevant. This quantity can be easily sampled in MD simulations and then used to test for instance a Markovian approximation on the coarse-grained scale. Second, if one has a theoretical guess for the functional form of the memory kernel, one can test it by constructing the leading Taylor coefficients. Third, and most important, an accurate sampling of $\omega_{n}$, actually allows to construct the entire “generalized GLE”. We now show this in a numerical example.

A numerical example

To illustrate the use of the method, we carried out MD simulations of a two-dimensional model system, defined by one heavy particle of mass $M$ that interacts with bath particles, each of mass $m$, via a potential $V(r)=V_{0}\exp\left(-r/r_{0}\right)$, where $V_{0}$ and $r_{0}$ define the units of energy and distance, respectively. The bath particles do not interact with each other and we set $M=10^{3}m$. The averaged quantity for which we construct an equation of motion is the x-component of the momentum $p$ of the heavy particle, i.e.  $C(t)=\left\langle p_{x}(t)p_{x}(0)\right\rangle$.

We initialize the system by placing the heavy particle in the center of box with a velocity drawn from a Gaussian distribution associated to a certain temperature $k_{B}T=10^{-2}V_{0}$. The particles of the bath are initially distributed homogeneously in the box, except in a circular region of radius $R=30r_{0}$ around the central heavy particle. Their velocities are also picked from a Gaussian distribution associated to the same temperature $T$. The boundary conditions are reflective. The system is initially strongly out-of-equilibrium, and it reaches an equilibrium state after going through a transient phase. Thus, it is a well-suited test case for our method.

The computation of the local correlations $\left\langle A^{*}(t)A^{(n)}(t)\right\rangle$ is the central operation to perform in order to reconstruct the kernel. As these functions will be obtained by simulation, they can in practice be noisy. To increase the numerical accuracy, we use the following relation (valid for real variables)

[TABLE]

where $\alpha_{n,0}=1/2$, $\alpha_{2n,n}=(-1)^{n}$ for all $n$, and the remaining elements are determined by

[TABLE]

Sampling the functions $\gamma_{p}(t)=\left\langle{A^{(p)}_{t}}^{2}\right\rangle$ yields much weaker fluctuations than sampling the functions $\left\langle A_{t}A^{(n)}_{t}\right\rangle$. Note that if we take $A(t)=B^{(m)}(t)$, we obtain

[TABLE]

In this example, we sampled $\gamma_{n}(t)=\left\langle\left|d^{n}p_{x}/dt^{n}\right|^{2}(t)\right\rangle$ in order to compute the functions $\omega_{p}(t)$ and from those the Taylor coefficients $\kappa_{m}(t)$. We used the approximations $\kappa_{2n}(t)=\omega_{2n+2}(t)$ and $\kappa_{2n+1}(t)=0$, which turn out to be very good in this case. In figure (2) we show as an example $\gamma_{4}(t)$ and the function we used to interpolate it. The functional form of $\gamma$ was similar for all orders that we computed (until order 10). We also plot the Taylor expansion of the memory kernel constructed from these measurements, until order 18, as function of $t-\tau$ at various times $t$. We then use use it to solve the equation of motion.

Fig. 3 shows the momentum auto-correlation obtained directly from the MD simulation (solid line), i.e. the data extracted from the full “microscopic dynamics”. The dashed line is the best approximation that one can get if making the assumption that the dynamics of the averaged observable is Markovian. This approximation is commonly used in coarse-graining procedures. It clearly fails here, as it does not capture the short time-plateau.

The dotted and dash-dotted lines have been obtained by constructing the memory kernel according to eqn. 48 up to order 18 and then solving the generalized Langevin equation. Both capture the initial plateau very accurately.

If one approximates the kernel until it changes sign and then simply truncates it (dotted line), the long time behaviour of $\left\langle p_{x}(t)p_{x}(0)\right\rangle$ is not reproduced convincingly, because the expansion diverges. This can be fixed by extrapolating a tail on $K$ for large values of $t-\tau$. From the equilibrium dynamics, it is well known that the VACF and its memory kernel exhibit a $t^{-d/2}$ long-time tail, where $d$ is the dimension pomeau:1975 ; corngold:1972 . This would suggest extrapolation by an algebraic tail proportional to $(t-\tau)^{-1}$. However, the transient VACF starting from the non-equilibrium configuration does not exhibit such a long-time tail. If we fit an exponential decay, as is indicated by the large $n$ behaviour of the observed $\kappa_{n}$, the coarse-grained description (dash-dotted line) accurately captures the features of the directly computed correlation function. We have thus succeeded in constructing a coarse-grained description of this non-equilibrium model system.

.6 Summary

We have introduced a time-dependent projection operator that is of practical use, if one whishes to study non-equilibrium trajectory averages of phase space variables. We showed that, in the case of non-stationary dynamics, the equation of motion for the trajectory averages, eqn. (24), resembles the Generalized Langevin Equation. The only difference is an explicit dependence on an additional time in the drift term, the memory kernel and the fluctuating force. For all practical cases of application, the memory kernel is a smooth function in the additional time. We also derived an equation of motion for the auto-correlation function of the observable, eqn. (29). Remarkably, as in the stationary Mori-Zwanzig case, this equation does not contain noise.

We also showed how to systematically construct the memory kernel of a non-stationary GLE using as input data from experiments or MD simulations of the underlying microscopic dynamics. We thus provide a general strategy to develop coarse-graining procedures in classical atomistic computer simulations. In particular, we Taylor-expand the kernel and express its coefficients in terms instantaneous correlation functions of the variable of interest with its consecutive time-derivatives. If one can accurately measure time-derivatives up to order $n$, one can Taylor-expand until order $2(n-1)$. This allows to infer how long the system keeps track of its history, and can thus be used to test approximations that are often made in simulations on the coarse-grained scale (as e.g. the assumption of Markovian dynamics). If those approximations fail, the method can be used to construct appropriate equations of motion.

I Acknowledgements

We thank T. Franosch, A. Kuhnhold, M. Dolgushev and G. Amati for useful discussions. This project has been financially supported by the National Research Fund Luxembourg (FNR) within the AFR-PhD programme. Computer simulations presented in this paper were carried out using the HPC facility of the University of Luxembourg.

Appendix

Link to microscopic dynamics and a semi-analytic example

Here, we show an example of how to link the functions $\omega_{n}(t)$ to the microscopic dynamics of a system in practice, and how to use them to reconstruct the memory kernel. We focus on the case of Hamiltonian dynamics in a microscopic many-particle system, i.e.

[TABLE]

where $q$, $p$ denote position and momentum, respectively. To compute $\omega_{n}(t)=\left\langle A_{t}A^{(n)}_{t}\right\rangle/\left\langle A_{t}^{2}\right\rangle$, the obvious first step is to calculate $A^{(n)}=(i\mathcal{L})^{n}A$, with

[TABLE]

with $V_{\alpha,\alpha^{\prime}}=V(q_{\alpha},q_{\alpha^{\prime}})$. Of course, such a calculation becomes very lengthy and untractable for a non-specified generic variable. Therefore, we will do it on a relevant example, to show how such a study can be carried out in a practical case.

We choose a variable $A$ which depends only on the positions and is of the form

[TABLE]

For this type of variables, one can show

[TABLE]

where $\mathcal{V}_{\alpha_{1},\cdots,\alpha_{n}}$ is a object involving various derivatives of the potentials with respect to the positions $q_{\alpha_{1},\cdots,\alpha_{n}}$, as well as the momenta $p_{\alpha_{1},\cdots,\alpha_{n}}$ and derivatives of $f(q_{\alpha_{1}})$. Thus, we obtain

[TABLE]

where $f_{\gamma}=f(q_{\gamma})$. As a specific case, we now consider a system of identical particles $(m_{\alpha}=m_{\alpha^{\prime}}=m,\left\langle p_{\alpha}^{n}\right\rangle=\left\langle p_{\alpha^{\prime}}^{n}\right\rangle=\left\langle p^{n}\right\rangle)$, and we compute the density fluctuations, i.e.

[TABLE]

We rewrite now equation (Link to microscopic dynamics and a semi-analytic example) as

[TABLE]

In a high-temperature regime (or for weakly interacting systems), the second term can be neglected with respect to the first one. Thus, we have

[TABLE]

Finally, we assume that the phase-space distribution remains symmetric with respect to momenta and is of the shape

[TABLE]

This approximation consists in assuming a slow relaxation of the system towards equilibrium. We obtain then

[TABLE]

The assumption (61) is valid only if the system evolves slowly. Thus, we assume that the derivatives of $\omega_{n}(t)$ involved in the Taylor coefficients of the memory kernel are also negligible, yielding $\kappa_{0}(t)=\omega_{2}(t)$, $\kappa_{2}(t)=\omega_{4}(t)-\omega_{2}^{2}(t)$, $\kappa_{4}(t)=\omega_{6}(t)-2\omega_{4}(t)\omega_{2}(t)+\omega_{2}^{3}(t)$, … Because of the scaling of $\omega_{2n}$ with $n$, we can write

[TABLE]

where we calculate the coefficients $f_{2n}$ analytically from our formalism for the first orders. As an example we have $\kappa_{2}=\omega_{4}-\omega_{2}^{2}=(2/3)\omega_{4}$, i.e. $f_{2}=2/3$. The Taylor expansion becomes then

[TABLE]

where $\mathcal{T}(t)^{2}=m\beta(t)/k^{2}$. This sum can be numerically computed until very large orders without effort. We show in figure (4) the resulting reconstructed kernel (until order 80), as well as the solution of equation (29) for $C(t)$ using the latter kernel with a constant $\mathcal{T}$. The agreement of the reconstructed correlation function with well-known result $C(t)=C_{0}\exp\left[-t^{2}k^{2}/2m\beta\right]$ at high temperature hansen:1990 ; wallace:2016 is perfect. Of course, this result holds only within the assumptions made for this specific case, but the calculation shows that it may be possible to retrieve useful information about the functions $\omega_{n}(t)$ from the microscopic dynamics, and hence to partially infer the friction kernel. One can attempt to apply this sort of method for other types of processes and variables.

Computation of the Taylor coefficients in practice

Here, we show how to compute $\kappa_{n}$ practically. We first define the quantity

[TABLE]

The index l stands here for the order of derivation and is such that $m_{0}+\sum_{i=p+1}^{n}m_{i}=l$. Now, we show that the derivative of this object with respect to $\tau$ obeys to the identity

[TABLE]

where $\delta$ is the Kronecker symbol. The way this relation is used to find the Taylor coefficients is quite easy then. The first thing is to write the kernel $K(t,\tau)$ as a sum of integrals $I$ :

[TABLE]

To obtain the Taylor coeffient of order m, we apply then recursively the identity (2) to the $m+1$ first terms of the sum, and then we keep only the terms with $p=0$ (corresponding to the limit $\tau\rightarrow t$). As an example, let us compute $\kappa_{1}$.

[TABLE]

Once one has the expression of $\kappa_{n}$ as a function of $k_{m}$, one can insert the expression of these function in terms of $\omega_{p}$ from eq. (.4). We show here as a an example $k_{2}(\tau,t_{1},t_{2})$ :

[TABLE]

One finally can express the sum as a function of $\gamma_{q}$ by using eq. (50). Again, as an example, we have for $\omega_{2}(t)$ :

[TABLE]

We show here the example of $\kappa_{2}$ as a function of $\gamma_{q}$

[TABLE]

In general, one can express $\kappa_{n}(t)$ as

[TABLE]

The sum runs over all possible terms for which the dimension is the same on both sides of the equal sign. The coefficients $\alpha_{k}$ are calculated by the method presented in the previous lines. For each term of the sum, one must thus have $\sum_{i,j}p_{i,j}^{(k)}=n+2$ and $\sum_{i,j}(2i+j)p_{i,j}^{(k)}=n+2$. One can show that the number of terms in the sum grows roughly exponentially with $n$.

Action of $\dot{P}_{\tau}^{b}$

Here, we show how the operator $\dot{P}_{\tau}^{b}$ acts on an arbitrary dynamical variable $F$. One has

[TABLE]

which can be written out as

[TABLE]

where we have defined $B_{\tau}=F_{\tau}-A_{\tau}\left<A_{\tau}^{*}F_{\tau}\right>/\left<|A_{\tau}|^{2}\right>$. But the last term in $B_{\tau}$ is simply $e^{i\mathcal{L}\tau}P_{\tau}^{b}F_{0}$, and thus we obtain

[TABLE]