Markovian equations of motion for non-Markovian coarse-graining and properties for graphene blobs

We obtain Markovian equations of motion for a many body system of interacting coarse-grained (CG) variables and additional fluxes. The investigated CG variables belong to the special family of linear combinations of atomistic degrees of freedom. The system of Markovian equations of motion approximates Mori's exact non-Markovian generalized Langevin equation and is easy to solve by computer simulation. All parameters of the equations can be obtained from equilibrium molecular dynamics simulations of the investigated microscopic system. These parameters are either equal to the famous static covariances from Mori's continued fraction or they represent generalized constant friction matrices. We propose two different ways to compute these friction matrices based on Mori's continued fraction. Finally, some of the parameters are computed numerically for the special case of centre of mass variables in the graphene lattice and it is found that the CG variables interact with their additional fluxes in a spatially very local way.


Introduction
The term 'coarse-graining' (CG) refers to a collection of methods used to reduce the number of degrees of freedom of complex atomistic systems such as polymers [1,2], proteins [3], fluids [4][5][6][7][8] or combinations thereof. The motivation for doing so is straightforward: computational time can be saved or larger systems over longer time-scales can be approached by keeping only the relevant part of the information. The hope is that a few relevant variables still reproduce the relevant predictions of the complex atomistic model.
Systematically coarse-graining Hamiltonian systems containing many microscopic variables with the Mori-Zwanzig projection operator formalism [9][10][11][12] leads to non-Markovian equations of motion of the type of a generalized Langevin equation (GLE) for so-called relevant variables. In many cases, by a clever choice of relevant variables, a Markovian approximation can be justified.
A prominent mesoscale method is dissipative particle dynamics (DPD) [13,14]. In the DPD-method it is usually assumed a priori that the equation of motion is Markovian and hence the noise is white. In general the noise is not white but coloured. Occurrences or applications of coloured noise include the bath of harmonic oscillators [15], extensions of the Kramers theory of chemical relaxation [16] or constant temperature molecular dynamics (MD) [17,18].
DPD is often used on an empirical basis with empirically fitted parameters and no rigorous justification of the Markovianity. For simple fluids this approach is very well justified, since there is a clear separation of timescales between slow macroscopic and fast microscopic variables [10]. The Navier-Stokes equations themselves are clearly Markovian. In lattice dynamics or for the CG of large molecules, the non-Markovianity deserves further investigations and it was argued that the Markovian assumption introduces significant errors into the coarsegrained description [19][20][21][22].
The work presented here is motivated by earlier work on the development of Markovian CG models for graphene and carbon nanotubes (CNTs). Early work was based on top-down CG CNT-models assuming a priori the DPD-form of the equations of motion for centre of mass (COM) variables of groups of carbon atoms [23,24]. In a systematic bottom-up approach for two-dimensional (2D) motion of graphene, a Markovian approximation was applied to GLEs obtained from Mori's or Zwanzig's projection operators [25][26][27]. COM-based relevant variables seemed to show not only less non-Markovian features than finite element based variables but also showed weaknesses in the reproduction of cross-correlation functions of the relevant variables or long-wavelength normal mode decay, while the reproduction of autocorrelation functions is nearly perfect. The finite element based variables obtained by the CG procedure 'CGMD' [28] clearly showed non-Markovian behaviour. Both choices of variables have in common that they belong to the family of linear combinations of the atomistic variables. Therefore, this work aims for the development of general CG non-Markovian equations of motion for this family of CG variables and for schemes for the approximate numerical computation of these equations.
It was already shown many times for a single scalar relevant variable that Mori's hierarchy of non-Markovian GLEs (or equivalently the continued fraction representation) [29] can be reformulated into a system of Markovian linear equations with auxiliary variables which is much easier to solve if truncated appropriately [30,31]. In the work of Darve et al [32], a Galerkin discretization of a generalized Fokker-Planck equation is introduced, using a set of basis functions of the relevant variables. Also this approach was applied to one single COM variable of a single protein or to a single reaction coordinate. Our aim is a formulation of a CG many body system in 2D or 3D with fewer degrees of freedom than in the microscopic description, but still with many coupled CG variables, where a direct discretization of the phase space of CG variables becomes impractical due to its high-dimensionality. For this purpose we find the equations of motion derived by Karasudani et al [33] a very useful starting point. They involve a change in perspective by extending directly the set of relevant variables entering the GLE. The additional variables that are used in this approach are exactly those that appear in Mori's hierarchy, but now propagated with a computable dynamics described by the nonprojected Liouvillian. In contrast to Mori's continued fraction in the Laplace domain, this exact transformation provides equations that allow one to combine the two steps of truncation of the hierarchy and to find an extended Markovian system in the time domain. It is a description of the CG dynamics in the time-domain, which we eventually require, in order to implement it into a practical computational tool.
In section 2, starting from Mori's GLE and its generalization to an infinite hierarchy of GLEs for the noises [29], we obtain useful approximations for the truncation of this hierarchy by constant generalized friction matrices and we derive a finite set of Markovian equations of motion for the main relevant variables and auxiliary fluxes. The latter is achieved by making use of Karasudani's approach to Mori's GLE [33]. Then, in section 3, we apply the method to relevant variables represented by the family of linear combinations of microscopic atomistic degrees of freedom. In this way we obtain a clear prescription for how to compute all coefficients in the resulting set of Markovian equations from equilibrium MD simulations. Since the subject is very algebra intensive, the main results of sections 2 and 3 are completed by a few more details provided in the appendix and all the details of the calculations provided as supplementary material (available from stacks.iop.org/NJP/15/125015/mmedia). Finally, section 4 gives some basic numerical results for the special case of hexagonal COM variables constructed from the 2D in plane motion of a simple atomistic graphene model. We observe that the interaction of the CG displacements and momenta with the auxiliary fluxes is very local in nature, i.e. the CG momentum of a COM-blob µ interacts mostly with its own auxiliary fluxes and with its nearest neighbours ν in the graphene lattice. This property is advantageous from the computational point of view.

