Lattice Gas with Molecular Dynamics Collision Operator

TL;DR

This paper presents a lattice gas model derived from molecular dynamics that accurately captures equilibrium behavior and improves the understanding of temperature and pressure in lattice Boltzmann simulations.

Contribution

It introduces a lattice gas implementation based on MD that serves as an optimal model for simulating systems with lattice Boltzmann methods.

Findings

Equilibrium behavior of lattice Boltzmann aligns with the optimal lattice gas.

Temperature and pressure expressions relate to physical and lattice parameters.

A specific temporal discretization is required for accurate equilibrium in lattice Boltzmann.

Abstract

We introduce a lattice gas implementation that is based on coarse-graining a Molecular Dynamics (MD) simulation. Such a lattice gas is similar to standard lattice gases, but its collision operator is informed by an underlying MD simulation. This can be considered an optimal lattice gas implementation because it allows for the representation of any system that can be simulated with MD. We show here that equilibrium behavior of the popular lattice Boltzmann algorithm is consistent with this optimal lattice gas. This comparison allows us to make a more accurate identification of the expressions for temperature and pressure in lattice Boltzmann simulations which turn out to be related not only to the physical temperature and pressure but also to the lattice discretization. We show that for any spatial discretization we need to choose a particular temporal discretization to recover the lattice Boltzmann equilibrium.

Equations

$$n_i^{\ast }(x,t)=n_i(x,t)+{\Xi }_i(\{n_j(x,t)\})$$

$$\rho$$

$$\rho u_{\alpha }$$

$$\sum _i{\Xi }_i$$

$$\sum _i{v}_{i\alpha }{\Xi }_i$$

$$n_i(x+v_i,t+1)=n_i^{\ast }(x,t)$$

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

$$\frac{d{x}_i(t)}{dt}$$

$$\frac{d{v}_i(t)}{dt}$$

$$\Delta_{x}(x^{\prime})=\left\{\begin{array}[]{lcl}1&&x_{\alpha}<x^{\prime}_{i\alpha}(t)\leq x_{\alpha}+\Delta x\;\forall\alpha\in\{x,y,z\}\\ 0&&\mbox{otherwise}\end{array}\right.$$

$$n_i(x,t)=\sum _j{\Delta }_x(x_j(t)){\Delta }_{x-v_i}(x_j(t-\Delta t)).$$

$${\Xi }_i=n_i(x+v_i,t+1)-n_i(x,t).$$

$$N(x,t)=\sum _j{\Delta }_x(x_j(t)).$$

$$N(x,t)$$

$$=\sum _i\sum _j{\Delta }_x(x_j(t)){\Delta }_{x-v_i}(x_j(t-\Delta t))$$

$$=\sum _j{\Delta }_x(x_j(t)).$$

$$\sum _i{\Delta }_{x-v_i}(x_j(t-\Delta t))=1,$$

$$\sum _i{\Xi }_i$$

$$=$$

$$-\sum _k{\Delta }_x(x_k(t)){\Delta }_{x-v_i}(x_k(t-\Delta t))]\}$$

$$=$$

$$N(x,t)u(x,t)=\sum _i{n}_i(x,t)v_i$$

$$\sum _i{v}_{i\alpha }{\Xi }_i\ne 0\text{ in general}$$

$$V(x)=4ϵ\left[{\left(\frac{\sigma }{x}\right)}^{12}-{\left(\frac{\sigma }{x}\right)}^6\right].$$

$$\tau =\sqrt{\frac{m{\sigma }^2}{ϵ}}.$$

$$f_i^{eq}=$$

$$=$$

$$\delta x=x(t)-x(t-\Delta t).$$

$$f_i^{eq}=$$

$$=$$

$$W(x)=$$

$$=$$

$$P(v)=\frac{1}{(2\pi k_{B}T{)}^{d/2}}\exp \left(-\frac{(v-u{)}^2}{2k_{B}T}\right).$$

$$⟨(\delta x_{\alpha }{)}^2{⟩}^{bal}=k_{B}T(\Delta t{)}^2.$$

Full text

Lattice Boltzmann, Molecular Dynamics, kinetic theory

I Introduction

Lattice Boltzmann methods are an important computational tool that is most commonly employed to simulate hydrodynamic systemsSucci (2001), and it has been adapted to address many complex phenomena from turbulence Teixeira (1998); Yu et al. (2005) over multi-phase and multi-componet flowSwift et al. (1995); Shan and Doolen (1995); He and Luo (1997); Briant et al. (2004); Yuan and Schaefer (2006) to pore-scale simulations of porous mediaTölke (2002); Kang et al. (2006) and simulations of immersed boundaries Dupin et al. (2003); Jansen and Harting (2011). It derives its power from an underlying mesoscopic description that ensures exact mass and momentum conservation. The exact physical meaning of the lattice Boltzmann densities, however, remains poorly understood.

