2.1 Memory approximation

The physical theory of Brownian motion assumes that thermal fluctuations which occur due to the collision of the small, light particles with the larger Brownian particle happen at a time scale which is much shorter than the one describing the evolution of the larger particle. The consequence of this assumption is that the force acting on the Brownian particle at each time is completely uncorrelated with the history of its evolution. Mathematically speaking, this corresponds to the idea that the random force is a stochastic process which obeys the Markov property. We note that in the theory of Brownian motion, there is a further important assumption, which is that the behaviour of the heat bath is spatially uniform. This assumption is clearly justified for a particle within a fluid which is at rest macroscopically. However, the spatial uniformity may not be generally valid, since many coarse-grained DoF are tightly bound within a molecule or bulk solid, and so are in close contact with their particular environment, rather than a homogeneous random one.

Moreover, for coarse-graining of molecular systems, it is common to make an analogous time-scale separation assumption about the action of the neglected DoF on the coarse-grained DoF, i.e. that the discarded DoF fluctuate at a faster time-scale and so act on the CG DoF randomly with no correlations in time.

Setting aside the justification for the use of time-scale separation and spatial homogeneity, under these two assumptions the fluctuating forces ${\boldsymbol{\mathcal{F}}}$ are usually modelled as a white noise process, with an ACF which is assumed to take the form:

(2.1) \begin{equation} \beta\,\mathbb{E}[{\boldsymbol{\mathcal{F}}}_K(0)\otimes{\boldsymbol{\mathcal{F}}}_J(t)] = 2\,\Phi_{KJ} \delta\left({{t}}\right)\quad\text{where}\quad \Phi_{KJ} = \beta\int_{0}^{\infty} \mathbb{E}[{\boldsymbol{\mathcal{F}}}_K(0)\otimes{\boldsymbol{\mathcal{F}}}_J(\xi)] \mbox{d} \xi,\end{equation}

with ${\boldsymbol{\mathcal{F}}}$ being the fluctuating forces which act in the system at full resolution. Since the fluctuating forces are assumed to be stationary as a random process (in fact they are defined as averages over equilibrium ensembles, see [Reference Berne and Pecora4]), the terms $\Phi_{KJ}$ do not depend on the particular time origin [Reference Lei, Caswell and Karniadakis39, Reference Kinjo and Hyodo33], i.e. there is no time dependence in the definition of $\Phi_{KJ}$ given in (2.1). One of the problems with this description is the so-called plateau problem, first noticed by Kirkwood and Buff [Reference Kirkwood and Buff34]. The plateau problem stems from the fact that the integral definition of $\Phi_{KJ}$ is zero if the upper boundary is $+\infty$ , because of the decay of the negative part of the correlation function at infinity [Reference Helfand26]. For this reason, one approach has been to cut the domain of integration for the last integral in equation (2.1) at a certain time $\tau$ which should be large enough to include the relevant information in the ACF, but small compared to the time scale of the evolution of the CG DoF. By cutting the integral at a finite time, we are assuming that it remains constant for $t > \tau$ (i.e. the integral reaches a plateau). The problem comes from the fact that the intermediate time $\tau$ is not well defined within the theory and there is no method to estimate it a priori.

In the following, we will show that our derivation overcomes this issue by proposing an asymptotic methodology which allows us to use data to fit a functional form for the memory, and thereby derive an explicit form for the Markovian approximation directly from data.

2.2 Memory ansatz

Our analysis starts in a similar fashion by making an ansatz on the form of the memory kernel $ \boldsymbol{\mathcal{M}}(\boldsymbol Z_{},t)$ , albeit a less restrictive one. We assume that the memory kernel for the generic pair of beads K and J takes the general form

(2.2) \begin{equation} \mathcal{M}_{KJ}(\boldsymbol Z_{},t) = \mbox{exp}\left[{-\left({\frac{t}{\tau_{KJ}}}\right)^\alpha}\right]f_{KJ}(\boldsymbol Z_{},t),\end{equation}

where $\tau_{KJ}$ is a characteristic time, $\alpha$ is a constant, and $f_{KJ}(\boldsymbol Z_{},t)$ is a sufficiently smooth function on the interval $0\leq t < \infty$ . In general, each entry of the memory kernel is specific to each pair of beads. However, in the following, in order to simplify the notation, we will write

(2.3) \begin{equation} \boldsymbol{\mathcal{M}}(\boldsymbol Z_{},t) = \textrm{e}^{- \left({\frac t\tau}\right)^\alpha}f(\boldsymbol Z_{},t),\end{equation}