Mori's generalized Langevin equation (GLE)
For a set of relevant variables A 0 (t) we obtain from Mori theory [12,34] the GLE, which is a set of exact equations of motion of the forṁ The need for the index 0 will become clear later. We use boldface notation because, in general, we are interested in a vector A 0 containing a set of say M relevant variables, which may be vectors themselves in a D-dimensional space with D > 1, i.e. we might write A 0 ≡ A 0µ ≡ A 0µα with µ = 1 . . . M and α = 1 . . . D. We will use Greek 'blob' indices such as µ, ν, σ and Cartesian indices α, β. In order to keep the notation as simple as possible we will use index-notation only if it facilitates readability. Hence, expressions such as BC, B µµ C µ ν and B µαµ α C µ α νβ all denote the same product of two matrices B and C, with contraction of the blob-index µ and the Cartesian index α . The different elements of the GLE (1) require the definition of an inner product of phase functions φ, ψ * . We use the classical mechanics limit of Mori's definition as an equilibrium average [12] with the equilibrium ensemble microscopic variables z, the partition function Z and the Hamiltonian H (z). The second definition in (2) was introduced for brevity of notation. It will be used most of the time, implicitly assuming the presence of the star which denotes the conjugate transpose in general and in particular the conventional transpose in all cases considered in this work.
In (1) the drift 0 , the memory function K 0 (t) and the random force F 0 (t) are given by where we use the short notation for initial values, e.g. A 0 ≡ A 0 (t = 0). L is the Liouville operator associated with the Hamiltonian H (z). With the inner product, the Liouville operator is assumed to be anti-Hermitian, i.e. Lφψ = − φ Lψ . Q is a projection operator projecting onto the space orthogonal to the space spanned by the relevant variables A 0 (z). It is defined by Q ≡ 1 − P where P is a projection operator projecting any phase space function ψ(z(t)) onto the space of relevant variables according to It should be noted that Mori theory is linear. One consequence is that all nonlinearities of the Hamiltonian H (z), leading to the coupling of eigenmodes, can only show up in the memory function K 0 (t) of the GLE (1). The linearity of the theory is essential for the manipulations presented in the following.