The lattice Boltzmann method was derived as a theoretical tool for the analysis of lattice gas methods McNamara and Zanetti (1988). Lattice gas methods consist of particles moving on a lattice with velocities that connect neighboring sites. After the particles have moved a stochastic collision step rearranges the particles. If these collisions conserve both the number of particles and the total momentum there will be a hydrodynamics limit for mass and momentum conservation equation. The introduction of a hexagonal instead of a square lattice by Frisch, Hasslacher and Pomeau Frisch et al. (1986) recovered the necessary isotropy to allow the momentum equation to be related to the Navier-Stokes equation.

These lattice gas models had some deficiencies, one unfavorable feature was a large and essentially uncontrolled amount of noise that required a significant amount of averaging. To derive the Navier-Stokes equations from the lattice gas dynamics a theoretical ensemble average was performed, leading to a lattice Boltzmann representation. Higuera then proposed to simulate the ensemble averaged lattice Boltzmann evolution equation directly, and thereby avoid the need to average results of the lattice gas equation Higuera and Jimenez (1989); Higuera et al. (1989). The collision operation of this first lattice Boltzmann method could be mapped one-to-one to the lattice gas and shared some of the positive features of the lattice gas, like the existence of an H-theorem with unconditional stability, and also some of its deficiencies like velocity dependent viscosities.

It was then realized that there existed much more freedom in the choice of the collision operator, and in particular the relaxation towards a local equilibrium function, often called the Bhatnagar-Gross-Krook (BGK) approach, allowed the full recovery of the Navier-Stokes equation to second order Qian et al. (1992).

At this time a second approach to derive the lattice Boltzmann equation directly from the continuous Boltzmann equation with a BGK collision operator gained popularity He and Luo (1997). Over the years several different local equilibrium distributions have been proposed, and currently the most popular method is a standard form of a second order expansion in terms of velocities.

Typically these lattice Boltzmann methods are validated by their ability to recover the Navier-Stokes equation. Here, however, we want to establish a relation to an underlying Molecular Dynamics simulation. For any Molecular Dynamics simulation, we can bin the particles into lattice cells corresponding to the lattice Boltzmann lattice. We can then observe where the particles in cell $x$ migrate to after a time $\Delta t$, and associate these particles with a lattice velocity vector $v_{i}=x(t)-x(t+\Delta t)$. These particles will collect at their new lattice cells. After another timestep $\Delta t$ these particles are re-distributed to new lattice sites, and can be associated with new lattice velocities. We call this representation of the MD simulation Molecular-Dynamics-lattice-gas (MDLG). This redistribution can be understood to be an effective MDLG-collision operator. In some very fundamental sense this is the collision operator that the lattice Boltzmann approach is trying to mimic. The purpose of this paper is to understand the physical meanting of the lattice Boltzmann densities in terms of this fundamental MDLG representation.

The paper is organized as follows: we first introduce a general idea of a lattice gas and then derive a new lattice gas which consists of a coarse-graining of an underlying MD simulation. We then apply this general idea to a specific MD simulation of a Lennard Jones gas in two dimensions. We analyse the equilibrium properties of the associated MDLG method and show that we are able to predict its mathematical form analytically. We then introduce the lattice Boltzmann method and compare the equilibrium properties of the MDLG method to the lattice Boltzmann equilibrium. We show that there are particular choices for the coarsgraining time and space discretization that lead to equilibria that are compatible with the lattice Boltzmann results.

II Lattice Gas

At its very basis a lattice gas consists of particles, all located on lattice points, that move with a lattice velocities $v_{i}$. What we mean by lattice velocity is that if $x$ is a lattice point, so it $x+v_{i}$. There are $n_{i}(x,t)$ particles at time $t$ at position $x$ moving with velocity $v_{i}$. The evolution consists of two steps. A collision step that redistributes the particles at the lattice point $x$ to different velocities:

[TABLE]

where the collision operator $\Xi_{i}$ is a function of all the particles and their velocities that are located at lattice point $x$ at time $t$. This collision operator will ensure that none of the locally conserved quantities are changed in the collision process. These locally conserved quantities will vary, depending on the desired physical system that one wants to model. In the majority of cases one will ensure mass and momentum conservation. Early lattice gases restricted the number of particles to at most one per velocity $v_{i}$ at a lattice site, and the velocity vectors all had the same length, ensuring that mass and energy conservation were synonymous. Each conserved quantity will lead to a corresponding hydrodynamic equation. Most applications focused on the fluid flow, and the key hydrodynamic equations to recover were the continuity and Navier-Stokes equations. Energy conservation is often abandoned in favor of an isothermal condition for many practical applications. Local mass and momentum densities are defined as

[TABLE]

and the conservation of these quantities then implies

[TABLE]

For the new kind of lattice gas collision operator proposed in section III we will see that mass conservation of Eq. (4) is indeed fulfilled, but the momentum conservation of Eq. (5) is not strictly obeyed. Because the new algorithm is based on MD, however, the algorithm conserves momentum rigorously. Its representation of momentum through Eq. (3), however, is inexact.

This collision is then followed by a streaming step

[TABLE]