i.e. we will consider the case where all the entries in the memory kernel are equal. The general case can then be easily obtained by repeating the same arguments for all the possible KJ pairs. It is interesting to note that a similar ansatz was considered when using a MZ-derived Generalised Langevin Equation to study the dynamical behaviour of amorphous solids [Reference Cui, Gebbia, Romanini, Rudić, Fernandez-Perea, Bermejo, Tamarit and Zaccone12, Reference Cui, Gebbia, Tamarit and Zaccone13, Reference Cui, Milkus and Zaccone14, Reference Cui, Yang, Qiao, Jiang, Dai, Wang and Zaccone15, Reference Zaccone54].

While equation (2.3) represents a somewhat strict ansatz for the form of the memory kernel, it is consistent with similar quantities calculated in other systems [Reference Li, Lee, Darve and Karniadakis42], and moreover, our approach includes a possibly non-uniform dependence on the spatial position, in contrast to assumptions made in other applications of CG theory discussed above. Moreover, given the connection between the memory kernel and the ACF of equation (1.9), both $\tau$ and $\alpha$ can be extracted from the ACF, for example by fitting to numerical data. In particular, $\tau$ represents the time decay of the ACF and we will show in Section 3.3 its explicit calculation in a test system. We will also make the natural physical assumptions that $\boldsymbol{\mathcal{M}}(\boldsymbol Z_{},t)$ is maximal when $t=0$ , and that the parameters satisfy $\tau>0$ and $\alpha > 0$ .

Under the assumptions that $\alpha>0$ and $\tau>0$ , the exponential factor makes it possible to expand the memory kernel as an asymptotic series from which the behaviour close to the origin can be explicitly obtained. Here, we argue that this asymptotic series approach leads naturally to a form of the Markovian approximation commonly considered in CG theory. We note that the case where $\alpha = 0$ corresponds to algebraic decay of the memory kernel, and will require a different approach which we leave for future work.

In the following, we consider only the first integral of the memory term (i.e., equation (1.7)) since $\boldsymbol{\mathcal{M}}$ does not depend on the momenta for an appropriate CG mapping, see [Reference Di Pasquale, Hudson and Icardi18]. Furthermore, we make the memory term dimensionless by multiplying it by suitable constants:

(2.4) \begin{align} &\frac{M L_c }{\tau_P^2}\int_0^{\tilde{t}}\, \boldsymbol{\widetilde{\mathcal{M}}}\big(\widetilde{\boldsymbol Z_{}}(\tilde{t}-\tilde{s}),\tilde{s}\big) \cdot \widetilde{\textbf{P}}(\tilde{t}-\tilde{s}) \,\mbox{d}{\tilde{s}} \end{align}

\begin{align} &\qquad = \int_0^{\tilde{t}} \, \textrm{e}^{-\left({\frac{\tau_P}{\tau}}\right)^\alpha \tilde{s}^\alpha} \tilde{f}\big(\widetilde{\boldsymbol Z_{}}(\tilde{t}-\tilde{s}),\tilde{s}\big) \widetilde{\textbf{P}}(\tilde{t}-\tilde{s}) \,\mbox{d}{\tilde{s}} \nonumber\end{align}

where:

• M is the mean effective mass of the coarse-grained DoF;

• $L_c$ is a characteristic length whose specific definition is not important for what is discussed here (and can be taken as the average bead size, for instance);

• $\tau_P$ is a characteristic time, which can be interpreted as the decay of the autocovariance function of the momenta (ACM) of the beads. This quantity can be interpreted also as the characteristic time scale of the macroscopic system.

We indicate with $\widetilde{\cdot}$ the dimensionless quantities and we note that we can always rewrite the different terms in equation (2.4) as function of the dimensionless time $\tilde{t} = t/\tau_P$ , by writing $\boldsymbol Z_{}(t-s)=\boldsymbol Z_{}(\tau_P(\tilde{t}-\tilde{s}))=\widetilde{\boldsymbol Z_{}}(\tilde{t}-\tilde{s})$ , where we absorb $\tau_P$ in the definition of the function. We will also simplify the notation in equation (2.4) by defining the dimensionless quantity $\lambda=\left({\frac{\tau_P}{\tau}}\right)^\alpha$ . We want to highlight here the first result of the formulation presented here. The concept of the time-scale separation can be precisely and quantitatively defined by measuring it with the parameter $\lambda$ . From the definition of $\lambda$ , we can distinguish two cases:

1. (i) $\lambda \gg 1$ , then $\tau_P \gg \tau$ , that is to say that the macroscopic time scale is much larger than the microscopic one. The Markovian approximation is expected to hold.

