Comments on the validity of the non-stationary Generalized Langevin Equation as a coarse-grained evolution equation for microscopic stochastic dynamics

We recently showed that the dynamics of coarse-grained observables in systems out of thermal equilibrium are governed by the non-stationary generalized Langevin equation [J. Chem. Phys. 147, 214110 (2017), J. Chem. Phys. 150, 174118 (2019)]. The derivation we presented in these two articles was based on the assumption that the dynamics of the microscopic degrees of freedom was deterministic. Here we extend the discussion to stochastic microscopic dynamics. The fact that the non-stationary Generalized Langevin Equation also holds for stochastic processes implies that methods designed to estimate the memory kernel, drift term and fluctuating force term of this equation as well as methods designed to propagate it numerically, can be applied to data obtained in molecular dynamics simulations that employ a stochastic thermostat or barostat.


I. INTRODUCTION
The concept of projection operators as a tool to reduce the dimension of physical systems goes back to the 1960s. Zwanzig showed how projection operators can be used to derive an equation of motion for almost arbitrary phase space observables [1,2]. Five years later, Mori derived an analogous result [3] employing a different kind of projection operator. In more recent works of e.g. Chorin et.al. [4,5], the topic is addressed in a more general mathematical setting abstracting from the context of physical observables and phase space.
Grabert showed in the 1970s how the projection operator technique can be extended to non-stationary processes by using time-dependent projection operators [6,7]. However, Graberts approach was very general and did not directly relate to either Zwanzig's or Mori's projection operator. Concrete applications of time-dependent projection operators to model colloidal suspensions and undercooled liquids were provided by Shea, Oppenheim and Latz in the 1990s [8][9][10][11]. In recent work Meyer et.al. [12,13] as well as Vrugt et al. [14,15] used timedependent projection operators similar to Mori's to derive a general equation of motion for corase-grained variables in non-equilibrium systems, including even systems under time-dependent external driving. They obtained a non-stationary version of the Generalized Langevin Equation (abbreviated as "nsGLE"). (A similar attempt was made by Kawai and Komatsuzaki via the Zwanzig projection operator, but it turned out to be more involved mathematically than the Mori approach [16].) Further, Meyer et.al. introduced a general and fast method to compute the memory kernels appearing in the ns-GLE from a set of time-resolved values of an observable for individual trajectories [17]. As such data is often accessible through molecular dynamics simulations, the question arises naturally if one can apply the same for-malism in the context of stochastic microscopic propagators (e.g. dynamics generated using thermostats and barostats). This is not directly clear because the derivation of the nsGLE demands a deterministic Liouvillian. To illustrate that the method can in fact be extended to stochastic processes, we will start with a brief recapitulation of the derivation of the nsGLE in section II before discussing the case of stochastic dynamics in section III.

II. RECAPITULATION AND NOTATION
Given an initial phase space distribution ρ(Γ, 0), that is neither necessarily stationary nor, in particular, the equilibrium distribution, and a (time-dependent) Liouvillian L t , the phase space distribution is determined for any future time by the Liouville equation. Here, Γ denotes the collective set of phase space coordinates. If one intends to apply a time-dependent projection operator to a time-depentend Liouvillian, it can be useful to switch from the usual phase space representation to an "augmented" phase space that includes one additional coordinate for the system time Γ := (Γ, τ ) [13]. Then a Liouvillian L in the augmented phase space can be defined generating the original dynamics without an explicit time-dependence (L does not depend on t) and projections can be carried out on this Liouvillian. Unless stated otherwise, the following calculations in this section are carried out in the augmented phase space and the prime in the notation is dropped. Now, one can write the time evolution of an observable A(Γ), that is a function of the augmented phase space coordinates, as Here, Γ 0 denotes the initial point in the augmented phase space of one trajectory, Γ t is the point in the augmented phase space reached by the same trajectory at time t, and A t is the value of the observable for a specific trajectory as a function solely of time. We implicitly assumed that the dynamics can be described by analytic functions on arXiv:2103.07112v1 [cond-mat.stat-mech] 12 Mar 2021 the whole observation interval. To allow for easy readability, we will omit spelling out the dependencies on the phase space coordinates of A(Γ) and the insertion of the initial point in phase space Γ 0 from now on. By taking the time derivative of eq. (1) and using Graberts approach [6,7] of applying time-dependent projection operators one obtains Here, P t is a time-dependent projection operator, Q t := 1 − P t is its orthogonal complement and where exp − denotes a negatively time-ordered exponential function. Next, we introduce a time-dependent product on the set of phase space observables by where ρ(Γ, 0) is the initial probability density in the augmented phase space obtained by multiplying the probability density of the initial ensemble with the term δ(τ −0) syncing the observation time and the augmented time coordinate. One can define a specific projection operator by where A is a specific observable of interest and X is any phase space function. Inserting this explicit projection operator, eq. (2b) can be rewritten as Here, the quantities ω(t), K(τ, t) and η 0,t are defined by Note that the functions ω(t) and K(t, τ ) are the same for different trajectories, if their initial conditions are drawn from one given initial probability density (i.e. for trajectories drawn out of one given non-equilibrium ensemble or "swarm" of trajectories). Hence, we denote their time dependencies using parenthesis, whereas quantities with time as a subscript, such as A t and η t , do depend on the individual trajectory. K(t, τ ) is the so called memory kernel describing non-local (in time) contributions for the equation of motion of A t . It can be shown that these quantities fulfill a relation that is similar in structure to a fluctuation-dissipation theorem [12], namely