where the particles move to the lattice site indicated by the velocity index $i$, i.e. they move from $x$ to $x+v_{i}$. The full evolution equation for these densities can then be written as

[TABLE]

Of course to make this description complete we need to define the collision operator. Originally lattice gases were defined such that there could be at most one particle for each $n_{i}$ Frisch et al. (1986). For the purpose of this paper, however, we will make no such restriction. Instead we investigate a collision operator that is defined by an underlying molecular dynamic simulation.

III Lattice Gas with Molecular Dynamics collision operator

In principle most systems of interest for a lattice gas (LG) simulation could be simulated using a Molecular Dynamics (MD) approach as well. MD is a standard tool which follows classical particle trajectories for particles interacting with a pair-potential by numerically integrating Newton’s equation of motion.

To construct a lattice gas method from a molecular dynamics simulation we overlay a lattice onto our MD simulation. The number of particles in each reciprocal lattice cell around the lattice position $x$ then corresponds to the lattice gas density $\rho(x)$, as shown in Fig. 1. If we then choose a time-step $\Delta t$ we can observe where particles ending up in cell $x$ came from. The number of particles moving from cell $x-v_{i}$ to cell $x$ then corresponds to the lattice gas occupation number $n_{i}(x,t)$. What is important to note is that this resulting lattice gas model is fundamentally correct in the sense that will obey the continuity and Navier-Stokes equations simply because the molecular dynamics simulation does so.

The initial condition for a set of $N$ particles in a finite container with periodic boundary conditions is given by their initial positions $x_{i}(0)$ and velocities $v_{i}(0)$. These particles then interact through an interaction pair-potential that we take here to only depend on the distance between the two particles: $V_{ij}=V(|x_{i}(t)-x_{j}(t)|)$. The MD simulation then provides (too good accuracy) the trajectories $x_{i}(t)$ which solve Newton’s second law:

[TABLE]

We now superimpose a lattice onto the computational domain of the MD simulation. For simplicity we can imagine a cubic lattice of lattice spacing $\Delta x$, although any lattice will do here. Let us define a function that determines if a particle resides in a specific cell of the reciprocal lattice associated with a lattice point $x$:

[TABLE]

Next we pick a time step $\Delta t$. We can now determine the lattice displacement each particle originally residing in a lattice point $x$ experiences. The set of all such displacements makes up the minimal set of lattice velocities $v_{i}$ for our lattice gas method, and the number of particles associated with this displacement makes up the lattice gas densities $n_{i}(x,t)$. We define

[TABLE]

This definition ensures that the particle numbers $n_{i}(x,t)$ will undergo a streaming step given by Eq. (6). For any given MD simulation we then know all $n_{i}(x,t)$. From Eq. (7) we see that the collision operator is then given by

[TABLE]

This fully defines the MDLG algorithm, a lattice gas with a collision operator that is defined through an underlying MD simulation. In some sense this is an ideal lattice gas model that can handle even the most complex situations, i.e. anything that can be addressed by MD, correctly. The key question is whether this collsion operator can be reduced to some stochastic collision operator that only depends on the local $n_{i}(x,t)$. Clearly this will only be the case for very simple systems since the MDLB collision operator contains information about temporal and spatial correllations of the underlying MD algorithm and can, in principle, deal with many complex phenomena like liquid-gas-solid coexistence, large varieties of transport parameters, including phenomena at high Knudsen, high Mach, and/or high Reynolds numbers, which we don’t expect to be accessible to a simple lattice gas algorithm of Eq. (7) with a local collision operator. Such extensions will be subject of future research, but are outside the scope of the current paper.

The local number of particles in lattice cell $x$ at time $t$ is give by

[TABLE]

This is consistent with the lattice gas definition of the local density because

[TABLE]

The last equality follows because

[TABLE]

i.e. every particle will be found somewhere on the lattice. Note that we have not yet restricted the velocity set. We will use as many velocities as needed. Mass conservation of Eq. (4) is clearly fulfilled since

[TABLE]

The definition of momentum in the lattice gas sense is typically defined as:

[TABLE]

However, relating this to the underlying momentum of the MD simulation is not exact, as can be seen in the example of a single MD particle moving with a velocity less that is not a lattice velocity shown in Fig. 2. The correspondence could be made exact if we were to introduced an average over all possible placements of the lattice. Such an average would make no difference to the global equilibrium distribution, which is the main focus of the remaining paper. We therefore avoid this additional complication for the current paper.

Similarly, momentum conservation of Eq. (5) is only exact if we introduce an average over lattice placements:

[TABLE]

Of course this does not mean that there is a problem with momentum conservation. Instead the problem arises due to the definition of momentum through measured mass transfer between sites for a fixed lattice.

Despite the apparent lack of momentum conservation the MDLG collision rules are still correct, even without the averaging, since the underlying MD simulation respects momentum conservation. As such the apparent violation of momentum conservation of the MDLG model is benign. We reserve a closer examination of this averaged lattice gas implementation for a followup paper.

