The interplay between memory and potentials of mean force: A discussion on the structure of equations of motion for coarse-grained observables

When a coarse-grained model for a polymeric system or a biological macromolecule is designed, groups of atoms are merged into larger units. Then approximate equations of motion for these units are solved to predict the evolution of the model [1–3]. Depending on context, the coarse-grained degrees of freedom can range from the positions of small molecules or chemical groups to reaction coordinates such as the relative orientation of structural motifs in a biomolecule. In principle, to obtain the exact equation of motion of a coarse-grained degree of freedom, one would need to integrate out systematically the atomistic degrees of freedom. However, as this is a very difficult task, researchers usually resort to effective models. For instance, the underdamped non-linear Langevin equation is frequently used

$$\begin{equation} m \ddot{x}_t = -\frac{\mathrm{d}W(x)}{\mathrm{d}x} \bigg\vert_{x=x_t} - \gamma \dot{x}_t + \sqrt{2\gamma k_{\rm B} T}\xi_t, \end{equation} \tag{ 1 }$$

where x_t is the position of a coarse-grained unit at time t (respectively, the value of a more general reaction coordinate), m is a generalized mass, γ is a friction coefficient, W(x) is an effective potential, $k_{\rm B} T$ is the thermal energy and $\xi_t$ is white Gaussian noise [4].

The notation $\ddot{x}$ and $\dot{x}$ for the time derivatives is frequently used in the physics literature to indicate that eq. (1) could be interpreted in analogy to the Newtonian equation of motion of a particle in a potential energy landscape. However, as $\xi_t$ is a stochastic process, the terms $\ddot{x}$ and $\dot{x}$ are stochastic derivatives, and W(x) is not an external potential but the potential of mean force

$$\begin{equation*} W(x):= - k_{\rm B}T \ln{\left(\rho_X^{\text{eq}}(x)\right)}. \end{equation*}$$

Here $\rho_X^{\text{eq}}(x)$ is the so-called "relevant density" of the coarse-grained observable X, i.e., the probability of the observable X having the value x in the equilibrium ensemble [5]. Thus the analogy to Newtonian mechanics can be misleading. For canonical dynamics one also often encounters the terms effective free energy and free energy landscape for $- k_{\rm B}T \ln{(\rho_X^{\text{eq}}(x))}$, denoted by $\Delta G(x)$.

The dynamics of the atomistic degrees of freedom which have been integrated out, in general, produce memory effects. Therefore integro-differential equations are also often used to model coarse-grained variables, such as, e.g.,

$$\begin{equation} m \ddot{x}_t = - \left. \frac{\mathrm{d}W(x) }{\mathrm{d}x}\right\vert_{x=x_t} - \int_0^t \mathrm{d}\tau \; K(t-\tau) \dot{x}_\tau + f_t, \end{equation} \tag{ 2 }$$

where K is the memory kernel and f the flucutating force, which is related to K by the second fluctuation-dissipation theorem (see refs. [6–14] for examples of recent work in which this equation is used to model coarse-grained dynamics). In this letter, we type-set times in parentheses if the time dependence is on the level of the ensemble average (as, e.g., in the memory kernel $K(t-\tau)$) and times as subscripts if the dependence is on the level of the individual trajectory (as, e.g., in x_t ).

These effective equations of motion are frequently used in the soft matter modelling community. They provide a practical pathway to coarse-grained modelling, because the functions W(x) and K(t) can be parameterised and then fitted to simulation data. Therefore it is interesting to check under which assumptions these equations can be derived from first principles. We are aware of only two publications, in which derivations for eq. (2) are shown [15,16]. Note that we are referring specifically to the form of the generalized Langevin equation, in which the variable x enters W(x) non-linearly, while the memory term is linear in $\dot{x}$, and $K(t-\tau)$ and f_t are related by the second flucutation-dissipation theorem. Other forms of the generalized Langevin equation, as discussed, e.g., in refs. [4,17–22], are not in the focus of our letter. We also acknowledge that, if one replaces the potential of mean force W(x) in eq. (2) with an external potential $V_{\rm ext}(x)$, one obtains the equation of motion of a specific, well-known model system: one particle linearly coupled to a bath of harmonic oscillators [23]. Also this is not the type of problem we are referring to in this letter. We are aiming at coarse graining the dynamics of complex systems such as polymers and biomolecules, in and out of equilibrium. The question we address in this letter is if an equation of motion with the structure of eq. (2) can be derived for a broad class of systems and observables, as claimed in refs. [15] and [16].