III. INCLUDING STOCHASTICITY
The derivation of the nsGLE presented in Refs. [12,13] requires taking time-derivatives of the Liouvillian. The derivation can therefore not be directly applied to stochastic dynamics. However, the idea of the augmented phase space can be used to to tackle this problem.
We construct an augmented phase such that we can keep track of the random numbers generated along a single trajectory, i.e. of the values taken by the noise function at each time step. Here, we assume that there is only a finite number of random numbers per trajectory, which is the case for computer simulations, in which time is discretized. If we refer to the set of random numbers per trajectory by R, then the augmented phase space is given by (Γ ) := (Γ, τ, R), where we have already included the time-coordinate τ for convenience. Further, we note that the Liouvillian does not affect the random variables (iLR i = 0 ∀R i ∈ R) and, thus, they may be regarded as integrals of motion. The values of the random variable taken along a single trajectory, when interpreted as a function of simulation time, may look like the bold red line in fig. 1. We labeled the line by pRN for "pseudo random number" as this is the case in a typical computer simulation, but the arguments hold equally well for numbers generated by a truly stochastic process.
Note that the pRN data has discontinuous jumps from one time step to another and, hence, naive coupling of physical degrees of freedom to the bold red pRN curve would introduce non-analytic behavior. However, we argue in the following that one can interpret this data as the limit of a series of analytic functions. Let us take a trajectory with M time steps in total and random numbers given by R 1 , · · · , R M . To keep the notation simple, we will assume that ∆t = 1 and that the jumps of the pRN data occur in the middle between consecutive time steps. Then, the value of the pRN as function of time can be expressed through Here, θ(t) denotes the Heaviside step function. In contrast to pRN(t), eq.9, the function defined through with n ∈ N is analytic and convergent on R 1 . Figure 1 shows two exemplary curves of s n (t) for two different values of n. s n (t) converges pointwise to pRN(t) in the limit n → ∞ if one neglects the null set of points where the jumps for the pRN(t) occur. Molecular dynamics simulations are aimed at a numerical integration and, hence, one needs to check that also the limit of the integrals over eq. (10) converges to the integral over eq. (9). Here, pointwise convergence is not a sufficient condition. However, if the values of the pRN are bounded, s n (t) is uniformly integrable ∀n ∈ N and the convergence of the integrals is given as well.
Finally, we need to check that also the value of the observable along the trajectory produced using s n (t) converges to the one obtained by using pRN(t). As our observable is an analytic function of the phase space coordinates, it suffices to check that the distance in phase space between trajectories produced using s n (t) and pRN(t) vanishes in the limit n → ∞. However, this point needs some additonal consideration.
The equations of motion for the two cases, i.e. dynamics with the pseudo random numbers generating the "noise" and dynamics with the analytic term 10, will generate different trajectories. In chaotic systems the distance between the trajectories is expected to grow in time. More precisely, if one defines a metric of phase space to quantify the deviation of two states, Oseledets's theorem states that this deviation grows (or shrinks) asymptotically exponentially (described by its Lyapunov exponents) for non-integrable systems [18]. Now assume that one intends to analyze molecular dynamics trajectories of some finite duration. The deviation between two final states after the complete simulation-time due to some "small perturbation" at the beginning will always remain finite -even if it is usually very large -and hence it can be scaled down by weakening the initial perturbation, which in turn can be achieved by increasing n. Using the Lyapunov exponents of the given system and demanding that the trajectories deviate less than some desired value (e.g. a value on the order of numerical precision), it is straight forward to find some threshold for the allowed deviation after a single integration time-step. However small this threshold may be, naturally one can choose a value of n large enough such that the two dynamics, the one with pRN(t) and one with s n (t) data, deviate less than that after a single integration time-step. Further, assuming a Liouvillian generating analytic dynamics, not only the deviation in the long-time limit but also short-time fluctuations between the two dynamics can be reduced to an arbitrarily small but nonzero value.
Formally, this can be expressed in the following way: Let ∆Γ(t) be the distance between two initially close trajectories in the physical degrees of freedom of phase space (e.g. ∆Γ(t) = (∆q 1 (t), · · · , ∆q N (t), ∆p 1 (t), · · · , ∆p N (t)) T with the generalized coordinates q i (t) and momenta p i (t)). Further, with a meaningful norm ∆Γ(t) , e.g. ∆Γ(t) = a 1 ∆q 1 (t) 2 + · · · + a 2N ∆p N (t) 2 where the a i have a strictly positive value and cancel the units of the corresponding factors, one gets ∆Γ(t) ∝ exp(λt) in the long time limit. Here, λ is the largest Lyapunov exponent. By demanding that at the end of the simulation time t = T the norm of the deviation takes some arbitrarily small but finite value ∆Γ(T ) , one can now calculate a threshold for the deviation after a single simulation-step, namely This value may become absurdly small but will always remain finite. Next, we assume that the flow fieldΓ(Γ) is bounded 2 and that the initial separation is small (∆Γ → 2Γ (Γ) as a function of the complete phase space will in most cases be unbounded. However, any finite subspace of finite trajectories will lie in some compact subset of phase space. By demandinġ Γ(Γ) to be analytic it must also be continuous and henceΓ(Γ) is bounded in this subspace which suffices for the following considerations.
dΓ). Then, we can write with some constant c. Hence, it is clear that also the short time fluctuations become arbitrarily small as the initial separation diminishes. Thus, by reinterpreting the noise as a function obtained as a limit of analytic functions, one can analyze data generated with stochastic propagators by applying the methods provided by the nsGLE. Note that we made no assumption on how the random variables R for each trajectory are obtained. Hence, the above reasoning holds for both pseudorandom numbers and truly random numbers. Further, we note that the probability distribution of the random numbers enters the initial phase space probability distribution ρ (Γ , 0). In the simplest case when the initial phase space probability distribution of the physical degrees of freedom and the random degrees of freedom are independent, which is by no means necessary nor demanded by the formalism of the nsGLE, one could write where ρ(Γ, 0) is the initial phase space density in the original (i.e. not augmented) phase space, δ(τ − 0) is the contribution syncing the augmented phase space coordinate τ with the observation time t (cf. [13]), and p(R) is the joint probability distribution of all degrees of freedom describing the random noise. The arguments above handled the special case where the noise is a step function. This case is quite common in the context of computer simulations, but we will generalize the approach in the following to allow for more general types of stochastic processes. First, we introduce a new space, which will replace the augmented phase space. Given that the initial phase space of our system is of type R N , the new space will be of the type R N ×R×F N , where F is the space of power series converging in the whole interval of observation. Here, R N contains again the physical degrees of freedom and the additional R contains, again, the degree of freedom for the trajectory time. Further, for every physical degree of freedom q i there is now an analytic function f i (x) ∈ F that will describe the realization of the stochastic properties. We will denote points in this space by (Γ ) = (Γ, τ, f ). The action of the Liouvillian in this space, L , on the physical degrees of freedom q i can be expressed as Here, L is the Liouvillian of the original system, acting and depending only on (Γ). Again, the f i (x) can be interpreted as integrals of motion (expressed though iL f i (x) ≡ 0).
If we consider an observable A (Γ ) which depends only on the physical degrees of freedom (Γ), we can write From this, we can derive the nsGLE analogously as before. However, for the application of the projection operators we need the probability density ρ (Γ , 0) and we need to carry out integrals over all degrees of freedom of the augmented phase space. To circumvent the problem of needing a measure for this infinite dimensional space, we point out that if the analytic functions can be specified by a finite number of real control parameters R M , we could introduce a mapping R M → F N , define the probability density in these coordinates, and use a simple measure of the R N +1+M , e.g. the Lebesgue measure.
An example for such a mapping could be a particle moving in three dimensions and getting random kicks. Assume that these kicks can be described by a force that is a Gaussian in time. So, every kick is determined by the time of occurrence (center of Gaussian), the width of the Gaussian, an amplitude, and a direction (e.g. specified by two real numbers). Hence, every individual kick can be specified by five real numbers.

IV. CONCLUSION
We have shown that the general framework of the nonstationary Generalized Langevin Equation can be applied to a wide range of processes with (pseudo)stochastic contributions. This includes particularly simulation data obtained using some kind of stochastic propagator, such as e.g. a thermostat or barostat. The generalization was achieved by further abstraction from the usual phase space towards a more convenient space that includes degrees of freedom capturing the stochastic contributions.

V. ACKNOWLEDGMENTS
The authors acknowledge funding by the Deutsche Forschungsgemeinschaft (DFG, German Research Foundation)-Project No. 430195928, and useful mathematical remarks from Fabian Coupette.

VI. DATA AVAILABILITY
Data sharing is not applicable to this article as no new data were created or analyzed in this study.