Mori hierarchy
By recursively repeating the procedure of splitting the dynamics into a relevant and an orthogonal part Mori constructed a recursive cascade of GLEs for projected forces [29]. The recursion relation readṡ for j 0 with the definitions The jth projection operator P j is defined by its action on a phase space function ψ(z) as in analogy to (8). This also gives P 0 = P. L is the conventional Liouville operator introduced earlier.

Continued fraction representation of the GLE
Transforming the hierarchy of GLEs (9) into the Laplace domain, one obtains the recursion relation [29] for the correlation function Ξ j (t) defined by with Laplace variable p and unit matrix 1. This can be rewritten as a continued fraction for any j, for example for j = 0 where we have expressed the continued fraction starting from level j + 1 by K j ( p) = Ξ j+1 ( p)∆ 2 j+1 according to (21). We can directly express the memory kernel in the same way, for example for K 1 (t) we have or, in an ambiguous but more graphic notation where the fractions must be interpreted as matrix products of the form used in (25) with prepended inverses (i.e. denominators). K j ( p) represents a kernel at level j and can be expressed by another continued fraction of the same form starting with ∆ 2 j+1 in the nominator. Now it can immediately be seen that the Markovian approximation of (18) at level j amounts to set K j ( p) ≈ Γ j , i.e. to a constant friction matrix in the Laplace domain. This corresponds to the well known horizontal truncation of the continued fraction [36] and it is the essential approximation that we have to make. We do not believe that the dynamics on level j is generally Markovian. But, nonetheless, the representation of the essential object, which is the memory function K 0 ( p), improves, the more we go down in the continued fraction expansion, by including more and more additional static equilibrium properties ∆ 2 j and j as determined from the microscopic system. Furthermore, as shown in the following section, we choose Γ j in two ways that are both consistent with the continued fraction representation of K 0 ( p).

Approximations for the friction Γ j based on continued fractions
In this work we focus on the family of relevant variables A(z) = {X(x), P( p)} where the X(x) are linear combinations of the microscopic displacements x only and the P( p) are linear combinations of the microscopic momenta p only. Since x and p are statistically independent, so are then, due to the linear dependence, also X and P. The same is true for variables proportional to higher even and odd time-derivatives of X in general. This is equivalent to saying that odd time-derivatives of X(t)X vanish at t = 0, i.e. that all relevant correlation functions are even functions of time. This is the reason why, for the chosen family of relevant variables, the j vanish for j 1 as shown in appendix B.
This leads to simplified continued fractions for Ξ 0 ( p) and K 0 ( p) which read where we have already made the Markovian approximation (18), i.e. K j ( p) ≈ Γ j . Now, making a long-time approximation in (27), i.e. setting p = 0 and solving for Γ j leads to with Ξ 0 ≡ Ξ 0 ( p = 0). Note that which is the only time-integral we have to compute, regardless of the level j. For the first few friction matrices the explicit expressions are Making the same long-time approximation in (28) gives is the standard Markovian approximation for K j (t) for the special case j = 0 which contains projected dynamics and is hence not easily computable from MD in contrast to (30). We immediately see that the expressions (36) become identical to the expressions (29) for the approximation K 0 ≈ Ξ −1 0 + 0 . We can improve on this 'diffusive limit' by using the approximation F 1 Then the matrix resulting from the integral is to be inserted into (36), giving a second way to approximate the friction matrix. Both equations (30) and (38) are closures for the determination of Γ j , which are consistent with the Mori theory of continued fractions.