The key question is then whether the collision operator (Eq. (12)) can take the form of Eq. (1), i.e. a stochastic collision operator that only depends on the current local occupation numbers $n_{j}(x,t)$. Since the is a whole ensemble of MD simulations that is consistent with a set of $n_{j}(x,t)$, and these different MD simulations will lead to different collision terms, it is clear that there can be no exact mapping. However, it is reasonable to hope that we will be able to construct a stochastic lattice gas collision operator that is statistically equivalent to the collision operators for the ensemble of corresponding molecular dynamics simulations. Establishing this is not a trivial task, and we will focus on the easier problem of showing that these collision operators are consistent with the equilibrium behavior of the lattice gas. In the next section we will present the lattice Boltzmann method which conceptually represents the ensemble average of a lattice gas method. Given the complexity of the task we focus in this paper on examining for which, if any, discretizations the MDLG and the standard lattice Boltzmann method give an equivalent equilibrium behavior.

IV MDLG for an two-dimensional Lennard-Jones gas

As a test case we use for our underlying MD simulations particles interacting with the standard Lennard-Jones interaction potential, which is given by

[TABLE]

The interaction strength is controlled by $\epsilon$ and the spatial scaling by $\sigma$. If $m$ is the mass of a particle we can construct the time-scale

[TABLE]

We performed a molecular dynamics simulation using LAMMPS for $100\,000$ particles in a two-dimensional box with the length of $1000\sigma$, corresponding to a nominal volume fraction of 0.0785 if we approximate the particles as circles with diameter $\sigma$. The system was initialized with a homogeneous distribution of particles with a kinetic energy corresponding to a temperature of 50 in the LJ units defined above. This corresponds to a gas at a high temperature, dense enough so that there are a significant number of collisions. The temperature is well above the critical temperature for a liquid-gas coexistence of $T_{c}=1.3120(7)$Potoff and Panagiotopoulos (1998). We ran an simulation of an equilibrium system with a time-step of $0.0001\tau$ for $3\,000\,000$ time-steps, i.e. up to $\tau=300$. Early time data of $1\,000\,000$ timesteps was discarded, to ensure that we were only probing the equilibrium dynamics.

We then analyse the resulting MD trajectories to obtain the resulting averaged MDLG occupation numbers $n_{i}(x,t)$ from Eq. (11). We should note here that results for different mean velocities $U$ can be obtained simply be using a lattice displaced by $-U\;\Delta t$ for the MDLG analysis. It is therefore not necessary to re-run the MD simulations to examine different mean velocities.

The first information to be gleaned from this is the resulting velocity set for the $v_{i}$. For small times $\Delta t$, only the nearest neighbors have non-negligible contributions, but as $\Delta t$ is increased more densities get populated. We identify the velocities $v_{i}$ using the numbering scheme shown in Fig. 3. So for $\Delta t\rightarrow 0$ only $v_{0}-v_{8}$ will have contributions. These velocities form a complete shell around the central point $v_{0}$. Most standard lattice Boltzmann methods work hard to make due with this minimal velocity set. This comes at some cost, the most important one is that only one temperature, $\theta=1/3$ in lattice unites, is allowable in lattice units to recover the correct viscous stress tensor (see Eq. (57)). For larger $\Delta t$ particle will travel further and a larger velocity set is required. The average occupation numbers are given by the global equilibrium distribution. The next subsection discusses how these equilibrium distribution can be obtained analytically.

IV.1 Global Equilibrium Distribution

For a system in thermal equilibrium, sufficient averaging will give an equilibrium distribution

[TABLE]

We can numerically approximate this equilibrium density by averaging the values of the $n_{i}$ from Eq. (11) over the whole lattice and for the duration of the simulation. For a given MD simulation the results will depend both on the lattice spacing $\Delta x$ and on the time step $\Delta t$. As mentioned above the MD simulations considered in this paper deal with fairly hot gases, that should be reasibably well approximated by an ideal gas.

To theoretically calculate this expectation value we assume that the particles are uniformly distributed, so only the displacement during the time-interval $\Delta t$ will enter the averaging:

[TABLE]

The key for our averaging will then be the probability of finding such a displacement $P(\delta x)$, and this allows us to write the average as

[TABLE]

where $W(x)$ is the d-dimensional wedge function defined as

[TABLE]

where $\Theta$ is the Heaviside function and $\alpha$ denotes cartesian coordinate index.

For very short times $\Delta t$, i.e. times shorter than the mean free time between two collisions, particles simply undergo ballistic motion. The velocity distribution of the particles is given by the Maxwell-Boltzmann distribution

[TABLE]

With this, and neglecting any collisions between the particles, we get for the mean squared displacement in one dimension

[TABLE]

The probability distribution for the displacement is then given by

[TABLE]

For times much longer than the mean free time particles undergo multiple collisions and instead of following a ballistic motion they will diffuse. If we call the self-diffusion constant $D$ we when have

[TABLE]

This implies that the probability of the discplacement is given by

[TABLE]