2. (ii) $\lambda \approx 1$ , then $\tau_P \approx \tau$ and there is no separation between the macroscopic and the microscopic time-scale. The Markovian approximation is expected to fail and the full memory term must be considered.

In the next section, we will present a result which will make more clear the qualitative explanation given here in terms of validity or not of the Markovian approximation. Then, the problem will become the explicit calculation of the two quantities, $\tau$ and $\tau_P$ , which will be presented in the following sections. In order to simplify our notation, we will drop tildes from now on, and assume all quantities are dimensionless.

2.3 Rigorous approximation result

Under the assumption that the ansatz equation (2.3) holds, we now derive an error estimate for an expansion of the memory integral which we can use to define proxy dynamics, providing a theoretical guarantee of accuracy in the case where the time-scale separation is significant.

Proposition 1. Suppose that the memory kernel $\mathcal{M}$ takes the form

\begin{equation*}\mathcal{M}(\boldsymbol Z_{}(t-s),s) = \textrm{e}^{- \lambda s^\alpha}f(\boldsymbol Z_{}(t-s),s),\end{equation*}

where $\lambda>0$ , $\alpha>0$ , and the function

\begin{equation*}g(t,s) = f(\boldsymbol Z_{}(t-s),s)\textbf{P}(t-s)\end{equation*}

is of class $\textrm{C}^{k,\gamma}$ in both t and s, for some non-negative integer k and Hölder exponent $\gamma\in(0,1]$ . Then the absolute error committed by replacing the memory

\begin{equation*}\int_0^t \mathcal{M}(\boldsymbol Z_{}(t-s),s)\textbf{P}(t-s) \mbox{d}{s} = \int_0^t \textrm{e}^{- \lambda s^\alpha} g(t,s) \mbox{d}{s}\end{equation*}

by a k-term approximation is bounded by

\begin{equation*} \Bigg|\int_0^t \textrm{e}^{- \lambda s^\alpha} g(t,s)\mbox{d}{s} - \sum_{n=0}^{k} \frac{G_n(t)}{n!}\frac{\Gamma\big(\frac{n+1}{\alpha}\big)}{\alpha \lambda^{\frac{n+1}{\alpha}}}\Bigg| \leq \sum_{n=0}^{k}\frac{C_n}{n!}\int_t^\infty \textrm{e}^{-\lambda s^\alpha} s^n \mbox{d}{s}+C_{k+1}\frac{\Gamma\Big(\frac{k+\gamma+1}{\alpha}\Big)}{k!\alpha \lambda^{\frac{k+\gamma+1}{\alpha}}},\end{equation*}

where

\begin{equation*}G_i(t) = \frac{\partial^i}{\partial s^i}g(t,0)=\frac{\partial^i}{\partial s^i}\Big(f(\boldsymbol Z_{}(t-s),s)\textbf{P}(t-s)\Big)\bigg|_{s=0},\end{equation*}

$\Gamma$ is the Gamma function, and $C_n$ are appropriate positive constants which depend upon the regularity of the trajectories.

A full proof of this proposition is given in Appendix A, and we focus on the physical interpretation of the result here. The error bound for the memory matrix is composed of two terms. The first of these is

(2.5) \begin{equation} \psi(t,\lambda) \overset{\operatorname{def}}{=} \sum_{n=0}^{k}\frac{C_n}{n!}\int_t^\infty \textrm{e}^{-\lambda s^\alpha} s^n \mbox{d}{s},\end{equation}

and it can be shown that $\psi$ decays exponentially in time for fixed $\lambda$ (see Appendix B). Physically, this decay represent transient behaviour as the system relaxes towards the equilibrium. We believe this transient behaviour is closely related to the plateau problem which has been observed in earlier works referred to above.

The second term

(2.6) \begin{equation} \zeta(\lambda) \overset{\operatorname{def}}{=} C_{k+1}\frac{\Gamma\Big(\frac{k+\gamma+1}{\alpha}\Big)}{k!\alpha \lambda^{\frac{k+\gamma+1}{\alpha}}}\end{equation}

does not decay with increasing time t. However, we can study its behaviour with respect to $\lambda$ . In particular, it follows that $\zeta(\lambda)$ decays to zero as $\lambda\to\infty$ because $\frac{k+\gamma+1}{\alpha} > 0$ , then

(2.7) \begin{equation}\zeta(\lambda) = \mathcal{O}\left({\lambda^{-\frac{k+\gamma+1}{\alpha}}}\right) \,, \quad \, \lambda\to\infty.\end{equation}