Equivalent system of linear equations
Besides an approximation for Γ j we need a convenient way to solve the truncated form of the GLE hierarchy (9) by time integration. For this purpose it is convenient to convert the integrodifferential GLE into a system of ordinary differential equations. Our starting point will be a reformulation introduced by Karasudani et al [33] which involves a change in perspective by considering a slightly different set of relevant variables A(z(t)). This set is with Note that, with this definition, contains the conventional Liouville operator L while the propagator of F j (t) contains the projected Liouville operator L j as defined in (16). Then, a new projection operator P = P is associated with this new extended set of relevant variables projecting now onto the subspace spanned by the extended vector A = A(t = 0). The projection operator is defined by The vector A represents an orthogonal set because it follows from the definition (12) With this orthogonality one can show that and obtains the GLĖ with the building blocks [33] = and with j , F j (t), K j (t) and ∆ 2 j given, respectively, by (10), (12), (21) and (23). Therefore, the GLE with these building blocks represents a linear system of equations of motion for A 0 (t) = {X(t), P(t)} and jth order fluxes A j (t) with j = 1..m, where only the mth order flux evolves according to a non-Markovian GLE. To this GLE and its memory kernel K m (t), both approximations can be applied which have been discussed previously in section 2.
is the very same mth memory kernel from the conventional Mori-hierarchy as described in section 2.2.
We approximate the constant friction matrix Γ m by (29), or by (36) (with (38)). The calculation of the first few ∆ 2 j can be found in appendix C.

Coarse-grained variables
We consider the atomic positions and momenta z ≡ {r i , p i } of i = 1..N atoms as microscopic variables. On the CG level we focus on relevant variables of the family A 0 (z) ≡ A µ0 (z) ≡ {X µ , P µ }. X µ is a generalized CG displacement and P µ is a generalized CG momentum defined by is the displacement of atom i from its equilibrium position r i0 and m i is its mass. The blob index µ runs over µ = 1..M, M < N so that f µi , f µi are rectangular matrices of constant weighting coefficients. Since we assume these matrices to behave like a diagonal matrix in the Cartesian indices α, β, we omit to write them in boldface if applied to a vector. Specific instances of these matrices are for example for COM variables [25] where the entries of the matrix θ µi take the value 1 if atom i is in blob µ and zero otherwise, and is the total mass of blob µ. One atom i may also contribute to more than one CG variable µ. This is for example the case in the CGMD method [27,28]. The most general requirements on In the special case of (49), f µi is just a scaled version of f µi , i.e.
That this relationship and (50) are also meaningful for all other physically useful linear combinations of microscopic variables (47, 48) follows from the additional requirement that the total COM of the microscopic and of the CG system should in general be consistent. Particularly COM-variables using (49) fulfil this requirement by construction. In general, if we write (47) as then multiplication with M µ , summation over µ and division by the total mass M tot gives exactly the desired relation for the displacements with and because of (50) and (51). For the total momentum of the system we sum (48) over all µ and get the desired relation µ P µ = i p i , again by using (51). Then (48) can be rewritten as where v i =ẋ i and (47) have been used. Here we have introduced the diagonal matrix allowing us to use the short matrix notation P = MV without indices if we wish.

Mori GLE for coarse-grained variables
By inserting those definitions we formally geṫ Equation (56) assures that the time derivative of the displacement is proportional to a relevant variable. Since F X 1 contains only the contributions from the irrelevant variables we have F X 1 = 0 and as a consequence of the fluctuation-dissipation relation (5) also K X X 0 (t) = K X P 0 (t) = K P X 0 (t) = 0. Further we get X X 0 = 0 and by using (47) and (48) and v i v j = k B T δ i j δ αβ /m i . Similarly we also get for the equation of motion for the CG momentum where we defined a generally non-diagonal mass matrix We also have P P 0 = L P P P P −1 = 0.
The only non-vanishing memory term reads for t t In summary the exact equations of motion arė where P X 0 is given by (63) and K P P 0 (t − t ) is given by (66).

