Response to “Comment on ‘Markovian approximation in a coarse-grained description of atomic systems’ ” [J. Chem. Phys. 128, 147101 (2008)]

We show within the Mori theory of projection operators that the Green-Kubo formula using microscopic rather than projected forces is only valid if the description is Markovian. Therefore, the only way to assess whether a description is Markovian is through examining the predictions of the theory under the Markovian assumption. Although, in principle, the blob description for a unidimensional chain is non-Markovian, in practice the Markovian approximation describes reasonably well the coarse dynamics of the chain.

In a series of papers^1–5 a debate about the validity of the Markovian assumption in a coarse-grained blob description of one-dimensional interacting particles has arisen. Cubero has detected a mathematical error in Ref. 5 that led him to conclude that our claim that a Markovian description is valid for a high temperature Lennard-Jones fluid is not correct. We have to admit that non-Markovian effects cannot be discarded in this model, and many of the observations made by Cubero are correct. However, although, in principle, the system may exhibit memory effects, in practice these effects may be small.

The crucial point of the discussion is whether the correlation of projected forces, i.e., the kernel, may be approximated with the correlation of microscopic forces [i.e., the corrected Eq. (3) of Cubero]. Rather than focusing the discussion in the Laplace transformed Eq. (3), we place directly the discussion of Markovian effects on the general projection operator theory of Mori and emphasize two aspects of it: (i) Being a linear theory, Markovian behavior is strictly linked to an exponential evolution of relevant variables (in general a complex exponential matrix) and (ii) the usual approximation of the projected forces with real forces in Green-Kubo formulas is strictly justified only when the behavior is Markovian and the decay is exponential. While the first point is well-known, the second point seems not to be sufficiently acknowledged in the literature and it is of relevance in the present discussion. Let us see in detail these points.

(i) The Mori theory of projection operators allows one to obtain the following exact linear equation for the matrix of equilibrium correlation functions $Aij=⟨AiAj(t)⟩$⁠, of relevant variables $Ai$

where $Ωij=⟨AjAk⟩−1⟨AkiLAi⟩$⁠, $Kij(t)=⟨AjAk⟩−1C̃ki$⁠, $F̃i(t)=Qexp{iLQt}iLAi$⁠. The projected forces $F̃i(t)$ satisfy $⟨F̃i(t)⟩=0$⁠, $⟨F̃k(0)F̃i(t)⟩=C̃ki(t)$⁠, $⟨AkF̃i(t)⟩=0$⁠, where ⟨⋯⟩ denotes an equilibrium average.

If the time scale of decay of the memory kernel is much shorter than the typical time scale of the relevant variables, then one is entitled to approximate the memory kernel as $Kij(t)=γijδ(t)$ where the friction matrix is given by a Green-Kubo formula,

If this separation of time scales exists, then the integrodifferential Eq. (1) becomes an ordinary differential equation,

which predicts an exponential matrix decay of the relevant variables $A(t)=exp{(Ω−γ)t}A(0)$⁠. This is an essential feature of the linear Mori theory.

(ii) Next we proof the following lemma: If the dynamics is Markovian [i.e., Eq. (3) holds and the macroscopic variables decay exponentially], and the microscopic force correlation matrix has a rapidly decaying part plus (small) terms that evolve in macroscopic time scales, then it is possible to calculate the friction matrix through the following Kirkwood expression:

where, as opposed to Eq. (2), the correlation of “microscopic forces” $Fi=iLAi$ appears and the upper limit of the integral is not infinite but a time $τ$ which is long with respect to the initial decay of the microscopic forces but small compared to the time scale of macroscopic variables. The proof is as follows. We introduce the matrix of correlation of microscopic forces $C(t)$⁠, defined by

This is Eq. (7) of Cubero that simply follows from the adjoint property of the Liouville operator. The correlation of microscopic forces $C(t)$ can be obtained from a molecular simulation, while the correlation of projected forces $C̃(t)$ cannot. The relationship between $C(t)$ and $C̃(t)$ is easily obtained by taking the time derivative of Eq. (1) and use of Eq. (5),

Note that $C(t)$ and $C̃(t)$ differ by the last two terms of the right hand side of Eq. (6). These two terms have been neglected in Ref. 5 and this has been questioned by Cubero in his Comment.

We will assume that the very short range decaying part of $C(t)$ may be modeled with a delta function. The time scale of this decaying part is much shorter than the time scale of the relevant variables. Obviously, $C(t)$ should have another contribution for two different reasons: First, it is required that $∫0∞C(t)dt=0$ (see Ref. 4) and, second, Eq. (5) implies that the time scale of the microscopic force is the same as the time scale of the relevant variables, as noted by Cubero. We can model $C(t)$ by assuming that the correlation $A(t)$ of relevant variables has the following structure, $A(t)=a(t)θ(t)+a(−t)θ(−t)$ where $a(t)$ is a matrix function with the property $ȧ(0)≠0$ and $θ(t)$ is the Heaviside step function. With this, we are assuming that $A(t)$ has a cusp at the origin. From Eq. (5) the correlation of microscopic forces is obtained from the second time derivative of $A(t)$⁠; this is

The first contribution is singular at the origin while the second contribution has the time scale of the relevant variables. Of course, in reality we never have a cusp in $A(t)$ and, consequently, never have a delta function in $C(t)$⁠. We rather have a tiny rounded cusp in $A(t)$ and a highly peaked decaying contribution to $C(t)$⁠, but the singular modelization serves for our purposes.

It is apparent from Eq. (6) that if the correlation of microscopic forces has a delta contribution as in Eq. (7), then $C̃(t)$ should have also the following structure:

where $ϵ(t)$ is a nonsingular contribution to the correlation which we expect to have a time scale comparable with the time scale of relevant variables. The substitution of Eq. (8) into Eq. (1) leads to the following equation for the nonsingular contribution $ϵ(t)$⁠:

Observe that if the Markovian property holds, then $a(t)$ obeys Eq. (3) and $a(t)$ decays exponentially. As a consequence, Eq. (9) predicts $ϵ(t)=0$⁠. We have thus demonstrated that $C̃(t)=2ȧ(0)δ(t)$ in the Markovian limit, and any memory contribution in the matrix of correlation of projected forces is consistently eliminated, even though $C(t)$ do has a memory. Because $C̃(t)=2ȧ(0)δ(t)$⁠, the Green-Kubo relation in Eq. (2) can be computed through Eq. (4). Note that the argument presented shows that the usual way of computing the friction matrix with the microscopic rather than the projected forces is, in strict sense, only valid if the macroscopic variables evolve exponentially and the Markovian approximation is valid. An explicit example of this situation is Brownian motion in the infinite mass limit.⁴ 

We are now in position to discuss the validity of the Markovian approximation for our problem of the one-dimensional chain of interacting particles. It is very difficult to say whether the macroscopic variables, like the momentum of the blobs, follow or not an exponential behavior, because the time dependence of an element of an exponential matrix may be very convoluted. Nevertheless, in the harmonic case the time dependence of the momentum is a piece-wise linear function and for the Lennard-Jones system at moderate temperatures it retains a similar shape. Strictly speaking, therefore, we would not be allowed to neglect the term $ϵ(t)$ in the projected force correlation. In fact, it is possible to show that the structure of this term is of the form $1∕nϵ*(t∕n)$⁠, where the functional form of $ϵ*(x)$ is independent of $n$⁠. Due to the fact that the amplitude is vanishingly small, we assumed in Ref. 5 that this term may be neglected. Rigorously, it cannot be neglected because this term appears within an integral and has a time scale which is comparable to that of the relevant variables. Cubero has correctly understood this point. Of course it is very difficult to asses quantitatively the importance of the term $ϵ(t)$⁠, as this amounts to solve the integrodifferential Eq. (9), which is not easy.

In order to quantitatively address this problem, we proceed through a direct route by assuming that a Markovian approximation can be taken. From Eq. (9) of Ref. 5, the Markovian approximation is

where the “friction matrix” $γμν$ is defined in terms of a Green-Kubo formula [Eq. (4)] which involves the microscopic forces and can be easily computed from the simulations. The results indicate that the friction matrix can be approximated by $γμν≈γAμν$⁠, where $Aμν$ is the Rouse matrix and $γ$ is the time integral of the microscopic force autocorrelation function evaluated at the plateau divided by $kBT$⁠. In Fig. 1 we show the autocorrelation function of the microscopic force, displaying the plateau that allows us to obtain the friction coefficient $γ=4.35(mϵ∕σ2)1∕2$⁠. The matrix $Xμν=⟨XμXν⟩$ has been computed from the simulations and numerically inverted. The resulting matrix has the following structure $Xμν−1≈kLJ(T)∕nkBTAμν$⁠, where $KLJ(T)$ is an effective spring constant for the Lennard-Jones system that depends slightly on the temperature but not on cluster size. For$kBT=0.6ϵ$ it takes the value $kLJ=142ϵ∕σ2$⁠.

By multiplying both sides of Eq. (10) by $Pσ$⁠, averaging, and taking a time derivative, one arrives at a dynamic equation for $C¯μσ(t)=⟨Pμ(t)Pσ(0)⟩$⁠,

The solution of Eq. (11) with initial conditions $C¯μν(0)=nmkBTδμν$⁠, $C¯̇μν(0)=0$ is

In this equation, we have introduced the following parameters $γμ≡γ∕nm(aμ)$⁠, $ω̃μ≡ωμ(1−γμ2∕4ωμ2)1∕2$⁠, and $ωμ≡[kLJ∕mn2(aμ)]1∕2$⁠. The orthonormal matrix $Λ$ diagonalizes the Rouse matrix and $aμ$ are the eigenvalues of the Rouse matrix. Their explicit expressions may be found in Eqs. (5) and (7) of Ref. 5. In Fig. 2 the comparison of the Markovian prediction in Eq. (12) with the simulation results for a cluster size of $n=10$ shows a very reasonable agreement, provided that there are no adjustable parameters (⁠$γ$ and $kLJ$ are extracted from the molecular dynamics simulation).

In summary, we have shown that the Green-Kubo formula using microscopic rather than projected forces can only be used if the description is Markovian. From a practical point of view this severely limits the usefulness of Mori theory to Markovian descriptions, otherwise it is not possible to extract transport coefficients from molecular dynamics. In addition, the only way to assess whether the Markovian approximation is a good one is through examining the consequences of the Markovian assumption, because the memory kernel is not accessible from simulations if the description is not Markovian. We have shown that even though, in principle, the arguments of Cubero hold in the coarse-grained blob system for a one-dimensional atomic system, the behavior of the macroscopic variables is rather well described under the Markovian assumption. In that respect, even though the time scale of the relevant variables is governed by sound propagation, a Markovian description is able to reproduce the appropriate time scales of the relevant variables.

The authors acknowledge financial support from the Ministerio de Educación y Ciencia, FIS2007–65869–C03–03. C.H. and M.S. thank “Programa Propio de la Investigación de la UNED” for financial support.