As already observed (see Section 2.2), $\lambda$ represents the scale-separation between the underlying atomistic model and the coarse-grained one, so as we will show below, the above result provides both a form for a Markovian approximation of the dynamics, and a theoretical guarantee backing up the commonly held belief that a Markovian approximation is more accurate when there is a large time scale separation between atomistic and CG models.

In particular, as a corollary of Proposition 1, we note that if the integrand function g(t,s) in the memory integral is Lipschitz, so that $k=0$ and $\gamma=1$ , then we have

\begin{equation*}\Bigg|\int_0^t \textrm{e}^{- \lambda s^\alpha} g(t,s)\mbox{d}{s} - G_0(t)\frac{\Gamma\big(\frac{1}{\alpha}\big)}{\alpha \lambda^{\frac{1}{\alpha}}}\Bigg|\leq C_0\int_t^\infty \textrm{e}^{-\lambda s^\alpha}\mbox{d}{s}+C_1\frac{\Gamma\big(\frac2\alpha\big)}{\alpha \lambda^{\frac2\alpha}}.\end{equation*}

We remark that the assumption that g(t,s) is Lipschitz is particularly mild. It is guaranteed to hold if trajectories of the system are twice continuously differentiable in time, and this will in turn be the case when the potential energy function is twice continuously differentiable in space. As such, the assumption is therefore satisfied for almost all Hamiltonian systems of interest.

We note that the function $G_0(t)$ is

\begin{equation*} G_0(t) = f(\boldsymbol Z_{}(t),0)\textbf{P}(t) = \beta\, \mathbb{E}\left[ {{ {\boldsymbol{\mathcal{F}}}(0,\cdot) \otimes {\boldsymbol{\mathcal{F}}}(0,\cdot) }} \middle\vert {{ \,\boldsymbol Z_{}(t) }} \right]\textbf{P}(t),\end{equation*}

and hence we are able to approximate equation (2.4), where we have dropped the tildes, as

(2.8) \begin{align} \int_0^t \mathcal{M}(\boldsymbol Z_{}(t-s),s)\textbf{P}(t-s) \mbox{d}{s} & \approx \frac{\Gamma\big(\frac{1}{\alpha}\big)}{\alpha \lambda^{\frac{1}{\alpha}}} \beta\, \mathbb{E}\left[ {{ {\boldsymbol{\mathcal{F}}}(0,\cdot) \otimes {\boldsymbol{\mathcal{F}}}(0,\cdot) }} \middle\vert {{ \,\boldsymbol Z_{}(t) }} \right]\textbf{P}(t), \nonumber \\[5pt] & \equiv \boldsymbol{\chi}(\boldsymbol Z_{}(t);\alpha,\lambda )\textbf{P}(t).\end{align}

or in other words, we can define a spatially varying friction cooefficient matrix $\boldsymbol{\chi}$ to replace the memory integral.

In this case, we can use this expression to show that the relative error of the approximation is bounded by

(2.9) \begin{equation} \Bigg|\frac{\int_0^t \textrm{e}^{- \lambda s^\alpha} g(t,s)\mbox{d}{s} - G_0(t)\Gamma\big(\frac{1}{\alpha}\big)\alpha^{-1} \lambda^{-\frac{1}{\alpha}}}{G_0(t)\Gamma\big(\frac{1}{\alpha}\big)\alpha^{-1} \lambda^{-\frac{1}{\alpha}}}\Bigg| \leq \frac{\int_t^\infty \textrm{e}^{-\lambda s^\alpha}\mbox{d}{s}}{\int_0^\infty \textrm{e}^{-\lambda s^\alpha} \mbox{d}{s}}+\frac{C_1}{\lambda^{\frac1\alpha}|G_0(t)|}\frac{\Gamma\big(\frac2\alpha\big)}{\Gamma(\frac1\alpha)}.\end{equation}

As we already observed for the absolute error bound, the former term exhibits decays to zero as $t\to\infty$ , regardless of the value of $\lambda$ , while the latter decays as $\lambda\to\infty$ for fixed t.

In the next section, we will show how to compute $\chi$ in practice and discuss the error committed in replacing the memory term with its approximation shown in equation (2.8).

3.1 Model

We assume the atoms are arranged in a repeating pattern, and a coarse-grained model is obtained by combining the single repeated units into beads. The CG variables are given by the centre of mass of each bead and by the corresponding momenta, and the repeating pattern consists of two atoms of mass $M_2$ and a single atom of mass $M_1$ arranged in the following configuration: $(M_2-M_1-M_2)$ . We considered two test cases: $M_1=M_2=1$ and $M_1=1$ , $M_2=100$ .