A useful framework to tackle this task is the projection operator formalism as originally introduced by Zwanzig [24,25] and Mori [26]. Let $\Gamma = \{q_i,p_i\}$ denote the phase space coordinates of the microscopic system and $\mathrm{i}\mathcal{L}$ the Liouvillian, which for now shall not be explicitly time dependent (we will come to time-dependent Liouvillians later). $\vec{\mathbb{A}}(\Gamma)$ shall denote a set of phase space fields, e.g., the coarse-grained observables for which we intend to derive an equation of motion. We use blackboard bold for phase-space functions and italics for the value they take at specific points in phase space, i.e., $\vec{A}_t=\vec{\mathbb{A}}(\Gamma_t)$. In the case of Hamiltonian dynamics, the action of the Liouvillian on the fields $\mathbb{A}_i(\Gamma)$ is given by

$$\begin{align} \dot{\mathbb A}_i&=\mathrm{i}\mathcal{L}\mathbb{A}_i = \left\{\mathbb{A}_i,\mathbb{H}\right\}=\sum\limits_j\frac{\partial \mathbb{A}_i}{\partial q_j}\frac{\partial\mathbb{H}}{\partial p_j}-\frac{\partial \mathbb{A}_i}{\partial p_j}\frac{\partial\mathbb{H}}{\partial q_j}. \end{align} \tag{ 3 }$$

The equation of motion of the observables can be integrated formally

$$\begin{align} \frac{\mathrm{d}\vec{A}_t}{\mathrm{d}t} &= \exp(t\mathrm{i}\mathcal{L})\mathrm{i}\mathcal{L}\vec{\mathbb{A}}(\Gamma)\big|_{\Gamma=\Gamma_0}. \end{align} \tag{ 4 }$$

The right-hand side is the time-evolution operator for a time span of t acting on the time derivative of $\vec{\mathbb{A}}$. (Note that the initial phase-space coordinates $\Gamma_0$ are inserted after performing this operation.) We introduce projection operators $\mathcal{P}$ which act on the space of phase-space fields. Using the Dyson-Duhamel identity, eq. (4) can be written as

$$\begin{align} \frac{\mathrm{d}\vec{A}_t}{\mathrm{d}t} &= \int\limits_0^t\mathrm{d} s\,\exp(s\mathrm{i}\mathcal{L})\mathcal{P}\mathrm{i}\mathcal{L}\mathcal{Q}\exp\left((t-s)\mathrm{i}\mathcal{L}\mathcal{Q}\right)\mathrm{i}\mathcal{L}\vec{\mathbb{A}}\nonumber\\ &\phantom{=}+\exp(t\mathrm{i}\mathcal{L})\mathcal{P}\mathrm{i}\mathcal{L}\vec{\mathbb{A}}+\mathcal{Q}\exp\left(t\mathrm{i}\mathcal{L}\mathcal{Q}\right)\mathrm{i}\mathcal{L}\vec{\mathbb{A}}, \end{align} \tag{ 5 }$$

where $\mathcal{Q}:= 1-\mathcal{P}$. (We dropped the explicit insertion of the initial point in phase-space $\Gamma_0$ on the right-hand side.) Next, we need to choose a specific projector. Two of the most prominent types of projectors are the Zwanzig and the Mori projector. As we will see in the following, the Zwanzig projector has the useful property that the second term in eq. (5) turns into a derivative of a potential of mean force under certain conditions. On the other hand, the Mori projector, which is linear in the observable(s), yields a fluctuation-dissipation relation.

We begin with a projection operator similar to Zwanzig's original one to bring the second term of the right-hand side of eq. (5) into the form of a derivative of a potential of mean force. We define

$$\begin{align} \mathcal{P}\mathbb{X}(\Gamma) &= \frac{1}{\rho_A(\vec{\mathbb{A}}(\Gamma))}\int\textrm{d}\Gamma'\rho^{\text{eq}}(\Gamma')\mathbb{X}(\Gamma') \delta(\vec{\mathbb{A}}(\Gamma')-\vec{\mathbb{A}}(\Gamma)) \end{align} \tag{ 6 }$$

with

