Large fluctuations in non-ideal coarse-grained systems

TL;DR

This paper investigates fluctuations in non-ideal systems using the MDLG approach, revealing that fluctuations do not diminish with coarse-graining as in ideal gases, impacting the accuracy of current hydrodynamics simulations.

Contribution

It demonstrates the limitations of ideal gas assumptions in coarse-grained systems and highlights the need to modify existing fluctuating hydrodynamics models for non-ideal systems.

Findings

Fluctuations persist in non-ideal systems upon coarse-graining.

MDLG provides insights into the nature of fluctuations in such systems.

Current models may require augmentation for accurate non-ideal system simulation.

Abstract

Using the recently introduced Molecular Dynamics Lattice Gas (MDLG) approach, we test fluctuations of coarse-grained quantities. We show that as soon as the system can no longer be considered an ideal gas fluctuations fail to diminish upon coarse-graining as is usually expected. These results suggest that current approaches to simulating fluctuating hydrodynamics may have to be augmented to achieve quantitative results for systems with a non-ideal equation of state. The MDLG method gives a guidance to the exact nature of the fluctuation in such systems.

Equations

$$n_i(x+v_i,t+\Delta t)=n_i(x,t)+{\Xi }_i$$

$$P(n_i)=\exp (-f_i^{eq})(f_i^{eq}{)}^{n_i}/n_i!$$

$$n_i(x,t)=\sum _{j=1}^N{\Delta }_x[x_j(t)]{\Delta }_{x-v_i}[x_j(t-\Delta t)]$$

$$a^2=\frac{⟨\delta (\Delta t{)}^2⟩}{d\Delta x^2}\approx 0.18\approx \frac{1}{6}$$

$$⟨(\rho -{\rho }^{eq}{)}^2⟩=\sum _i\sum _j⟨(n_i-f_i^{eq})(n_j-f_j^{eq})⟩.$$

$$⟨n_i⟩$$

$$=$$

$$\sum _{k=1}^N{\Delta }_x(x_k){\Delta }_{x-v_i}(x_k-{\delta }_k)$$

$$=$$

$${\Delta }_x(x_1){\Delta }_{x-v_i}(x_1-{\delta }_1)$$

$$=$$

$$P(n_i)$$

$$=$$

$$\Theta (n_i;\{x_1,{\delta }_1,\cdots ,x_N,{\delta }_N\})$$

$$=$$

$$\displaystyle\left(\begin{array}[]{c}N\\ n_{i}\end{array}\right)\Theta(n_{i};\{x_{1},\delta_{1},\cdots,x_{n_{1}},\delta_{n_{1}}\})$$

$$\Theta (n_i;\{x_1,{\delta }_1,\cdots ,x_M,{\delta }_M\})$$

$$=$$

$$P^{N,id}(x_1,{\delta }_1,\cdots ,{\delta }_N)=\prod _{k=1}^N{P}^1(x_k,{\delta }_k)$$

$$P^{id}(n_i)$$

$$\approx \exp (-f_i^{eq})\frac{(f_i^{eq}{)}^{n_i}}{n_i!}$$

$$⟨n_i(x)n_j(y)⟩$$

$$=$$

$${\Delta }_x(x_k){\Delta }_{x-v_i}(x_k-{\delta }_k){\Delta }_y(x_l){\Delta }_{y-v_j}(x_l-{\delta }_l)$$

$$=$$

$${\Delta }_x(x_1){\Delta }_{x-v_i}(x_1-{\delta }_1){\Delta }_y(x_2){\Delta }_{y-v_j}(x_2-{\delta }_2)$$

$$+N\int d{x}_1\int d{\delta }_1{P}^1(x_1,{\delta }_1)$$

$${\Delta }_x(x_1){\Delta }_{x-v_i}(x_1-{\delta }_1){\Delta }_y(x_1){\Delta }_{y-v_i}(x_1-{\delta }_1)$$

$$=$$

$${\Delta }_x(x_1){\Delta }_{x-v_i}(x_1-{\delta }_1){\Delta }_y(x_2){\Delta }_{y-v_j}(x_2-{\delta }_2)$$

$$+f_i^{eq}{\delta }_{ij}{\delta }_{xy}$$

$$⟨n_i(x)n_j(y){⟩}^{id}=(1-1/N)f_i^{eq}{f}_j^{eq}+f_i^{eq}{\delta }_{ij}{\delta }_{xy}$$