The bond stiffness between the two different species of atoms is described by means of an inter-atomic potential, which is a simple 12-6 Lennard–Jones potential:

\begin{equation*} U_{i,i+1}(r) = 4\,\epsilon_{i,i+1} \left[{\left({\frac{\sigma_{i,i+1}}{r}}\right)^{12} - \left({\frac{\sigma_{i,i+1}}{r}}\right)^6}\right] = \epsilon_{i,i+1} \left[{\left({\frac{r^*_{i,i+1}}{r}}\right)^{12} - 2\left({\frac{r^*_{i,i+1}}{r}}\right)^6} \right]\,.\end{equation*}

Here, $i \in \left\lbrace {1,2,\dots,N^{FG}}\right\rbrace$ and we choose $\sigma_{i,i+1}=2^{-1/6}$ so that the minimum value of the potential is obtained for $r^*_{i,i+1}=1$ . The value of the interaction strength, $\epsilon_{i,i+1}$ , is set to 1 and 10 for $M_1 - M_1$ and $M_1 - M_2$ pairs, respectively. In the following, all parameters, such as temperature and time step, are expressed in terms of Lennard–Jones reduced units [Reference Tuckerman48].

The simulations to obtain the fluctuating force (and therefore to determine the memory kernel) are implemented in the canonical (NVT) ensemble at fixed $k_B T = 1$ and friction parameter $\gamma=1$ . Furthermore, we assume periodic boundary conditions and nearest neighbour interactions.

The Hamiltonian of the atomistic system is given by

(3.1) \begin{equation} \mathcal{H}(\textbf{r},\textbf{p}) = \frac12 \,\textbf{p}^T\,\textbf{M}^{-1}\,\textbf{p}\,+\sum_{i=1}^{N^{FG}-1}U_{i,i+1}\left({r_{i+1}-r_i}\right)\,+ U_{N^{FG},1}\left({r_1-r_{N^{FG}}+r^*N^{FG}}\right), \end{equation}

where $\textbf{M}$ is the mass matrix (that is, a diagonal matrix whose entries are the masses of each atom) and $\textbf{r}$ and $\textbf{p}$ are the positions and momenta of the atoms, respectively. The final term takes into account that periodic boundary conditions are in place, i.e., that the two ends of the linear chain are connected to each other.

The results shown below were obtained through simulation of two different dynamics: (i) the Fine-Grained Dynamics (FGD) when solving the equations of motion with the Hamiltonian defined in equation (3.1),

\begin{equation*} \frac{\mbox{d}{\textbf{z}(t)}}{\mbox{d}{t}} = \mathfrak{L}\,\textbf{z} \,, \end{equation*}

where, $\textbf{z}=\left({ \textbf{z}_{ \textbf{r} },\textbf{z}_{ \textbf{p} } }\right)=( \textbf{r}_1,\ldots,\textbf{r}_{N^{FG}} , \textbf{p}_1,\ldots,\textbf{p}_{N^{FG}} )$ is a state of the system characterised by the instantaneous positions and momenta of the $N^{FG}$ particles that compose it; (ii) the Orthogonal Dynamics (OD) when solving the system of equations

\begin{equation*} \frac{\mbox{d}{r_k}}{\mbox{d}{t}} = \frac{p_k}{m_{k}}- \frac{P_I}{M_I} \,; \qquad \frac{\mbox{d}{p_k}}{\mbox{d}{t}} = - \frac{\partial{U(\boldsymbol\zeta_{\textbf{r}})}}{\partial{r_k}} + \frac{m_k}{M_I}\sum_{i\in S_I}\frac{\partial{U(\boldsymbol\zeta_{\textbf{r}})}}{\partial{r_i}}\,; \end{equation*}

for $k = 1,2,\dots,N^{FG}$ , where I is the index such that $k \in S_I$ and $S_I$ is the set of particles of the FG system included within the bead of index I. We note that both FGD and OD have roughly the same numerical performances, and that for the simple system investigated here the speed-up achieved by the coarse-graining procedure is $\approx\!5$ .

3.2 Numerical sampling of the memory kernel

All the simulations reported below were performed using the code HybridZwanzJulia [Reference Icardi, Hudson, Di Pasquale, Spinaci and Rovigatti29] written in Julia v1.4.1 [Reference Bezanson, Edelman, Karpinski and Shah5]. This code allows a relatively simple recalibration of the test system through Julia ‘types’, which are initialised by the user to fix all the parameters of the simulation, including the repeating pattern of the beads, the number of beads, the mass and stiffness of the inter-atomic potentials. Here, we give an idea of the sampling implemented in the code to compute the fluctuating force, given in equation (1.6), and therefore the memory term as given in equation (2.8).