Now since both limiting displacements are given by Gaussian distributions is it reasonable to expect that the intermediate probabilities are also well approximated by a Gaussian and, if we know the mean squared displacement in one dimension $\langle(\delta x_{\alpha})^{2}\rangle$ (and assume isotropy), we get for the probability

[TABLE]

In all of these cases the probability distribution factorizes

[TABLE]

where

[TABLE]

and we can write Eq. (30) as a product of Gaussian integrals:

[TABLE]

The solution is given by

[TABLE]

where

[TABLE]

where

[TABLE]

This is a lattice equilibrium distribution function derived from first principles. At first glance it looks different than other lattice equilibrium distributions, and we will examine its relation to know equilibrium distribution functions below.

First we need to fully define the lattice equilibrium distribution. To do so we need to obtain the mean square displacement $\langle(\delta x_{\alpha})^{2}\rangle$. In general the mean square displacement can be measured in our MD simulations, but this would require us to consider a whole function of $\Delta t$ as in input parameter. For our simple ideal gas system we can obtain a simpler dependence on a single parameter by expressing it in terms of the velocity corellation function as

[TABLE]

For gases this velocity correlation function is typically well approximated by a simple exponential decay. There is also a long range $1/t$ contribution to the velocity correlation function for our two dimensional system, but for the times $\Delta t$ that are of interest here this divergent contribution does not yet contribute noticeably.

[TABLE]

For our system we compare this prediction of an exponentially decaying velocity correlation function to the measured corellation function in Fig. 4(a). We see that for early times we see good agreement with this prediction for $\tau=0.5(0)$. We also see the long time-tale typcical for a two dimensional system, which we ignore here. This is justified below Uhlenbeck and Ornstein (1930); Tsang and Tang (1977); Green (1952, 1954); Kubo (1957); Weitz et al. (1989).

Then the mean squared displacement can be predicted according to Eq. (47) as

[TABLE]

We show that this prediction recovers the measured mean squared displacement well in Fig. 4(b). Deviations resulting from the long time tails of the velocity corellation function only show up for later times and larger displacements considered in this paper, which justifies our ignoring these long time tails here.

This fully completes the definition of the MDLG equilibrium function in the case of gases. To verify our results we compare a numerically measured equilibrium distribution with the theoretically predicted one for different discretizations. The results are shown in Fig. 5. The agreement between our theoretical results and the experimental ones is excellent.

In the next section we introduce the lattice Boltzmann method and then examine the relation of this MDLG equilibrium function with existing lattice equilibrium distribution functions derived for lattice Boltzmann methods.

V Lattice Boltzmann

Lattice Boltzmann methods were derived as ensemble averages of lattice Boltzmann methods. The variables in a lattice Boltzmann method are distribution functions

[TABLE]

where the $\langle\cdots\rangle_{neq}$ represents a non-equilibrium ensemble average over microscopic lattice gas states. Taking the same ensemble average, the evolution equation for these lattice Boltzmann densities derives from the underlying lattice gas evolution Eq. (7)

[TABLE]

where the collision operator $\Omega_{i}=\langle\Xi_{i}\rangle$ is a deterministic function of all the densities at lattice point $x$. We will investigate later if this collision operator can be, at least approximately and for some suitable discretization, be cast in the standard BGK form typically employed for lattice Boltzmann simulations. This question is a crucial first step if one wants to relate lattice Boltzmann to an explicit molecular system which will be represented by our MDLB algorithm.

This standard LB collision operator is a first order BGK approximation and can be written as

[TABLE]

where $f_{j}^{0}$ is a local equilibrium distribution which depends only on the locally conserved quantities

[TABLE]

although other collision operators are also being used Ansumali and Karlin (2002); Ansumali et al. (2003); Ansumali and Karlin (2005); Geier et al. (2015) and it is a longer term goal of the MDLG method to help identify which of these collision operators are most realistic.

To ensure that the lattice Boltzmann equation reproduces the continuity and Navier-Stokes equations in the hydrodynamic limit it is necessary that the equilibrium distribution matches the first four (apart sometimes from a $u^{3}$ term) velocity moments of the Maxwell-Boltzmann distribution:

[TABLE]

where $Q_{\alpha\beta\gamma}$ should be zero. For velocity sets including only one shell we have $v_{i\alpha}\in\{-1,0,1\}$. For these velocity sets these moments overconstrain the equilibrium distributions. In particular we have $v_{i\alpha}^{3}=v_{i\alpha}$, which couples the first and the third moment. This is a key source of Galilean invariance violations in lattice BoltzmannWagner and Li (2006). These moments can only be reconciled for the special choice of

[TABLE]

and the third order tensor $Q_{\alpha\beta\gamma}=\rho u_{\alpha}u_{\beta}u_{\gamma}\ll 1$ which is assumed to be small because $u<0.1$ in typical situations. The equilibrium distribution is typically given in terms of an expansion in terms of the local velocity $u$ up to second order as

[TABLE]

where the weights $w_{i}$ depend on the velocity set and summation over repeated Greek indices is implied. In this article we focus on the question whether this form of an equilibrium distribution is compatible with a concrete MDLB implementation.