$$W_{ii}=\frac{⟨(n_i(x)-f_i^{eq}{)}^2⟩}{{\rho }^{eq}{f}_i^{eq}}.$$

$$P^2(x_1,v_1,x_2,v_2)\propto P(v_1)P(v_2)g(|x_1-x_2|)$$

$${\delta }_k={\int }_0^{\Delta t}{v}_k(t)dt$$

$$\underset{|x_2-x_1|\to \infty }{\lim }{P}^2(x_1,{\delta }_1,x_2,{\delta }_2)=P^1(x_1,{\delta }_1)P^1(x_2,{\delta }_2)$$

$$\frac{⟨{\delta }_1{\delta }_2⟩(|x_1-x_2|)}{⟨{\delta }_1^2⟩}\approx 2(\Delta x{)}^{-1/2}{\rho }^{1/2}\exp \left(-\frac{|x_1-x_2|}{\xi \Delta x}\right)$$

$$P$$

$$\times \exp \left(-\frac{({\delta }_1+{\delta }_2{)}^2}{4{\sigma }_{+}^2(r)}\right)\exp \left(-\frac{({\delta }_1-{\delta }_2{)}^2}{4{\sigma }_{-}^2(r)}\right)$$

Full text

fluctuating hydrodynamics, lattice Boltzmann, Molecular Dynamics, kinetic theory

Predicting fluctuations at small scales for coarse-grained systems can be challengingHoyt et al. (2001); Kocher et al. (2016); Altaner and Vollmer (2012). Typically theoretical predictions exist for fluctuations in the hydrodynamics (i.e. long wavelength) limitGrosfils et al. (1992) that can be used to tune fluctuation terms, which is the basis of the Langevin approach to fluctuations Lemons and Gythiel (1997). Fluctuating terms in the governing equations are, however, usually assumed to be spatially uncorrelated, so that the hydrodynamic behavior has to emerge from the local fluctuations that are introduced into the system. To complicate the situation fluctuations are often most relevant at scales much smaller than the hydrodynamic length-scales. This makes it imperative to develop methods that are able to reproduce correct fluctuations at much smaller length-scalesKocher et al. (2016). It is, however, generally not known what form the small-scale fluctuations for such a coarse-grained system should take. This is a general problem that arises in most simulation contexts where fluctuations are important.

One obvious way to address this shortfall is to use information from a microscopic model, e.g. a Molecular Dynamics (MD) simulation, to measure the small scale fluctuationsHoyt et al. (2001). This information can then be used to tune the fluctuating terms in a coarse-grained simulation method Kocher et al. (2016). In this letter we show that it is possible to go a step further and use a coarse-graining mapping between a MD simulation and a coarse-grained mesoscopic model to directly measure the correct fluctuations of the model variables for different coarse-graining scales. Our results indicate that fluctuations of the model variables can be orders of magnitude larger than was previously expected.

In this letter we focus on predictions for fluctuating lattice Boltzmann approaches, although we would like to stress that the mapping approach is far more general and could be applied to most other coarse-grained models. For the special case of an ideal gas fluctuations are fairly well understood. This is why this is usually taken as a starting point, e.g for fluctuating lattice Boltzmann Ladd (1993); Adhikari et al. (2005); Dünweg et al. (2007); Kaehler and Wagner (2013). For other systems fluctuations arise from the discrete nature of the representation as in lattice gases Grosfils et al. (1992); Blommel and Wagner (2018) or Stochastic Rotation Dynamics Ihle and Kroll (2003); Tüzel et al. (2006). Other discrete versions like Dissipative Particle Dynamics Hoogerbrugge and Koelman (1992); Espanol (1995); Vázquez-Quesada et al. (2009) includes tunable fluctuating forces. For non-ideal systems, however, it is typically less clear what the correct fluctuations should look like Gross et al. (2010); Ollila et al. (2011); Thampi et al. (2011); Belardinelli et al. (2015). Because of this difficulty we developed a direct mapping from Molecular Dynamic (MD) onto a lattice gas (MDLG) Parsa and Wagner (2017) where fluctuations in a non-ideal coarse-grained system can be easily observed. In this letter we show the results of applying MDLG to analyzing equilibrium fluctuations at different densities of a system of Lenard-Jones particles.