As mentioned in the introduction, the idea is to sample using OD, since it is possible to determine the fluctuating force through constrained dynamics. First, the initial conditions for the simulation are fixed, so that the momenta are distributed according to the Maxwell–Boltzmann distribution and the positions lie in a neighbourhood of the equilibrium positions, determined by the equilibrium distance of the inter-atomic potential. Then the trajectory of the FGD is computed and samples of the positions and momenta are stored. These samples are used to determine the CG positions and momenta and are used as initial condition for computing the trajectories of the OD. Along these trajectories, time averages of the first and second moment of the effective force between beads are stored. The first moment gives the value of the mean force acting between adjacent beads given the position of the center of mass, while the second moment is used to compute the fluctuating force. Since the sampling steps of the OD can run at the same time without having to communicate between the processes, these procedures can be run in parallel.

In order to improve the time required to sample the fluctuating force, we use the translational symmetry of this simple system, the fact that the beads are identical and the assumption that the effective potential can be approximated as a sum of two-body interactions based on the distance to the nearest-neighbour beads only. We thus compute the interactions between any combination of beads and assume that the interactions computed are valid for any equivalent pair in the system. In this way, we reduce the dimensionality of the space to be sampled from the total number of CG variables to simply one, i.e., the inter-bead distance, which simplifies the exploration of the entire constrained phase space.

If we consider bead I, we can assume that the forces that act on it can be decomposed as the independent ‘stress’ contributions $\sigma_{I+1,I}$ and $-\sigma_{I,I-1}$ , which are the forces exerted on bead I by its two neighbours, $I-1$ and $I+1$ . Under these conditions, we write the conditional expectation in the equation (2.8) as

\begin{align*} \mathbb{E}\left[ {{{\boldsymbol{\mathcal{F}}}_I \otimes {\boldsymbol{\mathcal{F}}}_I}} \middle\vert {{\textbf{R}}} \right] &= \mathbb{E}\left[ {{ (\Delta\sigma_{I+1,I}-\Delta\sigma_{I,I-1}) \otimes (\Delta\sigma_{I+1,I}-\Delta\sigma_{I,I-1}) }} \middle\vert {{\textbf{R}}} \right] \\[5pt] &= \mathbb{E}\left[ {{ \Delta\sigma_{I+1,I} \otimes \Delta\sigma_{I+1,I} }} \middle\vert {{\textbf{R}}} \right] + \mathbb{E}\left[ {{ \Delta\sigma_{I,I-1} \otimes \Delta\sigma_{I,I-1} }} \middle\vert {{\textbf{R}}} \right] \,,\end{align*}

where $\Delta\sigma$ is the fluctuating part of the ‘stress’ contribution, $\textbf{R}$ is the inter-bead distance and we assume that the cross terms are negligible, i.e.,

\begin{equation*} \mathbb{E}\left[ {{ \Delta\sigma_{I+1,I} \otimes \Delta\sigma_{I,I-1} }} \middle\vert {{\textbf{R}}} \right]\approx 0 \quad\text{and}\quad \mathbb{E}\left[ {{ \Delta\sigma_{I,I-1} \otimes \Delta\sigma_{I+1,I}}} \middle\vert {{\textbf{R}}} \right]\approx 0\,.\end{equation*}

In other words, we assume that the fluctuating stresses acting between adjacent pairs of beads are not correlated, since we have made the assumption that pair interactions are sufficient to describe the behaviour of the system.

Under this assumption, we obtain

\begin{align*} \mathbb{E}\left[ {{{\boldsymbol{\mathcal{F}}}_I \otimes {\boldsymbol{\mathcal{F}}}_{I+1}}} \middle\vert {{\textbf{R}}} \right] &= \mathbb{E}\left[ {{ (\Delta\sigma_{I+1,I}-\Delta\sigma_{I,I-1}) \otimes (\Delta\sigma_{I+2,I+1}-\Delta\sigma_{I+1,I}) }} \middle\vert {{\textbf{R}}} \right] \\[5pt] &= -\mathbb{E}\left[ {{ \Delta\sigma_{I+1,I} \otimes \Delta\sigma_{I+1,I} }} \middle\vert {{\textbf{R}}} \right] \,, \\[-5pt] \\ \mathbb{E}\left[ {{{\boldsymbol{\mathcal{F}}}_I \otimes {\boldsymbol{\mathcal{F}}}_{I-1}}} \middle\vert {{\textbf{R}}} \right] &= \mathbb{E}\left[ {{ (\Delta\sigma_{I+1,I}-\Delta\sigma_{I,I-1}) \otimes (\Delta\sigma_{I,I-1}-\Delta\sigma_{I-1,I-2}) }} \middle\vert {{\textbf{R}}} \right] \\[5pt] &= -\mathbb{E}\left[ {{ \Delta\sigma_{I,I-1} \otimes \Delta\sigma_{I,I-1} }} \middle\vert {{\textbf{R}}} \right] \,.\end{align*}