This collision operator together with the local equilibrium distribution implies mass and momentum conservation

[TABLE]

which is consistent with the typical conditions for lattice gases of Eqs. (4) and (5).

In the following we will examine the MDLG method for the example of a hot, dilute gas. For this lattice gas we examine the resulting distribution functions and see under which circumstances this lattice gas can reproduce (to some approximation) the lattice Boltzmann method equilibrium distribution Eq. (60).

VI Relation of MDLG equilibrium functions to lattice Boltzmann equilibria

We are now in a position to predict for which set of parameters $\Delta x,\Delta t$, if any, we can recover the traditional form of the lattice Boltzmann equilibrium from our MDLB algorithm. Most lattice Boltzmann methods use a limited velocity set that corresponds to a single shell in Fig. 3. For our two dimensions this corresponds to 9 velocities. The corresponding equilibrium distribution is typically given as a second order polynomial in the velocities, as we have presented earlier in Eq. (60). For the two dimensional D2Q9 lattice Boltzmann method we consider here the weights $w_{i}$ in Eq. (60) are given by

[TABLE]

where the velocity indices correspond to the numbering of Fig. 3.

In Fig. 5 we show the lattice Boltzmann weights $w_{i}$ as horizontal lines. To match the MDLG and LB equilibria we require

[TABLE]

This is a over-determined system of equations.

Fortuitously the solutions for the three distinct $w_{i}$ for a D2Q9 lattice Boltzmann give the same value for $a^{2}\approx 1/6$ to very good approximation. This is shown as the green vertical line in Fig. 5.

This suggests that matching lattice Boltzmann and MD simulations would likely benefit from using the conditions where the $f_{i}^{eq}$ match up, and the methodology explained above would give guidance on the appropriate time step $\Delta t$ for a given spatial discretization $\Delta x$. Given a $\Delta x$ can numerically solve Eq. (66) for $\Delta t$ for a system with zero mean momentum. This is shown in Fig. 6. We find that there is close agreement between the solutions for different velocities $v_{i}$. Corresponding to the transition from ballistic to diffusive regime around $\Delta t=\tau$ that we saw in Fig. 4(b) we also see a transition here from a $\Delta t\propto t$ regime for $\Delta t\ll 1$ to a $\Delta t\propto t^{2}$ regime for $\Delta t\gg 1$. We expect this relation that gives $\Delta t$ in terms of $\Delta x$ to be valuable when one tries to generate a coarse-graining transition between an MD and LG region in a multi-scale numerical method.

So far we have only matched the equilibrium distribution at zero velocity. The theory contains the mean velocity as $u$ in Eq. (43). For the measurements we could set up different simulations for mean velocities using the algorithm described above in Eq. (11). It is more practical, however, to move the Grid instead or equivalently use Galilean transformed particle positions of $\hat{x}_{i}(t)=x_{i}(t)+ut$ instead. Using this approach to find the equilibrium distribution for different mean velocities $u$ we show our comparison between the measured discrete equilibria (Eq. (11)), their theoretical prediction (Eq. (42)), and the lattice Boltzmann equilibrium distribution (Eq. (60)) in Fig. 7. For small velocities $|u|<0.1$ we find good agreement between all three densities. This is the relevant range, as lattice Boltzmann is only considered reliable for small enough velocities. The agreement between the measured and predicted MDLG equilibrium distributions continues to be excellent throughout the whole regime. Note, that the agreement would continue to be excellent for larger velocities. We would only need to adapt the velocity set we consider as velocities with magnitude larger than 0.5 are the same as velocities with magnitude larger than 0.5 plus an additional integer lattice displacement.

VI.1 Moments of the equilibrium distribution

The key property of the equilibrium distribution in kinetic theory are the velocity moments. For the derivation of the Navier-Stokes equation moments up to third order are required. It is therefore helpful to examine the moments of the discrete MDLG equilibrium distribution, compare them to the expected moment for a lattice Boltzmann equilibrium distribution and examine how these moments relate to the continuous velocity distribution function for the MD simulation.

Let us spend a moment considering these different concepts, since they are not usually clearly separated in a LB derivation. Firstly we have the velocity distribution of the MD simulation, given by Eq. (33). The moments of this velocity distribution are usually used as a rational for constructing LB equilibrium distributions such that the relevant moments of the discrete LB equilibrium distribution function match those of the continuous Maxwell Boltzmann distribution. This is the rational behind the moment Eqs. (55–58) for the lattice Boltzmann equilibrium function.

The moments of the discrete MDLG method are not a priori constrained to obey such a constraint, and indeed we don’t expect such an agreement for two reasons. First the underlying displacement probabilities (Eq. (38)) are in general different from the displacement probabilities (Eq. (35)) directly related to the Maxwell-Boltzmann distribution (Eq. (33)). Second there is no reason to believe that the averaging procedure of Eq. (41) will preserve the moments in general.