The direct mapping between MD and a lattice gas, leads to an integer lattice gas. Integer lattice gases exist, but are somewhat rare Monaco (1988); Chen et al. (1997); Teixeira (1997); Boghosian et al. (1997); Chopard et al. (1998); Blommel and Wagner (2018). Blommel showed that an integer lattice gas can closely model the fluctuations of an ideal gas Blommel and Wagner (2018). Such lattice gases have occupation numbers $n_{i}(x,t)$ as their fundamental variables. They are defined on lattice points $x$ and are associated with lattice velocities $v_{i}$ such that $x+v_{i}$ is again a lattice position. The evolution equation of a lattice gas can be written as

[TABLE]

where the $\Xi_{i}$ is the lattice gas collision operator. We also define the number density $\rho=\sum_{i}n_{i}$. For the integer lattice gas of Blommel Blommel and Wagner (2018) the $n_{i}$ follow the expectation for an ideal gasLandau and Lifshitz (1969) and are Poisson distributed

[TABLE]

where $f_{i}^{eq}=\langle n_{i}\rangle$. The number density $\rho$ is similarly Poisson distributed with average $\rho^{eq}=\sum_{i}f_{i}^{eq}$.

To examine the fluctuation behavior of a non-ideal system (as represented by a Molecular Dynamics simulation) we employ the MDLG approach. Here a lattice is superimposed on the Molecular Dynamics simulation and the displacement of particles from one lattice cell to another a distance $v_{i}$ away is during a timestep $\Delta t$ is then identified with the lattice gas occupation number $n_{i}$. This provides us with a lattice gas that exactly represents the underlying Molecular Dynamics system Parsa and Wagner (2017). Explicitly, if $x_{j}(t)$ is the position of the $j$th particle at time $t$, we can write

[TABLE]

where $\Delta_{y}[x]$ is one if $x_{\alpha}\in[y_{\alpha},y_{\alpha}+\Delta x]$ and zero otherwise for all coordinates $\alpha$ and $N$ is the total number of particles.

We have previously shown that such an MDLG can have an equilibrium distribution that is equivalent to the standard lattice Boltzmann equilibrium distributions when the combination of the time step $\Delta t$ and lattice size $\Delta x$ is chosen such that

[TABLE]

where $\delta(\Delta t)$ represents the displacement of a MD particle in the simulation for a time-interval of $\Delta t$ and $d$ is the number of spatial dimensions Parsa et al. (2019). For this particular choice the equilibrium distribution of the lattice gas recovers the standard D2Q9 equilibrium distribution of lattice Boltzmann for small velocities Qian et al. (1992). While the MDLG approach will generate correct results for all values for $a^{2}$ the approach only recovers the standard lattice Boltzmann approach for $a^{2}\approx 1/6$. To ensure comparability of results to fluctuating lattice Boltzmann implementations Ladd (1993); Adhikari et al. (2005); Gross et al. (2010); Ollila et al. (2011) we restrict ourselves to this value here.

We now examine the equilibrium fluctuations of this lattice gas. We start with a dilute gas, for which we expect the assumptions of ideal gas fluctuations to hold to good approximation. We perform a Molecular Dynamics simulation with $99856$ particles on a $3000\times 3000$ lattice (length measured in LJ units) at a high temperature of 50, also in LJ units as in Parsa and Wagner (2017), to prevent phase-separation. If we take the effective radius for excluded volume at about $r_{c}=0.75$ [corresponding to $1/2k_{B}T=V(r_{c})$], this translates into a volume fraction of $\phi=4.9\;10^{-3}$. For this volume fraction we indeed find good agreement between the distribution of the lattice gas occupation numbers and the Poisson distribution, as is shown in Figure 1 for two different lattice spacings $\Delta x$.

For larger densities, where the excluded volume becomes important, the total density will no longer be Poisson distributed, but rather show a narrower distribution. We therefore would also expect the distribution of the lattice gas occupation numbers $n_{i}$ to similarly deviate from the Poisson distribution. We performed a similar simulation for a denser system of a nominal volume fraction of $\phi=0.49$. The results shown in Figure 2(a) show the narrowing of the $\rho$ distribution, but surprisingly the distribution of rest-particles ($n_{0}$) appears little changed where as the distributions of particle moving to nearest neighbors ($n_{1}$) as well as diagonal neighbors ($n_{5}$) is much wider.