Therefore, we compute a smooth approximation of the variance of the fluctuating force between beads which allows us to define

(3.2) \begin{equation} \chi_{I+1,I} = \frac{\tau}{\tau_P}\,\frac{\Gamma\left(\tfrac{1}{\alpha}\right)}{\alpha}\,\beta\,\mathbb{E}\left[ {{ \Delta\sigma_{I+1,I} \otimes \Delta\sigma_{I+1,I} }} \middle\vert {{ \textbf{R} }} \right] \,,\end{equation}

and the spatially varying friction coefficient matrix $\boldsymbol{\chi}$ (see equation (2.8)) as

\begin{equation*}\chi_{IJ}(\textbf{R}; \alpha, \tau) = \begin{cases} -\chi_{I+1,I} & I=J-1 \mod N^{CG} \,, \\[5pt] \chi_{I+1,I} + \chi_{I,I-1} & I=J \,, \\[5pt] -\chi_{I,I-1} & I=J+1 \mod N^{CG} \,, \end{cases}\end{equation*}

where $I,J \in \left\lbrace {1,2,\dots,N^{CG}}\right\rbrace$ . We note here that the memory term we consider depends on CG coordinates. The importance of this feature in this kind of models was shown in [Reference Hudson and Li28].

3.3 Fitting and results

In this section, we show the numerical results for the autocorrelation function (ACF) of the fluctuating forces and of the bead momenta.

Applying our ansatz (see equation (2.3)), we fit the decaying of the ACF of the fluctuating forces obtained from a fine-grained simulation using the function

(3.3) \begin{equation} F_F(t;\alpha_F,\tau) = \mbox{exp}\left[ -\left({\frac{t}{\tau}}\right)^{\alpha_F} \right].\end{equation}

The best fit values for $\alpha_F$ and $\tau$ , obtained with the XMGRACE software [Reference Turner49], are reported in Table 1 for the two systems investigated here. In order to check the consistency of the results, we also fit the initial decay of the ACF of the momenta and of the fluctuating forces using

(3.4) \begin{equation} F_P(t;\tau_P) = \mbox{exp}\left[ -\left({\frac{t}{\tau_P}}\right)^2 \right]\,.\end{equation}

For comparison, the data and the accompanying fitted curves are shown in Figure 1.

Table 1.

Results of fitting to the numerical data

Figure 1.

Time autocorrelation function of the (a) fluctuating force and (b) momenta for two values of $M_2$ . Symbols represents the simulations data whereas dotted lines and dashed lines are fits to equations (3.4) and 3.3, respectively. In the ACF of the momenta (b) only one fitting is reported (solid black line) because equations (3.4) and (3.3) coincide.

For the fluctuating force, for both systems ( $M_2 = 1$ and $M_2 = 100$ ), the initial behaviour can be described by a quasi-Gaussian behaviour with a characteristic decay time that does not depend on the exponent used, followed by damped oscillations that we do not take into account in our analysis. By contrast, the bead momenta ACF of both systems exhibit a nearly-perfect Gaussian decay. This distinction is likely due to the quadratic nature of the kinetic energy term in the Hamiltonian.

We can now use equation (3.2) to obtain the elements of $\boldsymbol{\chi}$ , which explicitly depends on the inter-bead distance $\textbf{R}$ and are used to build our approximation of the memory term equation (2.8). The results for the two systems investigated here are reported in Figure 2. The numerical data shows that, while the fluctuating forces ${\boldsymbol{\mathcal{F}}}(0,\cdot)$ (shown in the inset of Figure 2) depends only weakly on the atomic masses, the resulting friction term features a much stronger dependence on $M_2$ through the value of $\tau$ . This formulation demonstrates the possibility of using the result of Proposition 1 to provide a CG approximation of the full system with appropriate theoretical guarantees on the error committed in replacing the memory integral with a simpler friction matrix.

Figure 2.

The friction coefficients evaluated for two values of $M_2$ . The inset reports the fluctuating force for the same systems. The curves are spline interpolations of the numerical data, which is shown with symbols in the inset.

