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

The concept of projection operators as a tool to reduce the dimension of physical systems goes back to the 1960s. Based on the work by Nakajima,¹ Zwanzig showed how projection operators can be used to derive an equation of motion for almost arbitrary phase space observables.^2,3 Five years later, Mori derived an analogous result⁴ employing a different kind of projection operator. In more recent works of Chorin et al.,^5,6 the topic was 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.^7,8 However, Grabert’s approach was very general and did not include any explicit generalization of either Zwanzig’s or Mori’s projection operator to non-equilibrium dynamics. (Notably, a derivation of a non-stationary equation of motion using a time-dependent generalization of Mori’s projection operator formalism was presented by Nordholm in his Ph.D. thesis as early as 1972,⁹ but unfortunately, this was not taken up by the community as it was not published elsewhere.)

Concrete applications of time-dependent projection operators to model colloidal suspensions and undercooled liquids were provided by Shea, Oppenheim, and Latz in the 1990s.^10–13 In recent work, Meyer et al.^14,15 and te Vrugt et al.^16,17 used time-dependent projection operators similar to those of Mori to derive a general equation of motion for coarse-grained variables in non-equilibrium systems, including even systems under time-dependent external driving. (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.¹⁸) Furthermore, Meyer et al. introduced a general and fast method to compute the memory kernels appearing in their non-stationary Generalized Langevin Equation (abbreviated as “nsGLE” and always referring to the version of Meyer et al.) from a set of time-resolved values of an observable for individual trajectories.¹⁹ As such data are often accessible through molecular dynamics simulations, the question arises naturally if one can apply the same formalism 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.

Español and Vázquez have used a projection operator formalism of the Zwanzig-type in order to coarse-grain dynamics that are governed by the Fokker–Planck equation.²⁰ Importantly, they observed that under conditions of time-scale separation, the coarse-grained dynamics are again governed by a Fokker–Planck equation. A similar but more explicit route has been taken in the work of Kranz et al.²¹ In both cases, the average over the stochastic degrees of freedom was taken before the projector was applied. Here, we take a different route in order to obtain an equation of motion that treats as independent the trajectories that do not differ in their initial physical configuration but only in the explicit stochastic contribution. To illustrate that the framework can in fact be extended to stochastic processes in such a way, we will start with a brief recapitulation of the derivation of the nsGLE in Sec. II before discussing the case of stochastic dynamics in Sec. III.

Given an initial phase space distribution ρ(Γ, 0), which is neither necessarily stationary nor, in particular, the equilibrium distribution, and a (time-dependent) Liouvillian $Lt$⁠, 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-dependent 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, Γ′ ≔ (Γ, τ).¹⁵ 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 for now.

One can write the time evolution of an observable $A(Γ)$⁠, which is a function of the augmented phase space coordinates, as

Here, Γ₀ 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 in 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 Γ₀ from now on.

By taking the time derivative of Eq. (1) and using Grabert’s approach^7,8 of applying time-dependent projection operators, one obtains

Here, $Pt$ is a time-dependent projection operator, $Qt≔1−Pt$ 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 δ(τ) 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 η_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 parentheses, 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,¹⁴ namely,

The derivation of the nsGLE presented in Refs. 14 and 15 requires taking time-derivatives of the observable. The derivation can therefore not be directly applied to stochastic dynamics. However, the idea of the augmented phase space can be used to tackle this problem.

We construct an augmented phase space 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. Furthermore, we note that the Liouvillian does not affect the random variables (⁠$iLRi=0∀Ri∈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 “pseudorandom 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 have 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 these 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₁, …, 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 a 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$⁠. [Note that we can assume that the Liouvillian couples some physical degree of freedom to the random variables, as given in Eq. (10), because the augmentation of the phase space includes a time coordinate.] 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 whether 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 whether 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 whether 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 additional consideration.

The equations of motion for the two cases, i.e., dynamics with the pseudorandom numbers generating the “noise” and dynamics with the analytic term Eq. (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.²² 

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. 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 straightforward 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 the one with s_n(t) data, deviate less than that after a single integration time step. Furthermore, 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)=Δq1(t),…,ΔqN(t),Δp1(t),…,ΔpN(t)T$ with the generalized coordinates q_i(t) and momenta p_i(t)]. Furthermore, with a meaningful norm $‖ΔΓ(t)‖$⁠, e.g., $‖ΔΓ(t)‖=a1Δq1(t)2+⋯+a2NΔpN(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 initial separation is small (ΔΓ → dΓ) and that the flow field $Γ̇(Γ)$ is bounded. [Note that $Γ̇(Γ)$ 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 the phase space. By demanding $Γ̇(Γ)$ to be analytic, it must also be continuous, and hence, $Γ̇(Γ)$ is bounded in this subspace, which suffices for the following considerations.] 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. Furthermore, 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, δ(τ) is the contribution syncing the augmented phase space coordinate τ with the observation time t (cf. Ref. 15), and p(R) is the joint probability distribution of all degrees of freedom describing the random noise. Note that in the more general case where the probability distribution of the initial physical degrees of freedom and the stochastic ones does not separate, every initial physical configuration can have its own probability distribution of the stochastic degrees of freedom, allowing for stochastic contributions that differ for different trajectories.

The above-mentioned arguments 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 $RN$⁠, the new space will be of the type $RN×R×FN$⁠, where F is the space of power series converging in the whole interval of observation. Here, $RN$ contains, again, the physical degrees of freedom and the additional $R$ contains, again, the degree of freedom for the trajectory time. Furthermore, 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 the points in this space by (Γ′) = (Γ, τ, f). The Liouvillian in this space can then be written as

where $L$ is the Liouvillian of the original system, acting and depending only on (Γ), and q_i are the physical degrees of freedom. Hence, the action on the physical degrees of freedom is given by

whereas the f_i(x) remain unchanged under the action of the Liouvillian (⁠$iL′fi(x)≡0$⁠).

If we consider an observable $A′(Γ′)$⁠, which depends only on the physical degrees of freedom (Γ) for which $iL′qi=dqi/dt$⁠, 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 $RM$⁠, we could introduce a mapping $RM→FN$⁠, define the probability density in these coordinates, and use a simple measure of the $RN+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. Hence, every kick is determined by the time of occurrence (center of the 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.

Note that the arguments presented here apply to coarse-graining of explicitly time-dependent Liouvillians and to any kind of stochastic process that is bounded and either time-discrete or time-continuous, each realization of which can be approximated by an analytic function. This constitutes a generalization compared to previous work,²⁰ which applies to the case where the Liouvillian is not explicitly time-dependent and the stochastic microscopic process can be described by a Fokker–Planck equation.

We have shown that the general framework of the non-stationary 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 a thermostat or barostat. The generalization was achieved by further abstraction from the usual phase space toward a more convenient space that includes degrees of freedom capturing the stochastic contributions.

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

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