We stated above that we chose $a^{2}=1/6$, but for each lattice spacing $\Delta x$ there is a time step $\Delta t$ that corresponds to this value. It is therefore essential how the distributions depend on $\Delta x$ (under the assumption $a^{2}=1/6$), corresponding to a modeling of the system at different scales. Increasing $\Delta x$ will increase the average number of particles per cell. For an ideal system the width of the Poisson distribution will grow only as the square root of the number of particles, making the distribution more peaked for larger number of particles. This means that for larger $\Delta x$ the importance of fluctuations declines. This classical result is found to good approximation for our simulation of a dilute system, as shown in Figure 1(b).

Now for a dense system we show the results for the larger $\Delta x$ in Figure 2(b). The density does indeed show a sharpening of the normalized width, as one would expect. The normalized width of the distributions of the $n_{i}$, however, does not show the same amount of narrowing. This is a counter intuitive result, which suggests that the importance of fluctuations for the $n_{i}$ does not diminish as $1/\sqrt{N}$ with increasing $\Delta x$.

To gain an understanding of this phenomenon let us first consider how combining wider distributions for the $n_{i}$ can lead to a narrower distribution of $\rho$. We can write

[TABLE]

For an ideal system the distributions of the $n_{i}$ are independent. For ideal systems $\langle(n_{i}-f_{i}^{eq})(n_{j}-f_{j}^{eq})\rangle=f_{i}^{eq}\delta_{ij}$. But for non-ideal systems the wider distributions for the $n_{i}$ requires that at least for some $i$ and $j$ these correlations must become quite negative to cancel the widening of the distributions for $i=j$.

To look at this more quantitatively we can examine the particle displacement probability for different particles $\delta_{j}=x_{j}(t)-x_{j}(t-\Delta t)$. First we can derive the expectation value

[TABLE]

where we introduced the useful (but not very common) N-particle distribution function for the displacements of particles during a time-step $\Delta t$. We see that the expectation value of the distribution is entirely dependent on the one-particle distribution function. This is the reason that Parsa et al. Parsa et al. (2019) found that the equilibrium distribution depends on the non-dimensional mean squared displacement $a^{2}$ only.

The probability for finding a specific occupation number $n_{i}$ is then given by

[TABLE]

where we define a binary flag that is one if we have the occupation number $n_{i}(x)$ as

[TABLE]

In the special case of an ideal gas, where the $N$-particle distribution factorizes

[TABLE]

we have

[TABLE]

where the last step is the familiar transition from the binomial distribution to the Poisson distribution in the limit of large $N$. We can shed light on the expected width of the distribution for the non-ideal case by reducing the expression to one for the two-particle distribution function:

[TABLE]

where, for an ideal gas, the two particle probability factorizes (13). With this we obtain for an ideal gas

[TABLE]

Which is the result that one would expect for an independently binomial distributed $n_{i}$. In the limit of a large system with $N\rightarrow\infty$ the $1/N$ can be neglected, and we obtain the result predicted for a Poisson distributed $n_{i}$.

For dilute gases the representation of the particles as ideal particles works quite well. We define the scaled squared width of the actual distribution as

[TABLE]

For a Poisson distribution this squared width will decrease as $1/\rho^{eq}$. We plot $\sqrt{W_{ii}}$ in Figure 3 for an ideal gas, a real dilute gas and a dense gas as a function of the number of average number particles per lattice site. This is accomplished by increasing the lattice spacing $\Delta x$, and increasing the time step $\Delta t$ at the same time such that $a^{2}=1/6$ as in Eq. (4). Compared to the result expected for a Poisson distribution the dilute gas agree quite well up to about a thousand particles per lattice site, but after that we see that the width appears to fail to decrease as fast as expected for an ideal gas. This is in agreement with the results for a dilute gas in Figure 1, which showed that the distributions agreed reasonably well with the prediction for a Poisson distribution.

For the dense gas the scaled width of the distribution diverges from the ideal gas case already at about 10 particles per lattice site, and decays much more slowly for larger numbers of particles. This implies that the importance of fluctuations does no longer decay as $1/\sqrt{\rho^{eq}}$. Therefore increasing the coarse-graining by choosing a larger lattice spacing $\Delta x$ will not diminish as rapidly as would be expected in standard Statistical Mechanics. The difference is striking: instead of decaying as $(\rho^{eq})^{-1/2}$ the rest-particles scale approximately as $(\rho^{eq})^{-1/4}$ and moving particle density only as $(\rho^{eq})^{-1/8}$. This result may seem counter intuitive: basic Statistical Mechanics would seem to demand that if we continue to double the lattice spacing $\Delta x$, the added components should eventually become independent. This argument, however, misses the important point that we are keeping $a^{2}=1/6$, so that the time step $\Delta t$ also increases.