Linear system of equations for the CG variables and auxiliary variables
Applying the scheme described in section 2.5 to our choice of CG variables, the whole system of equations (43) can be written in the forṁ The vector of relevant variables contains 2D(m + 1) entries, i.e.
The X-components of A 1 to A m vanish because the X-component of A 1 projects the time-derivativeẊ with Q 0 . But for the chosen family of CG variablesẊ is always proportional to P according to (67). P is another CG-variable, i.e. lying in the subspace spanned by A 0 . Then, any projection of P by a projector Q j is orthogonal to A 0 and must vanish. The matrix is given by (44) where each of the blocks j , ∆ 2 j , −1 or 0 is a square matrix of rank 2 × D × M, and there are (m + 1) × (m + 1) such blocks. We have The block P P j vanishes for our choice of A as shown in appendix C. The other blocks vanish due to the vanishing A X j . Note that the vanishing of j does not imply that the rows of the matrix in equation (44) For the first few levels j an explicit calculation in terms of covariances of CG variables and their derivatives at t = 0 gives (for details see appendix C) with F ≡Ṗ, and B 2 ≡ F F P P −1 + 2 2 P P .
The matrix Γ is the constant friction approximation of (46), i.e.
where we can obtain Γ m from one of the previously introduced approximations (29), or (36) (with (38)). Again, only the P P-component does not vanish and the random forceF(t) is given by (45) where only the D P-components of F m+1 (t) are non-zero. The non-trivial part of the Markovian fluctuation-dissipation relation is then where we can model the stochastic increment produced by the random force F P m+1 (t) as white noise, i.e. as a linear combination of Wiener processes such as where dW P m+1 (t) is a D × M vector of independent Wiener increments with the property dW P m+1 dW P m+1 = δ µν δ αβ dt ≡ 1 P P dt (83) and C P P m+1 is a constant D × M × D × M square matrix. Insertion of (82) into (81), time integration with Ito calculus, and using (83) gives which is a fluctuation-dissipation theorem relating the noise amplitudes C P P m+1 to the friction matrix Γ m . The second equality follows from the recursive relation (81) and where the mass matrix M was defined in (64). If C P P m+1 is symmetric then C P P m+1 C P P m+1 T = C P P m+1 C P P m+1 and (84) requires the computation of a matrix square root. In summary the equations of motion can be simplified and written in the form of stochastic differential equations as follows:

Application to centre of mass blobs in graphene
We particularise the general equation (86) and consider here as an example the in-plane 2D motion of the graphene honeycomb lattice (cf figure 1). Before writing down the specific instances of the matrices in (86) and computing some of them by MD simulation we first have to define the atomistic model.

Microscopic model
The dynamics of the graphene model is described by Hamilton's equations of motion for the set of microscopic variables {r i , p i }, where m i is the mass, r i the position and p i the momentum of C-atom i. For the sake of simplicity, instead of more realistic potentials such as the Brenner bond order potential [37], we use a simple 2D potential of the same form as already used in [27]. The interaction potential between C-atoms is assumed to have the form U = U s + U a . The pair potential U s models a harmonic spring and is defined as where k s > 0 is a stiffness constant, r 0 is the equilibrium bond-length, r i j = |r i j |, e i j = r i j /r i j , r i j = r i − r j and P i j denotes all those and only those pairs of nearest-neighbour C-atoms i and j that are bonded and assumed to interact by the given potential. The angular potential U (α i jk ) is defined as where T i jk denotes all those and only those triplets of neighbouring C-atoms {i, j, k} that are connected as shown in figure 1 and assumed to interact by the angular potential. The equilibrium angle is α 0 = 2π/3 and cos(α i jk ) = r i j · r k j r i j r k j . (89) Here, k a is a second stiffness constant for angular changes. The simulation units and parameters of the microscopic model are chosen in the same way as in [27]. The unit of length l * is fixed by setting l * ≡ r 0 ≈ 0.142 nm in graphene. The unit of mass is m * ≡ m C where m C ≈ 2 × 10 −26 kg is the mass of one C-atom. The unit of time is chosen as t * ≡ (m C /k s ) 1/2 , which amounts to set the value of k s = 1 in the selected units. The angular spring constant k a = 0.12 is found by requiring the model to reproduce the ratio of transverse to longitudinal in-plane speeds of sound c T /c L ≈ 0.633 based on measured values from [38,39]. We then get t * ≈ 5.18 fs. The unit of temperature is T * ≡ m C r 2 0 /k B t * 2 ≈ 1.09 × 10 6 K. MD-simulations are performed in the (N, V, E)-ensemble at a temperature of T = 2.74 × 10 −4 (corresponding to T ≈ 298K). Typically, the system is equilibrated for a sufficiently long time, of the order of 10t * and data are collected subsequently. The used time step in all MDsimulations is t = 0.005t * .