$$\begin{align} \rho_A(\vec{A}) &= \int\mathrm{d}\Gamma\,\rho^{\text{eq}}(\Gamma)\delta(\vec{\mathbb{A}}(\Gamma)-\vec{A}), \end{align} \tag{ 7 }$$

where $\mathbb{X}$ is an arbitrary phase-space field. The normalization factor $\rho_A(\vec{A})$ is the relevant density of the variables $\mathbb{A}_i$. In the next steps, we use the canonical equilibrium density $\rho^\text{eq}(\Gamma)\propto\exp(-\beta\mathbb{H}(\Gamma))$, where $\beta=1/k_{\rm B} T$, to project $\mathbb{X}$ onto a set of phase space fields $\{\mathbb{A}_i\}$. (This is not mandatory, a similar derivation can also be carried out for other ensembles.)

To derive an equation of a structure similar to eq. (2), we assume that the Hamiltonian can be split into a simple kinetic contribution and a potential that depends only on the generalized coordinates

$$\begin{align} \mathbb{H}(\Gamma) &= \sum\limits_{j} \frac{p_j^2}{2m_j} + \mathbb{V}(q_1,\ldots,q_N). \end{align} \tag{ 8 }$$

Now we project onto a set of two rather specific observables, $\vec{\mathbb{A}}=(\mathbb{A}_1,\mathbb{A}_2)^\top$. The first one takes the form $\mathbb{A}_1=\alpha_0+\sum_j\alpha_jq_j$ with constant prefactors $\alpha_i$. An example for such an observable would be a center of mass of a set of atoms. The second observable is the time derivative of the first one, namely $\mathbb{A}_2 = \mathrm{i}\mathcal{L}\mathbb{A}_1 = \sum_j \alpha_jp_j/m_j$.

To obtain an equation similar to eq. (2), where the left-hand side is the second time derivative of the observable, we consider the evolution of $\mathbb{A}_2$. Thus, the second term in the second component of eq. (5) reads