Now let us consider the first three velocity moments specifically. The zeroth moment relates to the total mass. Since the algorithm conserves mass exactly we expect all three approaches to agree on this moment. Indeed, as we saw in Eq. (22) mass is clearly conserved, and consequently this moment will agree for all of the above approach.

The first moment relates to the local momentum. Even if the MDLG approach does not locally conserve the momentum, the averaged momentum of the equilibrium distribution remains exact. This simply follows from the fact that this discrete moment corresponds to the net mass flow through the lattice. Even though this flow can be inexact at any instance in time (since particles may not cross a boundary despite the fact that they are moving) on average the count of particles crossing boundaries has to give the exact mass current.

Let us next consider the second moment. This second moment in the lattice Boltzmann approach (Eq. (57)) is related ideal gas equation of state $p=\rho\theta$. We calculate the second moment of our discrete MDLG equilibrium distribution

[TABLE]

For an equilibrium distribution that obeys the lattice Boltzmann moment (Eq. (57)) with a temperature $\theta=a^{2}$ this expression would give exactly one. For the MDLG equilibrium distribution this second moment is shown in Fig. 8. This shows that the discrete second moment does only agree with the MD temperature for large $a^{2}\gg 1$. As we see in Fig. 5 this corresponds to a situation where the populated set of velocities encompasses several shells in Fig. 3. For lower values of $a^{2}$ we find that the second moment diverges. The reason lies in the way we define the discrete equilibrium distribution. Even for very small $\langle(\delta x)^{2}\rangle$, corresponding to very small $\Delta t$, we will identify a fraction of particles that happen to cross from one lattice point to the next and are therefore assigned a lattice velocity $v_{i}$ of order one. This appearance of apparent large displacements causes the divergence of the discrete second moment. This effect is significantly enhanced by imposed velocites $u$.

For $a^{2}<0.1$ we see that we get significantly diverging values of $\Psi$ for different velocities $u$. This implies that this discrete second moment is not Galilean invariant. It is important to note that despite such violations of Galilean invariance of the discrete moments the full MDLG algorithm does not suffer from a Galilean invariance problem. Instead this is an indication that the collision operator Eq. (12) must exactly compensate this apparent Galilean invariance violation. For lattice Boltzmann methods one typically tries to avoid Galilean invariance violations by ensuring that both the local equilibrium distributions and the collision operator independently obey Galilean invariance.

For the lattice Boltzmann method we expect this second moment to be $\Psi_{\alpha\beta}=\theta/a^{2}\delta_{\alpha\beta}$. As we saw above [see Eq. (59)]lattice Boltzmann methods which have a velocity set consisting of a single shell in velocity space require $\theta=1/3$. This value of $\theta$ is consistent with the moment of the MDLG equilibrium for a value of $a^{2}\approx 1/6$. We find $\theta\approx 2a^{2}=\langle(\delta x)^{2}\rangle/(\Delta x)^{2}$. This corresponds to about the lowest value for $a^{2}$ where $\Psi$ does not strongly depend on $u$ and would therefore violate Galilean invariance.

The condition that $\theta=1/3$ came out of a consideration for the third moment for a minimal velocity set with $v_{ix}=v_{ix}^{3}$. We can define a third moment as

[TABLE]

We examine the behavior of $\Phi_{xxx}$ in Fig. 9. In fact this third moment should be zero, and for sufficiently large $a$ it converges to zero exponentially. However, if we artificially restrict our velocity set to a single shell, neglecting the small densities for discrete velocities outside the first shell, we find that there is a collapse of $\Phi_{xxx}$ to zero for the same $a^{2}\approx 1/6$ that we found in Fig. 8. For the full velocity set, however, nothing special occurs at $a^{2}=1/6$. This initially surprised us because we see in Fig. 5 that the densities associated with the second shell are still a factor of 30 smaller than the densities of the first shell and might therefore be taken to be negligible. The second shell velocities that are about a factor 2 larger than the first shell velocities, so these densities get multiplied by a factor of $2^{3}$ which allows these densities to contribute enough that we now have

[TABLE]

and they are different enough that the cancellation of terms that is supposed to lead to for $\Phi_{xxx}=0$ at $a^{2}=1/6$ no longer exists. Instead $\Phi_{xxx}$ monotonously approaches zero.

The second and third moments of the distribution functions both relate to the lattice Boltzmann temperature of Eq. (57). We can derive two expressions for the temperature:

[TABLE]

The dependence of these two quantities on the mean-square displacement is shown in Fig. 10. We see that these two definitions only agree with each other for large $a^{2}$, where agreement is facilitated by utilizing a large set velocity for the discretization of $f_{i}^{eq}$. For small $a^{2}$, and therefore a minimal velocity set, we get $\theta_{3}=1/3$, for the same reason that was discussed before, i.e. $v_{ix}=v_{ix}^{3}$.

If we artificially restrict the velocity set to a single shell (e.g. 9 velocities in our two-dimensional example) we affect both definitions of the LB temperature. In this case the agreement for large $a^{2}$ disappears, and we have exactly one point for $\theta=1/3$ corresponding to $a^{2}=1/6$ for which both expressions for the temperature agree. This corresponds to the special point in Fig. 9 where the results for the one-shell velocity set become independent of $u$ (for small enough $u$).