Coefficient matrices
In the following we list the specific forms of the matrices 2 j P P for j = 1, 2, 3 for COMblobs of equal mass M cg . They are obtained by noticing that the momentum covariance becomes diagonal, i.e. P µ P ν = M cg k B T δ µν . Then, also the mass matrix M µν is diagonal and identical to M µν , i.e. M µν = M µν = M cg δ µν . These replacements have to be performed in the general results given earlier and in the appendix and supplementary material (available from stacks.iop.org/NJP/15/125015/mmedia) where the detailed calculations can be found. Then we obtain 2 P P x-axis, xx-entry x-axis, yy-entry y-axis, xx-entry y-axis, yy-entry 0.
x-axis, xx-entry x-axis, yy-entry y-axis, xx-entry y-axis, yy-entry with Similarly for each of the two approximations for the friction matrices Γ P P j , the COM-blob expressions can be found, which also require the 2 j P P given above.

Numerical computation of some coefficient matrices
We show numerical results for the P P-components of ∆ 2 1 and ∆ 2 2 for a system of 10368 C-atoms that we coarse-grain into 108 hexagonal COM-blobs with N cg = 96 C-atoms and 432 COM-blobs with N cg = 24 C-atoms (cf figure 1). ∆ 2 1 is a crucial part of all approximations of the non-Markovian hierarchy (86) with m > 0. It couples the momentum of blob µ to the first auxiliary variable of blob ν. ∆ 2 2 couples the first to the second auxiliary variable. Therefore ∆ 2 1,2 are good representatives of the characteristic properties of the approximated memory in the equations of motion of the COM-blobs. A full analysis of the CG equations of motion will be presented in a forthcoming paper. Figure 2 shows plots of the Cartesian x xand yy-entries of the matrices ∆ 2 µν1 and ∆ 2 µν2 for pairs of blobs µ and ν aligned along the xor y-axis versus their distance R µν , for N cg = 96 (left) and for N cg = 24 (right). At a coarse-graining ratio of N cg = 96, the resulting periodically arranged 108 blobs have neighbours at four different non-zero distances on the x-axis and three on the y-axis. In the system of 432 blobs with N cg = 24, we have neighbours at nine different non-zero distances on the x-axis and six on the y-axis. It can be observed that the coupling of the momentum P µ of a blob µ to the auxiliary fluxes A ν1 (t) (by ∆ 2 1 ) and of A µ1 (t) to A ν2 (t) (by ∆ 2 2 ) is very local, i.e. the coupling to A ν=µ, j (t) is very dominant and only the absolute values of the coupling coefficients to the nearest blobneighbour are significantly larger than zero. One difference between ∆ 2 µν1 and ∆ 2 µν2 is the sign of the yy-entry of the nearest neighbour coupling (along the x-axis). At the moment we can not provide any explanation for this observation, but, as can be seen from the plots, this change in sign occurs consistently for both N cg = 96 and 24. In general, neighbours in y-direction do not seem to influence the blobs significantly. We attribute this behaviour to the nature of the lattice of the hexagonal blobs, which are directly aligned in the xbut not in the y-direction (cf figure 1). The local behaviour suggests a nearest neighbour approximation in order to gain computational efficiency when computing the equations of motion (86) numerically.