$$\begin{align} \exp(t\mathrm{i}\mathcal{L})\mathcal{P}\mathrm{i}\mathcal{L}\mathbb{A}_2 &= -\exp(t\mathrm{i}\mathcal{L})\frac{1}{\rho_A(\vec{\mathbb{A}}(\Gamma))}\int\mathrm{d}\Gamma'\,\rho^\text{eq}(\Gamma')\nonumber\\ &\quad\times\delta(\vec{\mathbb{A}}\left(\Gamma'\right)-\vec{\mathbb{A}}(\Gamma))\sum\limits_{j}\frac{\alpha_j}{m_j}\frac{\partial \mathbb{V}(\Gamma')}{\partial q_j'}. \end{align} \tag{ 9 }$$

To relate this expression to a potential of mean force, we take the derivative of the relevant density with respect to the first coordinate

$$\begin{align} \frac{\partial \rho_A\left(\vec{A}\right)}{\partial A_1}&= \int\mathrm{d}\Gamma'\,\rho^\text{eq}\left(\Gamma'\right)\delta\left(\mathbb{A}_2\left(\Gamma'\right)-A_2\right)\nonumber\\ &\quad\times\mu\sum\limits_j\frac{\alpha_j}{m_j}\frac{\partial}{\partial q_j'}\delta\left(\mathbb{A}_1\left(\Gamma'\right)-A_1\right), \end{align} \tag{ 10 }$$

where $\mu=\sum_i m_i/\alpha_i^2$. Next, we carry out an integration by parts and use the fact that

$$\begin{align} \frac{\partial }{\partial q_j}\rho^\text{eq}(\Gamma) &= -\beta \frac{\partial\mathbb{V}(\Gamma)}{\partial q_j}\rho^\text{eq}(\Gamma) \end{align} \tag{ 11 }$$

to obtain

$$\begin{align} \frac{\partial \rho_A(\vec{A})}{\partial A_1} &= \beta\mu\int\mathrm{d}\Gamma'\,\rho^\text{eq}\left(\Gamma'\right)\delta(\vec{\mathbb{A}}\left(\Gamma'\right)-\vec{A})\nonumber\\ &\quad\times \sum\limits_j \frac{\alpha_j}{m_j} \frac{\partial\mathbb{V}(\Gamma')}{\partial q_j'}. \end{align} \tag{ 12 }$$

Comparing eq. (9) and eq. (12), we see that

$$\begin{align} \exp(t\mathrm{i}\mathcal{L})\mathcal{P}\mathrm{i}\mathcal{L}\mathbb{A}_2 &= -\frac{1}{\mu}\frac{\partial}{\partial A_1}\underbrace{(-k_{\rm B} T\ln(\rho_A(\vec{A})))}_{=:W(\vec{A})}\Big|_{\vec{A}_t}, \end{align} \tag{ 13 }$$

where $W(\vec{A})$ is the potential of the mean force. Note that $\partial W(\vec{A})/\partial A_1$ does not depend on A₂, because the A₂-dependent factors in the numerator and denominator cancel each other.

A similar expression can be derived for a projection on multiple observables, e.g., for a projection on the centers of mass of several "blobs" (united atoms in a polymeric system). In this case one replaces the first observable $\mathbb{A}_1$ by a set of observables $\mathbb{A}_{1i}$, where each of these observables has the form $\mathbb{A}_{1i}=\alpha_{0i}+\sum_j \alpha_{ji}q_j$, with constant prefactors $\alpha_{ji}$. Then the set is extended by the time derivatives $\mathbb{A}_{2i}=\mathrm{i}\mathcal{L}\mathbb{A}_{1i}$. Similar as before, the effective masses $\mu_i=\sum_jm_j/\alpha^2_{ji}$ are defined. However, there is one additional restriction in this case: if a microscopic coordinate q_j enters one A_1k it must not enter any other A_1l . Otherwise, the partial integration performed to obtain eq. (12) will yield additional terms. If the coarse-grained observables are the centers of mass of different blobs, this means that a single particle must not be attributed to more than one blob.

In this case the derivation can be carried out as before and we obtain

$$\begin{align} \exp(t\mathrm{i}\mathcal{L})\mathcal{P}\mathrm{i}\mathcal{L}\mathbb{A}_{2i} &= -\frac{1}{\mu_i}\frac{\partial}{\partial A_{1i}}\underbrace{(-k_{\rm B} T\ln(\rho_A(\vec{A})))}_{=:W(\vec{A})}\Big|_{\vec{A}_t}. \end{align} \tag{ 14 }$$

Again, $\partial W(\vec{A})/\partial A_{1i}$ does not depend on A_2i .

Thus, the second term in eq. (5) takes the form of a derivative of a potential of mean force under the following conditions:

• The Hamiltonian is of the form given in eq. (8).
• A Zwanzig-type projector (cf. (6)) onto two observables is used.
• The first set of observables of the projector is of the form $\mathbb{A}_{1i}=\alpha_{0i}+\sum_j \alpha_{ji}q_j$ with constant $\alpha_{ji}$.
• If a coordinate q_j enters one A_1k it must not enter any other A_1l .
• The second set of observables of the projector is the time derivative of the first set $\mathbb{A}_{2i}=\mathrm{i}\mathcal{L}\mathbb{A}_{1i}$.

(We note that this derivation does also hold for a time-dependent Hamiltonian if it can be expressed in the form of eq. (8) with time-dependent masses and/or a time-dependent potential $\mathbb{V}(q_1,\ldots,q_N;t)$. In this case, we would need a time-dependent projector where the equilibrium density in eq. (6) is replaced by the equilibrium density with respect to the current Hamiltonian $\rho^\text{eq}(\Gamma;t)\propto\exp(-\beta\mathbb{H}(\Gamma;t))$. See ref. [22] for a suitable projection operator approach.)

Next, we consider the first and third term of eq. (5). Inserting the Zwanzig projector, eq. (6), into the first term of the right-hand side of eq. (5) we obtain a term which is in general non-linear in $\dot{A}_\tau$. Based on Zwanzig's work [25], Hijón et al. showed in ref. [20] that the memory term for the second component can be written as

$$\begin{align*} \int\limits_0^t\mathrm{d}s\!\sum_{i=1,2}\! \left(M_i(\vec{A}_s,t-s)\frac{\partial W(\vec{A}_s)}{\partial A_{i,s}}+k_{\rm B} T \frac{\partial M_i(\vec{A}_s,t-s)}{\partial A_{i,s}}\right) \end{align*}$$

with

$$\begin{equation*} M_i(\vec{A},t):= \frac{1}{k_{\rm B} T}\mathcal{P}\left(\left[\mathrm{i}\mathcal{L}\mathbb{A}_i\right]\left[\mathcal{Q}\exp(t\mathrm{i}\mathcal{L}\mathcal{Q})\mathrm{i}\mathcal{L}\mathbb{A}_2\right]\right). \end{equation*}$$

(The sum runs over the number of components of $\vec{\mathbb{A}}$.) This expression differs considerably from the desired one in eq. (2). Even if we assume that there is time-scale separation between the observables $\mathbb{A}_i$ and the other degrees of freedom, such that $M_i(\vec{A}_s,t-s)\propto \gamma(\vec{A}_s)\delta(t-s)$, we do not recover eq. (1), because in general, $\gamma(\vec{A}_s)$ is not a constant.

In order to see the structure of the first and third term of eq. (5) more clearly, we now use a mapping of the Zwanzig projector to a Mori projector as proposed in refs. [27,28]. (This mapping holds for a general set $\{\mathbb{A}_i\}$, not just for the specific one used in the previous paragraphs.) We define a scalar product between square-integrable phase-space fields by

$$\begin{align} \left(\mathbb{X},\mathbb{Y}\right) &= \int\mathrm{d}\Gamma\,\rho^{\rm eq}(\Gamma)\mathbb{X}(\Gamma)\mathbb{Y}(\Gamma), \end{align} \tag{ 15 }$$

and express the projector in eq. (6) as

$$\begin{align} \mathcal{P}\mathbb{X}(\Gamma) &= \left.\frac{(\delta(\vec{\mathbb{A}}-\vec{A}),\mathbb{X})}{(\delta(\vec{\mathbb{A}}-\vec{A}),1)}\right|_{\vec{A}=\vec{\mathbb{A}}(\Gamma)}. \end{align} \tag{ 16 }$$

The set of phase-space functions that depend on Γ solely through $\vec{\mathbb{A}}(\Gamma)$ are a closed subset of all phase-space functions. Thus, we can define a complete (infinite) set of phase-space functions $f_i(\vec{\mathbb{A}}(\Gamma))$ such that

$$\begin{align} (f_i(\vec{\mathbb{A}}),f_j(\vec{\mathbb{A}})) &= \delta_{i,j} \end{align} \tag{ 17 }$$

and

$$\begin{align} \sum\limits_{i=1}^\infty f_i(\vec{\mathbb{A}}(\Gamma))f_i(\vec{\mathbb{A}}(\Gamma')) &= \frac{\delta(\vec{\mathbb{A}}(\Gamma)-\vec{\mathbb{A}}(\Gamma'))}{(\delta(\vec{\mathbb{A}}-\vec{A}),1)\big|_{\vec{A}=\vec{\mathbb{A}}(\Gamma)}}. \end{align} \tag{ 18 }$$

These functions form a basis for the subspace of phase-space functions that depend on Γ solely through $\mathbb{A}_i$. In practice, such a set of basis functions can be obtained by means of a Gram-Schmidt process starting from monomials in $\mathbb{A}_i$. Note, that the denominator in eq. (18) is a mere consequence of the normalization of the basis functions. Thus, we can write the Zwanzig projector in eq. (16) as

$$\begin{align} \mathcal{P}\mathbb{X}(\Gamma) &= \sum\limits_{i=1}^\infty (f_i(\vec{\mathbb{A}}),\mathbb{X})f_i(\vec{\mathbb{A}}(\Gamma)). \end{align} \tag{ 19 }$$

However, this is nothing but a Mori projector on the infinitely many observables $f_i(\vec{\mathbb{A}})$. Using this expression, we can write the second term in eq. (5) as

$$\begin{equation*} \sum\limits_{i=1}^\infty\int\limits_0^t\mathrm{d}s\, K_{2,i}(t-s)f_i(\vec{A}_s), \end{equation*}$$

with

$$\begin{equation*} K_{j,i}(t-s) = (f_i(\vec{\mathbb{A}}),\mathrm{i}\mathcal{L}\mathcal{Q}\exp((t-s)\mathrm{i}\mathcal{L}\mathcal{Q})\mathrm{i}\mathcal{L}\mathbb{A}_j). \end{equation*}$$

Using the shorthand notation

$$\begin{align} \vec{\varepsilon}_t &= \mathcal{Q}\exp\left(t\mathrm{i}\mathcal{L}\mathcal{Q}\right)\mathrm{i}\mathcal{L}\vec{\mathbb{A}}, \end{align} \tag{ 20 }$$

we obtain the equation of motion

$$\begin{align} \frac{\mathrm{d}^2A_{1,t}}{\mathrm{d}t^2}&=\frac{\mathrm{d}A_{2,t}}{\mathrm{d}t}\nonumber\\ &=\sum\limits_{i=1}^\infty\int\limits_0^t\mathrm{d}s\, K_{2,i}(t-s)f_i(\vec{A}_s)\nonumber\\ &\quad-\frac{1}{\mu}\frac{\partial}{\partial A_1}W(A_1)\Big|_{A_1=A_{1,t}}+\varepsilon_{2,t}. \end{align} \tag{ 21 }$$

If, in particular, $\mathbb{A}$ is the center-of-mass position of a set of atoms, eq. (21) reads

$$\begin{align} \dot{x}_t & = -\frac{1}{\mu}\frac{\partial}{\partial x}W(x)\Big|_{x=x_t}\\ & +\sum\limits_{i=1}^\infty\int\limits_0^t\mathrm{d}s\, K_{2,i}(t-s)f_i\left(x_s,\dot{x}_s\right) +\varepsilon_{2,t}, \end{align}$$

which differs considerably from eq. (2). As the Liouvillian is anti-self-adjoint we can write the memory kernels as

$$\begin{align} K_{j,i}(t) &= -(\mathrm{i}\mathcal{L}f_i(\vec{\mathbb{A}}),\mathcal{Q}\exp\left(t\mathrm{i}\mathcal{L}\mathcal{Q}\right)\mathrm{i}\mathcal{L}\mathbb{A}_j). \end{align} \tag{ 22 }$$

If the phase-space distribution at time t equals $\rho^\text{eq}$, the scalar product can be interpreted as a correlation,

$$\begin{align} K_{j,i}(t) &= -\left\langle \left(\left.\frac{\mathrm{d}f_i(\vec{A}_t)}{\mathrm{d} t}\right|_{t=0}\right)\varepsilon_{j,t} \right\rangle\nonumber\\ &= -\left\langle \left(\vec{\varepsilon}_0\cdot\left.\frac{\mathrm{d}f_i(\vec{A})}{\mathrm{d} \vec{A}}\right|_{\vec{A}=\vec{\mathbb{A}}(\Gamma_0)}\right)\varepsilon_{j,t} \right\rangle, \end{align} \tag{ 23 }$$

i.e., there is a relation between correlations in the fluctuating force and the memory kernel, but all functions f_i (A_t ) enter this relation. This is as close as we get to a fluctuation-dissipation relation. Thus we conclude, if we enforce the drift term in the equation of motion to be a derivative of a potential of mean force, then the memory term and the fluctuating force term are not related by a fluctuation-dissipation relation.

Via the functions f_i the memory term of eq. (21) contains all powers and combinations of the variables $\mathbb{A}_i$, not just linear terms. To obtain an expression which is closer in structure to eq. (2) (i.e., one in which the integrand is linear in the observable), we begin the Gram-Schmidt procedure with the linear term $\mathbb{A}_2$ and ensure $f_1(\vec{A}_t)=\phi A_{2,t}$ where ϕ is the normalization factor. Then eq. (2) follows if we truncate the sum in eq. (21) at i = 1. However, this sum is in general not an expansion in a small parameter and, hence, we have no information on the magnitude of the other terms. Thus they should not be dropped without verifying that this constitutes a reasonable approximation for the specific system at hand. In the case where the fluctuations of the coarse-grained variables around their means are small, Kauzlarić et al. showed that the Zwanzig projector can be approximated by a "Mori-like" projector [29]. Under these conditions, the memory term can be approximated as one that is linear in the observable. However, the same applies to the drift term, therefore one then does not obtain an equation of motion that contains a non-linear generalized drift (which would allow to define a potential of mean force). A special exception to this is the case where the potential of mean force is quadratic such that its derivative is linear and coincides with the Mori drift term. However, this is certainly not the general case.

We conclude eq. (2) is neither exact nor, in general, the result of a controlled approximation. The authors of ref. [15] came to a different conclusion, because they switched between a Zwanzig projector and a Mori projector for a single variable (rather than the inifinitely many variables needed for eq. (19)) half-way through their derivation. In the work of Kinjo et al., the time dependences in eq. (17) of ref. [16] and eq. (26) of ref. [16] do not match up and, hence, the projector is implicitly switched as well. Unfortunately, in eq. (26) of ref. [16] the time dependences are not given expicitly.

Finally, we note that $\vec{\varepsilon}_t$, eq. (20), is orthogonal to any phase-space field $g(\vec{\mathbb{A}}(\Gamma))$ which depends on the phase-space coordinates solely through $\vec{\mathbb{A}}$. Thus, if the phase-space distribution equals $\rho^\text{eq}$ at all times,

$$\begin{align} \left\langle g_0 \vec{\varepsilon}_t\right\rangle &= \vec{0}\qquad\forall t. \end{align} \tag{ 24a }$$

Note, that also $\left\langle\vec{\varepsilon}_t\right\rangle=0\;\forall t$ as we could choose $g(\vec{\mathbb{A}})=\text{const}$.

Now we extend the discussion to full non-equilibrium, i.e., we allow for an explicit time dependence of the Liouvillian. To simplify the projection operator formalism, we "augment" the phase-space by one additional dimension (time) [22,27]. The new coordinates are $\Gamma^{\rm a} = (\Gamma,\tau)$, where the superscript a stands for "augmented phase-space". We denote the observable fields on the augmented phase-space by $\mathbb{A}^{\rm a}(\Gamma^{\rm a})$. However, $\mathbb {A}^\text{a}$ shall not depend on τ explicitly. The equivalent to the Liouville operator is

$$\begin{equation*} \mathrm{i}\mathcal{L}^{\rm a}\circ = \dot{\Gamma}^{\rm a}(\Gamma^{\rm a})\cdot\frac{\partial}{\partial \Gamma^{\rm a}}\circ \end{equation*}$$

and observables evolve according to the equation

$$\begin{equation} A_t(\Gamma^{\rm a}) = \exp\left(\mathrm{i}\mathcal{L}^{\rm a}t\right)\mathbb{A}^{\rm a}(\Gamma^a). \end{equation} \tag{ 25 }$$

We introduce an inner product on the augmented phase-space

$$\begin{eqnarray*} & \left(\mathbb{X}^{\rm a},\mathbb{Y}^{\rm a}\right)^{\rm a}_t = \int\mathrm{d}\Gamma^{\rm a} \;\rho^{\rm a}(\Gamma^{\rm a},t)\; \mathbb{X^{\rm a}}(\Gamma^{\rm a})\mathbb{Y^{\rm a}}(\Gamma^{\rm a})=\\ & \int \mathrm{d}\Gamma^{\rm a} \;\rho^{\rm a}(\Gamma^{\rm a},0) \left(\exp\left(\mathrm{i}\mathcal{L}^{\rm a}t\right)\mathbb{X}^{\rm a}(\Gamma^{\rm a})\right)\\ & \times\left(\exp\left(\mathrm{i}\mathcal{L}^{\rm a}t\right)\mathbb{Y}^{\rm a}(\Gamma^{\rm a})\right), \end{eqnarray*}$$

where the notation $\rho^{\rm a}(\Gamma^{\rm a},t)$ indicates that we synchronized the phase-space distribution such that $\rho^{\rm a}(\Gamma^{\rm a},t)= \rho(\Gamma,t)\delta(\tau-t)$. As above, we define an orthonormal basis $\{\phi^{\rm a}_i (\mathbb{A}^{\rm a},\tau)\}$ such that

$$\begin{equation*} \left(\phi^{\rm a}_i (\mathbb{A}^{\rm a},\tau),\phi^{\rm a}_j (\mathbb{A}^{\rm a},\tau)\right)_t=\delta_{i,j} \quad \forall t. \end{equation*}$$

Note that we will, in general, need a different set of basis functions for each time t. As a generalized version of the Zwanzig projector, we define

$$\begin{align*} \mathcal{P}(t)\mathbb{X}^\text{a}(\Gamma^\text{a})\! =\! \frac{\int \mathrm{d}{\Gamma^{\rm a}}' \rho^\text{a}({\Gamma^\text{a}}',t) \delta(\mathbb{A}^\text{a}({\Gamma^\text{a}}')\!-\!\mathbb{A}^\text{a}(\Gamma^\text{a}))\mathbb{X}^\text{a}({\Gamma^\text{a}}')}{\int \mathrm{d}{\Gamma^\text{a}}' \; \rho^\text{a}({\Gamma^\text{a}}',t) \delta(\mathbb{A}^\text{a}({\Gamma^\text{a}}')-\mathbb{A}^\text{a}(\Gamma^\text{a}))}. \end{align*}$$

Using the basis set, this projector can be brought into the form

$$\begin{align*} \mathcal{P}(t)\mathbb{X}^{\rm a}(\Gamma^{\rm a}) = \sum_{i=1}^\infty (\phi^{\rm a}_i (\mathbb{A}^{\rm a}(\Gamma^{\rm a}),\tau), \mathbb{X}^{\rm a}({\Gamma^{\rm a}}'))_t\, \phi^{\rm a}_i (\mathbb{A}^{\rm a}(\Gamma^{\rm a}),t). \end{align*}$$

In contrast to eq. (19), this expression is not a Mori projector on the augmented space. However, as it is linear in the functions $\phi^a_i$, it can still be inserted straight-forwardly into the Dyson-Duhamel identity. We obtain the equation of motion

$$\begin{equation} \frac{\mathrm{d}A_t}{\mathrm{d}t} = \sum_{i=0}^\infty\left(\omega_i(t)\phi^{\rm a}_i (A_t,t) + \int_{0}^t\mathrm{d}s\,K_i(t,s)\phi^{\rm a}_i (A_s,s)\right) + f_{t}, \end{equation} \tag{ 26 }$$

with

$$\begin{align*} &\omega_i(t) = (\phi^{\rm a}_i(\mathbb{A}^{\rm a}({\Gamma^{\rm a}}),\tau),\mathrm{i}\mathcal{L}^{\rm a}\mathbb{A}^{\rm a}({\Gamma^{\rm a}}))_t,\\ &K_i(t_1,t_2) = \\ &(\phi^{\rm a}_i(\mathbb{A}^{\rm a}({\Gamma^{\rm a}}),\tau),(\mathrm{i}\mathcal{L}^{\rm a}\!-\!\dot{\mathcal{P}}(t_2))\mathcal{Q}(t_2) \mathcal{G}(t_2,t_1)\mathrm{i}\mathcal{L}^{\rm a}\mathbb{A}^{\rm a}({\Gamma^{\rm a}}'))_{t_2}, \end{align*}$$

and

$$\begin{align} f_{t}=\mathcal{Q}(0)\mathcal{G}(0,t)\mathrm{i}\mathcal{L}^{\rm a}\mathbb{A}^{\rm a}({\Gamma^{\rm a}}'), \end{align} \tag{ 27 }$$

where $\mathcal{Q}(t) = 1-\mathcal{P}(t)$ and the negatively time-ordered exponential $\mathcal{G}(t_1,t_2)=\exp_-(\int_{t_1}^{t_2}\mathrm{d}s\,\mathrm{i}\mathcal{L}^a\mathcal{Q}(s))$.

If we again impose the condition that $\phi_1^{\rm a}(\mathbb{A}^{\rm a})\propto \mathbb{A}^{\rm a}$, we could, in principle, truncate the sum in eq. (26) at i = 1 in order to obtain a memory term linear in A_s . However, as above this sum is not an expansion in a small parameter thus the truncation does not produce a well-controlled approximation.

In summary, we discussed the structure of the non-linear, generalized Langevin equation for a set of coarse-grained observables. By means of a projection operator formalism, we showed that the widely used equation of motion, which consists of a derivative of a potential of mean force as the organized drift, a memory term which is linear in the observable and a fluctuating force which obeys a fluctutation-dissipation relation with respect to the memory kernel, is in general not exact.

To expand the memory kernel in a set of orthogonal polynomials and to then truncate this expansion after the linear contribution does not constitute a general pathway to a controlled approximation to the exact dynamics. Whether or not the combinaton of a potential of mean force with a linear memory term serves as a suitable approximation to a system's coarse-grained dynamics, therefore needs to be tested case by case.

The authors acknowledge funding by the Deutsche Forschungsgemeinschaft (DFG, German Research Foundation) –Project No. 430195928 and No. 431945604 (project P4 in FOR 5099). Further, the authors thank the Erwin Schrödinger Institute (ESI).

Data availability statement: No new data were created or analysed in this study.