This suggests that there is a serendipitous agreement between the MDLG equilibrium distribution and the standard LB equilibrium distribution for one-shell velocity sets. It allows us to recover velocity moments up to order 3, which is exactly the order required by kinetic theory to recover the continuity and Navier-Stokes equations. As seen in Fig. 8 and 9, this value is just large enough to avoid apparent Galilean-invariance violations in the moments of the MDLG equilibrium distribution, and just small enough so that the values of the $f_{i}$ for velocities with $v_{i\alpha}>1$ are small enough not to contribute considerably to the moments of the equilibrium distributions. These effects can only be noticed for the third moment, where they conspire to induce agreement between the measures of temperature implied by the second and third moments, as seen in Fig. 10.

VII Consistent discretizations

Up to now we have seen how the moments of the equilibrium distribution depend on $a^{2}$ which is a measure of the mean squared displacement. We saw that the discrete moments of $f_{i}^{eq}$ differ from the continuous moments of the Maxwell-Boltzmann distribution even in the ballistic regime. The underlying reason for this disagreements result from two conspiring effects. Firstly we only know the position of our particles to lie somewhere within their assigned lattice cell. This uncertainty enters the definition of the $f_{i}^{eq}$ in Eq. (30). Secondly we use a discrete second moment. Here we show how both of these give an offset of $1/12$ in $a^{2}$ giving rise to a total shift of $1/6$ observed in the previous numerical results.

For simplicity let us consider large times away from the ballistic regime. This occurs without loss of generality if our assumption of a Gaussian distribution in Eq. (38) remains correct. We can then assume the motion of the particles is entirely diffusive with some diffusion constant $D$:

[TABLE]

If we had known the position of the particle initially, i.e. $\rho_{c}(x,0)=\delta(x)$, then the particle probability density would evolve as

[TABLE]

which has a second moment of

[TABLE]

If we only know that at time $t=0$ the particle is inside a lattice cell centered around the origin with width $\Delta x$, then the density $\rho$ as a function of time is given by

[TABLE]

which for $t\rightarrow 0$ gives a density that is $1/\Delta x$ inside the interval $[-\Delta x,\Delta x]$ and zero outside it. This probability distribution then spreads out and approaches a Gaussian at late times. The second moment of the position is then given by

[TABLE]

i.e. a simple offsett of $(\Delta x)^{2}/12$ that does not depend on time. We can identify

[TABLE]

from the definition of Eq. (44). We then we see that the results are expected to be shifted by $1/12$. However, in Fig. 10 we clearly see that this only represents half of the observed shift $1/6$.

The second effect relates to discretizing the position into displacement bins. We now calculate the discrete second moment (normalized by $(\Delta x)^{2}$)as

[TABLE]

In Fig. 11 we show that this discrete moment quickly converges to $a^{2}+1/6$. We clearly see that the missing additional offset of $1/12$ that we observed in Fig. 10 is the result of taking the discrete moment.

This shows that there are two effects of discretization. Firstly the initially broader distribution of particles confined to a lattice site rather than a point shifts the second moment by exactly 1/12. The second effect of discretization is more complicated, particularly for small $a^{2}$. But for $a^{2}>0.2$ this effect quickly converges to another offset in $a^{2}$ of $1/12$. Together they make up the offset of $1/6$, seen repeatedly in our numerical results of Figs. 5, 8, 9, and 10.

VIII Outlook

In this paper we have introduced a new tool for comparing the results of Molecular Dynamics simulations with those of coarse-grained lattice gas or lattice Boltzmann methods. It consists of re-interpreting the MD results as a lattice gas which we call the MDLG. The dynamics of this special kind of lattice gas is entirely given by the MD simulation, and therefore will be able to give a coarse-grained picture for any results that are obtainable with MD simulations.

Here we focus entirely on the averaged equilibrium behavior and show that there exists a close connection between the equilibria of lattice Boltzmann methods and the equilibrium for the MDLG method, when applied to a hot dilute gas. We were able to determine this equilibrium distribution analytically and were able to verify this analytical solution with the results of the MDLG method. Importantly there is a surprisingly good agreement between our equilibrium distribution and the standard lattice Boltzmann result for carefully chosen (and analytically known) pairs of time and space discretizations $\Delta t$ and $\Delta x$. We were able to understand the observed offset of $1/6$ in the dimensionless measure of the mean squared displacement $a^{2}$ in terms of our discretization procedure.

This opens the way for a more careful analysis of the fundamental underpinnings of lattice gas and lattice Boltzmann methods. We intend to utilize our MDLG method to investigate the fundamental properties of the collision operator, including its fluctuating properties. Further down the line we hope to investigate how the behavior of liquids alters the behavior of the MDLG method and examine if MDLG can also be matched with lattice Boltzmann methods. We anticipate that this method will also be instrumental in putting lattice Boltzmann methods for non-ideal and multi-component systems on a firmer footing.