To understand why the correlations do not vanish more quickly when scaling up the coarse-graining of the system let us consider the probability of the particle displacements. Firstly it is important to realize that instantaneous velocities in equilibrium remain uncorrelated

[TABLE]

where $g(r)$ is the pair correlation function. The same, however, is not true for displacement probabilities $P_{2}(x_{k},\delta_{k},x_{l},\delta_{l})$. Since these displacements can be written as

[TABLE]

the probability contains information about velocity cross correlations that are considered in the context of cross-diffusion constants. To understand this better we have to consider $P(x_{1},\delta_{1},x_{2},\delta_{2})$ of Eq.(18). For two-dimensional systems considered in this letter this probability depends on 8 variables. Because of overall translational and rotational invariance this could be reduced to 5 variables. However, this is still too large to be efficiently probed in our simulations. It is reasonable to expect that for particles sufficiently far apart the displacement probability will factorize:

[TABLE]

For particles closer together, however, the probability no longer needs to factorize. It is well known that, particularly for longer times the time-correlator for velocities decays only algebraically Dorfman et al. (1994). It is reasonable to expect a similar effect for velocities of particles that have a small spatial separation. A comprehensive numerical evaluation of these correlations is outside the scope of this letter, but we did examine $\langle\delta_{1}\delta_{2}\rangle$ as a function of the initial displacement of the particles $|x_{2}-x_{1}|$ Parsa (2018). We found that this correlator does appear to decay exponentially with a correlation length $\xi\approx\langle\delta^{2}\rangle$ with a prefactor that varies remarkably little with density (from 1 to 1.25 for the two extreme densities studied in this letter).

We now pick a length scale for our lattice $\Delta x$ and a corresponding time scale $\Delta t$ such that $a^{2}=1/6$. This means that $\xi$ measured in terms of $\Delta x$ is nearly constant. We can also determine the dependence of the amplitude on $\Delta x$ and the overall density to get the approximate relation

[TABLE]

where $\xi$ now varies slightly from 1 to 1.25 for our density range. There are two caveats to this approximate result: firstly at short distances the normalized correlation of Eq. (24) has to be less than one. Secondly the total expectation of $\langle\delta_{1}\delta_{2}\rangle=0$ because of momentum conservation. In numerical experiments we find that for large $\Delta x$ there is actually a negative correlation, but the value of this correlation function for large $\Delta x$ is not of interest here. Even with these caveats this result implies that the non-ideal part of the two-particle probability also decays exponentially with correlation length $\xi$. We can now make the following Ansatz that will recover both the factorization of Eq. (21) and the correlation of Eq. (24):

[TABLE]

with

[TABLE]

This is the simplest two particle distribution function that recovers both Eq. (23) and Eq. (24). We can use numerical integrations to estimate the occupation number correlators in Eq. (18). The numerical integration was performed using Mathematica, using the rough two particle correlation function

[TABLE]

and a notebook is included in the supplemental material. The results of this numerical integration is shown in Fig. 3, and it qualitatively recovers the result of the direct simulations.

In conclusion we have found that apparent equal-time fluctuations in coarse grained models contain time-correlations, which can significantly alter the scaling of the fluctuations. It is well known that while equal-time fluctuations in equilibrium systems decay exponentially, different-time fluctuations only decay algebraically Dorfman et al. (1994) so that the contribution of different-time correlations can be significant. This does not change the short-ranged nature of these fluctuations, as can be seen in Eq. (24) but the effects are nevertheless significant. We find that the magnitude of fluctuations can be orders of magnitude larger for the coarse-grained system variables when compared to quantities that are obtained from equal-time correlators like density or momentum fluctuations. This property of mesoscopic methods has not been fully appreciated to date. A long history of previous lattice gas approaches made a Markov approximation, which implies that occupation numbers were viewed as instantaneous quantities. Such a view is in contrast to our lattice gas which is a true coarse-graining of a reality (as represented by a MD simulation), and has fundamentally different properties. This interpretation of a lattice gas as a coarse-grained model represents a shift in perception, and we believe that our results will facilitate the development of more realistic coarse-grained fluctuating methods.

Lastly we would like to emphasize that the large fluctuations observed in this letter should not be thought of as restricted to lattice gas approaches, but that these will appear in all coarse-grained approaches, like those mentioned in the introduction.