3.4 Discussion

We can now use the numerical results reported in Table 1 to estimate the the error committed in replacing the full memory term with its approximated value based on the ansatz reported in equation (2.3). In particular, this error is controlled by the two terms $\psi(t,\lambda)$ and $\zeta(\lambda)$ (see equations (2.5), (2.6)). Using the results obtained for equations (B2), (2.7), we can write for the two terms of the boundary of the relative error in equation (2.9):

(3.5a) \begin{align} \psi^\prime(\tilde{t},\lambda) = \frac{\int_{\tilde{t}}^\infty \textrm{e}^{-\lambda \tilde{s}^\alpha}\mbox{d}{\tilde{s}}}{\int_0^\infty \textrm{e}^{-\lambda \tilde{s}^\alpha} \mbox{d}{\tilde{s}}} = \mathcal{O}\left({\textrm{e}^{-\lambda \tilde{t}^{\alpha}}}\right) = \mathcal{O}\left({\mbox{exp}\left[{-{\frac{t}{\tau}}^{\alpha}}\right]}\right) \,, \quad \tilde{t}\to\infty\, \end{align}

(3.5b) \begin{align} \zeta^\prime(\lambda) = \frac{C_1}{\lambda^{\frac1\alpha}|G_0(\tilde{t})|}\frac{\Gamma\big(\frac2\alpha\big)}{\Gamma(\frac1\alpha)} = \mathcal{O}\left({\lambda^{-\frac{1}{\alpha}}}\right) = \mathcal{O}\left({\frac{\tau}{\tau_P}}\right) \,, \quad \lambda\to\infty\,. \end{align}

In order to check the magnitude of the leading term in equation (3.5a), we need the time t, which represents the total simulated time of our simulation. However, for the sake of giving an estimate of the order of magnitude of $\psi^\prime(\tilde{t},\lambda)$ , we can simply use a time of the order of time scales considered in a standard CG simulation. For this reason, it seems reasonable to assume for it the typical time scale of the bead dynamics given by $\tau_P$ . Using the numbers reported in Table 1, the leading term in equation (3.5a) is $\mbox{exp}\left[{-\left({\frac{\tau_P}{\tau}}\right)^\alpha}\right] \approx 10^{-21}$ for the system with $M_2 = 1$ , and $\mbox{exp}\left[{-\left({\frac{\tau_P}{\tau}}\right)^\alpha}\right] \approx 10^{-8}$ for the system with $M_2 = 100$ . Although the specific values of these numbers (and even their orders of magnitude) depend on the functional form used to fit the simulation data, they are in all cases always much smaller than unity, demonstrating that the portion of the relative error bounded by $\psi^\prime(t)$ is extremely small. Moreover, the simulated time is usually well beyond $\tau_P$ , that is to say that in a standard simulation we can expect an even lower value of $\psi^\prime(\tilde{t},\lambda)$ .

The term $\psi^\prime(t)$ is roughly related to the replacement of the integration boundary in equation (2.4) from (0,t) to $(0,\infty)$ . It is therefore expected that the error due to this replacement decreases with the increase of t and, given the argument of the integral, it decreases with exponential behaviour.

Let us now consider the portion of the relative error bounded by $\zeta^\prime(\lambda)$ , equation (3.5b). With the same values reported in Table 1 we obtain $\frac{\tau}{\tau_P} = 0.225$ for the system with $M_2 = 1$ and $\frac{\tau}{\tau_P} = 0.29$ for the system with $M_2 = 100$ . The term $\zeta^\prime(\lambda)$ represents the relative error committed by replacing the memory term by its Taylor expansion stopped at a some k. In this case, we considered $k=0$ , i.e. we replaced the whole memory term with the leading term of its Taylor expansion. We see that, in this case, the bounding value for the percentage error is $\approx$ 23% for the system with $M_2 = 1$ and $\approx$ 30% for the system with $M_2 = 100$ , meaning that the results of the simulations should be carefully checked as it is not guaranteed that the error committed would be negligible.

The analysis just showed puts on a more rigorous ground the qualitative statement that the characteristic time of the decay of the momenta ACF, being of the same order of magnitude of the decay time of the fluctuating-force ACF, makes resorting to equation (3.2) questionable. Moreover, our results show that including more terms in the approximation (see Section 2.3) should improve such model. Indeed, since in equation (2.7) $\zeta(\lambda)$ approaches zero faster upon the increasing of k, this error bound can be reduced, making the results of such simulations more robust. We plan to test this approach in future work.