Summary and discussion
We obtained Markovian equations of motion for a general many body system of interacting CG displacements X(t), momenta P(t) and additional fluxes A j (t). As CG variables we considered the family of linear combinations of atomistic degrees of freedom. The goal was to coarse-grain the system as a whole instead of picking out individual variables of interest such as distinguished atoms or modes. The system of Markovian equations of motion approximates Mori's exact non-Markovian GLE and is easier to solve by computer simulation since it is linear and memoryless. All parameters of the equations can be obtained from equilibrium MD simulations of the investigated microscopic system. These parameters are the static covariances ∆ 2 j known from Mori's continued fraction and generalized constant friction matrices Γ j . We proposed two different ways to compute the friction matrices, which require the computation of relatively few static covariances (mostly the ∆ 2 j up to the desired level j = m) and only one and the same time integral, independently of the level j = m. A detailed comparison of these approaches will be presented in a forthcoming paper.
Additionally, suitable criteria for the truncation of the hierarchy must be numerically investigated. On a general ground, the most obvious criterion is whether and how fast timecorrelation functions computed from the CG model, such as the CG momentum correlation function, converge to their counterparts computed from the microscopic dynamics. More specifically, we can subdivide the behaviour of time correlation functions into short time behaviour and long time behaviour. For the former, by including a certain number of ∆ 2 j parameters, the presented approach allows a direct control of the number of frequency moments (or 'sum rules' [36]) that should be reproduced by the CG model. The essential criterion with respect to the long time behaviour is to test the reproduction of decay rates of eigenmodes, i.e. of macroscopic transport coefficients.
In this work, we presented numerical results for the special case of COM variables in the graphene lattice that indicate that the CG variables interact with their additional fluxes in a spatially very local way, following the local definition of the COM-variables themselves. This property is useful for efficient numerical computation. We also know from previous work that, on one hand, the time evolution of COM-variables seems to be more Markovian than for example for finite element (FE) based variables [25,27]. On the other hand, COM-variables show weaknesses in reproducing dispersion relations. Therefore we are also investigating the application of the presented general non-Markovian coarse-graining scheme to the special case of FE variables which also represent linear combinations of atomistic degrees of freedom. Therefore, the general results presented here allow one to address a huge variety of molecular systems such as atomic lattices, proteins, polymers or fluids, with one or more of the preferred linear CG variables in that field.
A true separation of timescales at level j, i.e. the strict validity of a Markovian approximation at this level, can be expected for example for lattices with a phononic bandgap. In this case we might wish to represent the bands below the bandgap in a CG model. Either we directly choose the respective normal modes as the relevant variables, which could lead to a purely Markovian model, or alternatively to these global collective variables, we might stick to a smaller number of spatially local COM-blobs with auxiliary variables as internal degrees of freedom, and search for a level j which separates the fast from the slow dynamics. It should be kept in mind that the presence of nonlinearities in the lattice is mandatory. Otherwise, the modes do not couple and there is no relaxation to equilibrium, which renders Mori-Zwanzig theories non-applicable. Nonlinearity is thereby necessary but not sufficient [40]. Other systems where a separation of timescales might be identified, are large proteins coarse-grained into a few still large substructures. The auxiliary variables could again represent relevant internal degrees of freedom of these substructures. But even for systems without any separation of time-scales, we may expect that higher approximation levels j reproduce more features of the memory function, and hence give a better approximation to the non-Markovian behaviour.
A further issue to be investigated is the requirements on the accuracy of time-integration. The described methods require the calculation of correlations of higher-order time derivatives of the relevant variables, which translates into the computation of higher-order time derivatives of the microscopic variables. This may pose a problem for standard second order accurate integrators, e.g. of Verlet-type, as we have used up to now. For the time being we have partially circumvented this potential numerical problem by calculating the required derivatives analytically, which was possible for the simple microscopic interaction potential used here. But this is of course not an option for general potentials and for arbitrary truncation levels. In addition it does not improve the time-integration itself. Therefore it is worth while to compare the results with those obtained with higher order symplectic integrators (cf e.g. [41]).
An alternative to Mori's projector based on equilibrium averages is Zwanzig's general nonlinear theory where the projectors are defined by constrained averages. The drift term or the memory function are then objects depending on the whole phase space of relevant variables. For practical computations, this makes modelling assumptions indispensable. In the work by Hijón et al [42], it was reasonable to assume pairwise interactions for a molecular liquid. This can not be directly transferred to lattices that are the main focus of the present work. In addition, the connection to the Mori-hierarchy is not obvious. A promising starting point that begins from the Mori hierarchy and goes a step into Zwanzig's direction seems to be the work of Dygas et al [43] and of Grigolini [44]. Dygas et al [43] use a GLE where they multiply the memory kernel with a state-dependent factor. The authors show a Markovian system of equations with auxiliary variables derived from this extended GLE. Their approach is based on Grigolini's work [44] who derived a GLE for non-additive fluctuations which allows one to describe a multiplicative process. It is worth investigating whether this approach can lead us to a useful nonlinear Markovian system of equations with auxiliary variables.
Some of the equations we refer to in this appendix can be found in appendix D which we provide as supplementary material (available from stacks.iop.org/NJP/15/125015/mmedia). From [29] we know since P j P k = 0 for all j = k. The first few projections P 0 , P 1 and P 2 are defined by their action on a phase space function G µ according to with A j = F j . By this we obtain The last equality is shown in (C.5). For A 2 = F 2 we get since all j vanish for j 1 as shown in appendix B. We also need with We already took into account that all j vanish and used (C.9) and (D.81).
For A 3 = F 3 we get where j = 0 for j 1 was used. We further need 12) where the vanishing of γ 1 and δ 0 is shown in (D.85) and (D.90), respectively, and β 2 is given in (D.71).

Appendix B. The matrices 1 , 2 and 3 vanish
Some of the equations we refer to in this appendix can be found in appendix D which we provide as supplementary material (available from stacks.iop.org/NJP/15/125015/mmedia). To keep the notation simple we will omit the index 0 for the CG variables from here on and write A(t) ≡ A 0 (t) = {X(t), P(t)} and L A(t) = {Ẋ(t),Ṗ(t)}. We will also write ≡ 0 . There is where is calculated in detail in (D.12) and L F 1 F 1 vanishes as shown in (D.14). Then, 1 = 0 under the following interpretation: in fact F 1 F 1 = 0, hence it is not invertible and 1 not computable. Recalling that we have already shown in (66) that ultimately the only relevant nonzero part of the memory kernel is the P P-component, we can interpret the inverses of the form F j F j −1 as a pseudoinverse, resulting in inverses of zero X X-subblocks to be set to zero. In this way, 1 is well defined and indeed 1 = 0. This procedure will therefore be applied here and for similar cases in the following. One such case is also 2 where F 2 F 2 has been calculated in (D.28) and L F 2 F 2 vanishes as shown in (D.30). This leads to 2 = 0 if F 2 F 2 −1 is computed as a pseudoinverse as discussed above for F 1 F 1 −1 . Further where F 3 F 3 has been calculated in (D.56) and L F 3 F 3 vanishes as shown in (D.70). This leads to 3 = 0 if F 3 F 3 −1 is computed as a pseudoinverse as discussed above for F 1 F 1 −1 .
Appendix C. Calculation of the ∆ 2 j Some of the equations we refer to in this appendix can be found in appendix D which we provide as supplementary material (available from stacks.iop.org/NJP/15/125015/mmedia).

C.1. Calculation of ∆ 2 1
Mori's definition is given in (D.2) and F 1 F 1 calculated in (D.12). Inserted into (C.1) this gives and F ≡Ṗ. The following relation is also worth noting: with given in (D.93).

C.2. Calculation of ∆ 2 2
Mori's definition is where F 1 F 1 has been calculated in (D.12) and F 2 F 2 in (D.28). Then we get It can also be shown that Mori's definition is where F 2 F 2 is given by (D.28) and F 3 F 3 in (D.56), which is