Integer Lattice Gas with Monte Carlo collision operator recovers the entropic lattice Boltzmann method with Poisson distributed fluctuations

TL;DR

This paper introduces a Monte Carlo lattice gas model that mimics entropic lattice Boltzmann methods, providing insights into collision operators and fluctuations with Poisson statistics, advancing understanding of fluctuating lattice Boltzmann techniques.

Contribution

It presents a novel Monte Carlo lattice gas model that recovers the entropic lattice Boltzmann method and elucidates the origin of multi-relaxation time collision operators and Poisson fluctuations.

Findings

Derives equilibrium distribution for entropic lattice Boltzmann.

Shows Galilean invariant fluctuations follow Poisson statistics.

Provides insight into properties of fluctuating lattice Boltzmann methods.

Abstract

We are examining a new kind of lattice gas that closely resembles modern lattice Boltzmann methods. This new kind of lattice gas, that we call a Monte Carlo Lattice Gas, has interesting properties that shed light on the origin of the multi-relaxation time collision operator and it derives the equilibrium distribution for entropic lattice Boltzmann. Furthermore these lattice gas methods have Galilean invariant fluctuations given by a Poisson statistics, giving further insight into the properties that we should expect for fluctuating lattice Boltzmann methods.

Tables

Table 1. Table 1: Representation of the correlator of Eq. ( 117 ) for pairs of ⟨ n i ​ n j ⟩ delimited-⟨⟩ subscript 𝑛 𝑖 subscript 𝑛 𝑗 \langle n_{i}n_{j}\rangle for the worst statistical outlier for u x = − 0.68 subscript 𝑢 𝑥 0.68 u_{x}=-0.68 of Figure 10 .

	$n_0$	$n_1$	$n_2$	$n_3$	$n_4$	$n_5$	$n_6$	$n_7$	$n_8$
$n_0$	0.998863	-0.000028	-0.000482	-0.000170	-0.000343	0.000001	-0.000114	-0.000082	0.000010
$n_1$	-0.000028	0.999171	0.000001	0.000143	-0.000072	-0.000127	-0.000047	0.000013	0.000057
$n_2$	-0.000482	0.000001	0.999809	-0.000112	0.000177	-0.000515	-0.000110	-0.000022	0.000482
$n_3$	-0.000170	0.000143	-0.000112	1.000230	-0.000079	0.000014	0.000008	-0.000048	-0.000067
$n_4$	-0.000343	-0.000072	0.000177	-0.000079	0.999648	0.000492	-0.000023	-0.000103	-0.000480
$n_5$	0.000001	-0.000127	-0.000515	0.000014	0.000492	-0.999567	-0.000065	-0.000046	0.000593
$n_6$	-0.000114	-0.000047	-0.000110	0.000008	-0.000023	-0.000065	0.999801	-0.000085	0.000086
$n_7$	-0.000082	0.000013	-0.000022	-0.000048	-0.000103	-0.000046	-0.000022	0.999687	-0.000117
$n_8$	0.000010	0.000057	0.000482	-0.000067	-0.000480	0.000593	0.000086	-0.000117	0.999215

Equations

$$f_i(x+v_i,t+1)=f_i(x,t)+{\Omega }_i$$

$${\Omega }_i=\sum _j{\Lambda }_{ij}[f_i^0-f_i(x,t)]+{\xi }_i$$

$$\rho$$

$$\rho u_{\alpha }$$

$$\rho$$

$$\rho u_{\alpha }$$

$$\sum _i{f}_i^0(v_{i\alpha }-u_{\alpha })(v_{i\beta }-u_{\beta })$$

$$\sum _i{f}_i^0(v_{i\alpha }-u_{\alpha })(v_{i\beta }-u_{\beta })(v_{i\gamma }-u_{\gamma })$$

$$f_i^0(\rho ,u)=\rho w_i\left(1+\frac{v_{i\alpha }{u}_{\alpha }}{\theta }+\frac{v_{i\alpha }{u}_{\alpha }{v}_{i\beta }{u}_{\beta }}{2{\theta }^2}-\frac{u_{\alpha }{u}_{\alpha }}{2\theta }\right).$$

$$\sum _i{w}_i$$

$$\sum _i{w}_i{v}_{i\alpha }{v}_{i\beta }$$

$$\sum _i{w}_i{v}_{i\alpha }{v}_{i\beta }{v}_{i\gamma }{v}_{i\delta }$$

$$M^a=\sum _i{m}_i^a{f}_i$$

$$\sum _a{w}_i{m}_i^a{m}_j^a={\delta }_{ij}.$$

$$f_i=\sum _a{w}_i{m}_i^a{M}^a$$

$$\sum _i{w}_i{m}_i^a{m}_i^b={\delta }^{ab}.$$

$${\Lambda }^{ab}=\sum _{i,j}{m}_i^a{\Lambda }_{ij}{w}_j{m}_j^b=\frac{1}{{\tau }^a}{\delta }^{ab}.$$

$$v_{ix}\in \{-1,0,1\}.$$

$$w_0$$

$$w_{-1}=w_1$$

$$M=\left(\begin{array}[]{rrr}1&1&1\\ -\sqrt{3}&0&\sqrt{3}\\ \sqrt{2}&-\sqrt{\frac{1}{2}}&\sqrt{2}\end{array}\right)$$

$$v_i=\left[\begin{array}{c}{v}_x\\ v_y\end{array}\right]$$

$$w_{(v_x,v_y)}=w_{v_x}{w}_{v_y}$$

$$w_0$$

$$w_{1-4}$$

$$w_{5-8}$$

$$M=\left(\begin{array}[]{ccccccccc}1&1&1&1&1&1&1&1&1\\ 0&\sqrt{3}&0&-\sqrt{3}&0&\sqrt{3}&-\sqrt{3}&-\sqrt{3}&\sqrt{3}\\ 0&0&\sqrt{3}&0&-\sqrt{3}&\sqrt{3}&\sqrt{3}&-\sqrt{3}&-\sqrt{3}\\ 0&\frac{3}{2}&-\frac{3}{2}&\frac{3}{2}&-\frac{3}{2}&0&0&0&0\\ 0&0&0&0&0&3&-3&3&-3\\ -1&\frac{1}{2}&\frac{1}{2}&\frac{1}{2}&\frac{1}{2}&2&2&2&2\\ 0&-\sqrt{\frac{3}{2}}&0&\sqrt{\frac{3}{2}}&0&\sqrt{6}&-\sqrt{6}&-\sqrt{6}&\sqrt{6}\\ 0&0&-\sqrt{\frac{3}{2}}&0&\sqrt{\frac{3}{2}}&\sqrt{6}&\sqrt{6}&-\sqrt{6}&-\sqrt{6}\\ \frac{1}{2}&-1&-1&-1&-1&2&2&2&2\end{array}\right)$$

$$M^{a}=\sum_{i}m_{i}^{a}f_{i}=\left(\begin{array}[]{c}\rho\\ j_{x}\\ j_{y}\\ \sigma_{-}\\ \sigma_{xy}\\ \sigma_{+}\\ q_{x}\\ q_{y}\\ t\end{array}\right)$$

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

$$f_i^{eq}(\rho ,u\equiv 0)=\rho w_i.$$

$$w_i{w}_j{P}_{ij\to kl}=w_k{w}_l{P}_{kl\to ij}.$$

$$\frac{P_{ij\to kl}}{P_{kl\to ij}}=\frac{w_k{w}_l}{w_i{w}_j}.$$

$$P_{ij\to kl}\propto \min \left(1,\frac{w_k{w}_l}{w_i{w}_j}\right){\delta }_{(v_i+v_j),(v_k+v_l)}$$

$$\sum _{kl}{P}_{ij\to kl}=1$$

Full text

Lattice Gas, Lattice Boltzmann, Monte Carlo, Fluctuatons

I Introduction

The origin of fluctuations can be traced back to the discrete nature of matterEinstein (1906). When one examines a small enough systems, these fluctuations, originating from the stochastic nature of the dynamics of the discrete particles, manifest themselves. With few exceptions, e.g. the dynamics close to a critical pointKawasaki (1970), these fluctuations become irrelevant for larger systems, and the dynamics of this large-scale system becomes deterministic, allowing us to describe the system as a continuum. These continuum equations, like the continuity and Navier-Stokes equations for fluids, are enormously successful at describing macroscopic phenomena. So it is tempting to continue to utilize these continuous equations of motions down to scales where fluctuations become important. In these cases fluctuations, that were eliminated, have to be somehow re-introduced. One way to consistently introduce fluctuations is the Langevin approachLemons and Gythiel (1997), which consists of introducing a fluctuating term in the equations of motion and then adjusting the amplitude of the fluctuating term to give correct fluctuations in equilibrium. This is easily done since it is known that the states of the system obey a Boltzmann statistics.

This general narrative is the template for the development of fluctuating lattice Boltzmann methods. The original lattice gases introduced by Frisch, Hasslacher and Pomeau Frisch et al. (1987) consisted of a hexagonal lattice and particles moving along the links between nearest neighbor sites. The model was restricted to having at most one particle moving along each link. While it could be shown that this model had enough symmetry to recover the Navier-Stokes equations in the hydrodynamic limit, the model contained very large fluctuations, that needed to be averaged out to examine hydrodynamic phenomena. Part of the derivation of the Navier-Stokes equations required taking a formal ensemble average of the lattice gas, leading to a Boltzmann equation for the lattice gasFrisch et al. (1987). Mc Namara McNamara and Zanetti (1988) then realized that one could use this Boltzmann averaged lattice gas as a method in its own right. This method became known as the lattice Boltzmann method. The original lattice Boltzmann method was an exact ensemble average of the lattice gas, retaining the unconditional stability of the lattice gas as well as several of its flaws. At the time a key advantage of the new lattice Boltzmann method was the elimination of all noise.

Soon it was discovered that the lattice Boltzmann method could be simplified by using a Bhatnagar-Gross-Krook (BGK) approximation for the collision operator Higuera et al. (1989). The idea here is that the distribution function will approach the local equilibrium distribution, consistent with the locally conserved density and momentum. This significantly simplified the collision operation and at the same time removed some undesirable artifacts that had survived the transition from the lattice gas to the lattice Boltzmann methods, like the velocity dependence of the viscosityFrisch et al. (1987). While the improvements of these novel lattice Boltzmann methods made them extremely popular any application that needed fluctuations was stuck with lattice gas methods. Tony Ladd had utilized the fluctuations of lattice gases to simulate the Brownian motion of colloids Ladd et al. (1988). He then developed a fluctuating version of the lattice Boltzmann method by including a fluctuating stress tensor that would recover the fluctuating hydrodynamic equations Ladd (1993). While this recovered the correct hydrodynamic limit, the fluctuations could be seen to be only correct for wavelength of the size of the system. Fluctuations at smaller scales were suppressed. The reason for this became clear when Adhikari et al. noticed that all kinetic modes, not only those related to the hydrodynamic quantities, needed to have fluctuations addedAdhikari et al. (2005). The theoretical framework for this was work on generating a fluctuating version of the linearized Boltzmann equation by Bixon and ZwanzigBixon and Zwanzig (1969). Much work has gone into examining the correct form of these flucuations for lattice Boltzmann since thenDünweg et al. (2007); Kaehler and Wagner (2013a); Wagner and Strand (2016). The remaining difficulties with defining consistent fluctuations for lattice Boltzmann systems also hinder the progress for extended version of the fluctuating lattice Boltzmann method for non-ideal systems Gross et al. (2010); Ollila et al. (2011); Thampi et al. (2011); Gross and Varnik (2012); Belardinelli et al. (2015).

One obvious extension to lattice gases that maintains their discrete properties, is to extend the occupation numbers from booleans to integers. In 1988, just one year after the FHP paper Frisch et al. (1987), Kim Molvig et al. pioneered lattice gases that were able to remove lattice artifacts in the pressure and viscosity Monaco (1988), although the results are more sketched than presented in this article. This article was later credited by Chen et al. to have introduced an approach called “Digital Physics”. This approach was the idea behind the founding of the EXA company in 1991. There are some indications of further developments of this approach, notably the transition from boolean to integer occupation numbers, in subsequent papers Chen et al. (1997); Teixeira (1997); Anagnost et al. (1997); Halliday et al. (1999); Succi et al. (1999) but most of the details of the approach remain unpublished since they formed the basis for proprietory technology. The paper by Chen Chen et al. (1997) sketches an approach of extending boolean lattice gases to lattice gases with integer occupation numbers reminiscent of the one presented here, although it appears that it concerns an algorithm that appart from mass and momentum also conserves energy, rather than the approach presented here that, like most standard LB approaches, replaces conservation of energy with an isothermal approach. We will comment further on the similarities and differences later in the paper.

The integer lattice gas approach was independently re-introduced by Boghosian et al. in 1997 Boghosian et al. (1997). The focus of this paper is the thermodynamics behavior of integer lattice gases with the finite number of bits. While the theoretical part focuses on energy conserving systems, the explanation of the sampling method for the collision considers only mass and momentum conservation. The collision process then consists of a non-trivial sampling of points from a polyotope of allowed states, the details of which remain unpublished.

One year later Chopard and Masselot Chopard et al. (1998); Masselot and Chopard (1998) describe a different integer lattice gas approach: they impose a collision operator that closely mimicks a standard lattice Boltzmann BGK collision operator with an equilibrium distribution that is quadratic in the velocities. The result of this LB inspired collision process is then a continuous probability distribution that is sampled to obtain a new discrete distribution. The sampling process, however, will lead to a violation of momentum conservation, so an additional random walk process is required until a distribution is found that is consistent with the original momentum. This method shows a surprising density dependence of the viscosity Masselot and Chopard (1998) that the authors attribute to fluctuation effects.

Because of the rise of lattice Boltzmann, however, the development of lattice gases was then significantly curtailed. Even the EXA company switched to lattice Boltzmann implementations, and little has been done on the development of integer lattice gases since these early days. During his Diplomarbeit, Marin Geier worked on integer equivalents of the lattice Boltzmann approach, but he judged this early efforts as mostly unsuccessful. Some remnants of the integer representation can be seen in Geier et al. (2006), but here the appearance of integer rather than continuous distributions is incidental rather than central to the approach and was given up in further developments of the method in favor of continuous distribution functions.

Our renewed interest in integer lattice gas methods stems from our interest in fluctuations. We wanted to get away from the approach of introducing fluctuations in a continuous system through a Langevin approach, and back to a more natural approach of obtaining fluctuations as a direct result of the discrete nature of the fluid we are examining. The idea of using a lattice gas as an starting point to derive fluctuating lattice Boltzmann methods was pioneered by Duenweg Dünweg et al. (2007) who proposed a lattice gas with integer (rather than boolean) occupation numbers and Monte Carlo collision operator based on a zero-velocity equilibrium distribution as a starting point to derive fluctuating lattice Boltzmann methods. But this approach remained entirely theoretical without any attempt at an actual implementation, similar to the original introduction of lattice Boltzmann as a theoretical tool to derive the hydrodynamic limit of lattice gas methods.

The aim of this paper is to derive a discrete lattice gas version that corresponds to current state of the art lattice Boltzmann methods and which would have (at least in close approximation) these lattice Boltzmann methods as its ensemble average and its integer implementation shows Galilean invariant fluctuations of an ideal gas. This means the occupation numbers should be Poisson distributed in equilibrium. The results are very encouraging: we show below that the Monte Carlo Lattice Gas introduced in this paper recovers the equilibrium distribution of the entropic lattice Boltzmann method Ansumali et al. (2003); Chikatamarla et al. (2006) as its ensemble average, and recovers Galilean invariant fluctuations at even better accuracy than the approach by Kaehler et al. Kaehler and Wagner (2013a). The collision operator is to first order a multi-relaxation time BGK approach, but it has additional second order terms.

This paper is structured as follows: we first introduce the basic lattice Boltzmann method and then use these result to propose a lattice gas method with a Monte Carlo collision operator. We then show how this is practically implemented in one and two dimensions. The Boltzmann limit of this lattice gas implementation is derived and we show that it recovers the entropic lattice Boltzmann method (rather than a lattice Boltzmann method with a quadratic equilibrium that inspired it) in the limit of small deviations from local equilibrium. The last section shows that this implementation indeed recovers the independent Poisson statistics for the densities exactly; this is the statistics of an ideal gas that we previously tried to impose on lattice Boltzmann methods, but only recovered approximately Kaehler and Wagner (2013a).

II Basic lattice Boltzmann

The lattice Boltzmann algorithm consists of continuous densities $f_{i}$ that are associated with velocities $v_{i}$ which are lattice velocities. This means that if $x$ is a lattice position, $x+v_{i}$ will also be a lattice position. These densities are defined on each lattice point $x$ at discrete times $t$ which are taken to be integer values. They evolve in time according to the lattice Boltzmann equation

[TABLE]

where $\Omega_{i}$ is the collision operator. In its multi-relaxation-time (MRT) BGK form this collision operator can be written as

[TABLE]

where $f_{i}^{0}$ is the local equilibrium distribution and $\xi_{i}$ is a noise term. Most lattice Boltzmann methods are conserving mass and momentum, but instead of conserving energy they are coupled to a heat bath. We define the local density $\rho$ and momentum density $\rho u$ through the velocity moments of the $f_{i}$ as

[TABLE]

Mass and momentum conservation is ensured by the requirement that the local equilibrium distribution also obey

[TABLE]

For the recovery of the Navier-Stokes equations we also require higher order moments for the equilibrium distribution. From analogy to the velocity moments of the Maxwell-Boltzmann distribution we demand

[TABLE]

where $Q_{\alpha\beta\gamma}$ should be zero. Unfortunately for small velocity sets typically used which have $v_{i\alpha}^{3}=v_{i\alpha}$ because $v_{i\alpha}\in\{-1,0,1\}$ this is not possible Wagner and Li (2006). But for the special choice of $\theta=1/3$ it can be reduced to $\rho u_{\alpha}u_{\beta}u_{\gamma}$, which is assumed to be small for $u\ll 1$ found in lattice Boltzmann simulations. An expansion of the Maxwell Boltzmann distribution to second order in velocities then gives

[TABLE]

The $w_{i}$ are obtained by matching the discrete velocity moments to the continuous velocity moments of the equilibrium distribution. To recover the continuity and Navier-Stokes equations we require the first four moments. All odd moments are zero and the even moments obey

[TABLE]

The collision matrix $\Lambda_{ij}$ is constructed so that the stress moments are relaxed at a rate that determines the viscosity, and other moments can be relaxed at different rates d’Humieres (1994); Lallemand and Luo (2000), e.g. to optimize the stability of the method. In general we can write the moments as

[TABLE]

where the first moments will be related to the hydrodynamic moments. For completeness we require as many moments $M^{a}$ as we have densities $f_{i}$ so we can have a one-to-one correspondence between moment and velocity space. It is often useful to require that the square matrix generating the moments be orthogonal with respect to some measure. Particularly for fluctuating applications it is often Adhikari et al. (2005); Kaehler and Wagner (2013a, b) found to be advantageous to require

[TABLE]

This implies the backtransform

[TABLE]

as well as

[TABLE]

This implies that moments $m_{i}^{a}$ are constructed startig from the hydrodyanmic moments of interest (density, momentum, stress tensor), complemented by a set of ghost-modes, which are then orthonormalized with a Gram-Schmidt orthonomalization scheme using the scalar product implied by (14). In this representation the collision matrix is designed to be diagonal:

[TABLE]

Setting the $\tau^{a}$ then fully determines the algorithm. Next we will briefly discuss two common lattice Boltzmann velocity sets for one and two dimensions.

II.1 D1Q3

The minimal lattice lattice Boltzmann method in one dimension consists of only three velocities corresponding to particle moving to a lattice site to the right, the lattice site to the left or a particle remaining at its lattice site. This is referred to as the D1Q3 model. Our velocity set is given by

[TABLE]

The weights are given by

[TABLE]

The moment matrix is given by

[TABLE]

The first two rows correspond to the conserved mass and momentum moments respectively. So there is only one relaxation time for the moment related to the third row for the D1Q3 model.

II.2 D2Q9

In the D2Q9 implementation, our velocity set consists of 9 velocities given by the product set of the D1Q3 velocities of Eq. (18):

[TABLE]

where $v_{x},v_{y}\in\{-1,0,1\}$. The weights in the two-dimensional case are just products of the one-dimensional weights

[TABLE]

or, using the numbering of Fig. 1, the more standard

[TABLE]

The moment matrix is given by

[TABLE]

where the numbering of the rows and columns is given by the numbering of the velocities in Figure 1. We then get the moments

[TABLE]

This concludes our very brief recap of the lattice Boltzmann algorithm (in its most common form) and we will next present the Monte Carlo lattice Gas algorithm introduced in this paper.

III The Monte Carlo Lattice Gas algorithm

The main idea behind defining a lattice gas that behaves equivalently to a lattice Boltzmann method is that we will use a Monte Carlo collision operator that explicitly conserves mass and momentum and will recover the correct equilibrium distribution for $u=0$. At each lattice site $x$ the number of particles that streamed in from a lattice position $x-v_{i}$ is denoted as $n_{i}(x,t)$. This number is an integer, in contrast to the real number $f_{i}(x,t)$ in the lattice Boltzmann approach. These particles then are re-distributed due to collisions, and this re-distribution will be denoted as $\Xi_{i}$. We then have the evolution equation We write the evolution equation as

[TABLE]

which looks equivalent to the lattice Boltzmann equation (1), except that the occupation numbers $n_{i}$ are integers, and that the collision term is an inherently probabilistic term, not a deterministic one as in the lattice Boltzmann approach. We consider that this collision term is the accumulated effect of many two-particle collisions. This may not be the most numerically efficient algorithm, but it avoids complexities encountered by Boghosian Boghosian et al. (1997) and Duenweg Duenweg and Bellas (2016) because it is conceptually clean. It may be similar to the approach of Digital Physics, as hinted at in Chen et al. (1997), but we have been unable to uncover the details of this approach.

To define the collision term $\Xi_{i}$ we demand that it recover the zero-velocity equilibrium distribution of lattice Boltzmann. From Eq. (9) we see that it is given simply by

[TABLE]

Let us now consider a collision of two particles with velocities $v_{i}$ and $v_{j}$. We denote the probability of colliding these two particles and ending up with two particles with velocities $v_{k}$ and $v_{l}$ with $P_{ij\rightarrow kl}$. We assume detailed balance, which means that the forward and backward collisions times the probabilities of finding these pairs in equilibrium have to be equal. Let us a assume that the number of particles at a lattice site is $N$. At one lattice site the probability of picking a particle with velocity $v_{i}$ is then $w_{i}$.

We can then write the detailed-balance condition for the equal probability for the forward and backward collisions as

[TABLE]

This fixes the ratios of the forward and backward collisions to

[TABLE]

Assuming that we do not have any additional conserved quantities this condition will ensure that the system will approach the lattice Boltzmann equilibrium distribution if the system has a mean velocity of zero. For actual collisions, i.e. $ij\neq kl$, we pick the transition probability

[TABLE]

where $\delta$ is the Kronecker delta. Note that this ensures the probability ratio requirement of Eq. (31), which is always guaranteed for the case $ij=kl$. Next we have to ensure that the probabilities add up to one:

[TABLE]

and of course we also require $P_{ij\rightarrow kl}\geq 0$. Formally we can achieve this by introducing a proportionality factor $\lambda_{ij,kl}$ to get

[TABLE]

and requiring $\lambda_{ij,kl}=\lambda_{kl,ij}$ to ensure Eq. (31). The $\delta$-function ensures that only collisions that conserve momentum (mass conservation is trivially true) have a non-zero probability. Thus the collision will conserve mass and momentum. Note that the lower condition, together with the requirement of positive probabilities, puts an constraint for the largest achievable $\lambda_{ij,kl}$.

The lattice gas algorithm then consists of moving particles according to their associated velocities $v_{i}$ and then picking $C$ pairs of particles at random and colliding them according to the transition probabilities given above. In order to maintain conceptual clarity we use the naive, straight forward, non-optimized algorithm. As explained above, we denote the number of particles at position $x$ and time $t$ moving in direction $v_{i}$ with $n_{i}(x,t)$. The local density is then given by

[TABLE]

At each lattice site we then pick a pair of particles to collide by selecting two evenly distributed random numbers between 1 and $N$. Given these two random numbers $r_{1}$ and $r_{2}$ we identify the identity of these particles by

[TABLE]

This then implies two different random numbers $s_{1}$ and $s_{2}$, corresponding to the number of the two selected velocities. The two random numbers then imply picking two particles with velocities $v_{s_{1}}$ and $v_{s_{2}}$. For a collision of this pair of particles we then pick two outgoing velocities $v_{s_{3}}$ and $v_{s_{4}}$ with probability $n_{s_{1}}n_{s_{2}}P_{s_{1}s_{2}\rightarrow s_{3}s_{4}}/N^{2}$. The effect of the collision number $c$ is then given by a random variable $\vartheta_{i}(s_{1},s_{2},s_{3},s_{4})$ that takes on the value

[TABLE]

which corresponds to the change in the number of particles corresponding to velocity $v_{i}$ after the collision $(v_{s_{1}}v_{s_{2}}\rightarrow v_{s_{3}}v_{s_{4}})$ has been chosen. This process is repeated $C$ times, and we get the full collision operator as the sum of these random variables:

[TABLE]

where it is understood that each of these sub-collision operators $\vartheta_{i}(s_{1},s_{2},s_{3},s_{4})$ is the result of a collision for the state given by the previous sub-collisions, according to Eq. (36). Note that this collision operator explicitly conserves mass and momentum because of Eq. (32).

Next we will show two explicit implementations of this algorithm in one and two dimensions for standard velocity sets.

III.1 A D1Q3 implementation

A state is given by the numbers of particles $n_{i}$ moving in direction $v_{i}$. For this velocity set the only collision that conserves momentum and changes the state of the system consists of two particles moving in opposite directions that come to rest, or the inverse process of two rest particles that will move appart in opposite directions.

[TABLE]

For the D1Q3 lattice Boltzmann method the weights are given by Eq. (19). We get for the transition probabilities

[TABLE]

where we have set some irrelevant switching probabilities to zero since they do not change the state.

As a consequence of this small collision set, an accidentally conserved quantity is introduced. The number of rest particles at each lattice site will stay either even or odd throughout the entire simulation, as they can only be created or destroyed in pairs. This conserved quantity prevents us from exactly recovering the equilibrium fluctuations, particularly the Poisson distributions. This is shown in Fig. 2 for a worst-case scenario where all lattice sites have an even number of rest-particles. For this simulation, we used $L_{x}=100$, $C=10$, $\rho=12$, allowed the simulation to equilibrate for 20000 time steps, and averaged over 80000 time steps.

III.2 A D2Q9 implementation

The velocities of Eq. (21) allow for many more collisions, and these collisions can be grouped into equivalence classes, in which each collision within a class is simply a rotation of the other collisions. In order to ensure isotropy, we must ensure that the $\lambda_{ij,kl}$ for all collisions in an equivalence class are the same.

Some collisions in our D2Q9 model break the conservation of ‘even or odd-ness’ of $n_{0}$ at each lattice site. However, in the projection of our simulation along the x or y direction, we still have that each row and column is constant in it’s ’even or odd-ness’. This has a much smaller impact on the simulation as our lattice size grows. Let us assume the lattice has $L_{x}$ lattice sites in the x-direction and $L_{y}$ lattice sites in the y-direction. Then the number of conserved quantities grows only as $L_{x}+L_{y}$, but the number of lattice sites grows as $L_{x}L_{y}$, so the fraction of spuriously conserved degrees of freedom becomes small for large lattices.

The equivalence classes are shown schematically in Fig. 3. As a small side note, the the collision 9 was first proposed by Molvig Monaco (1988). This allows us to write down all possible binary collisions and their probabilities in principle, exactly in the same way as we did for the D1Q3 model. We make the simplification that the relaxation rates $\lambda_{ij,kl}$ are instead referred to as $\lambda_{i}$, where each index $i$ refers to the number of the equivalence class identified in Fig. 3. We now have many more collisions, which we will not explicitly list their probabilities here. It should also be noted that there are three and more particle collisions that could be included. For simplicity those are neglected here.

IV Boltzmann average of the lattice gas

To compare the lattice gas results to the lattice Boltzmann method we will now examine an non-equilibrium ensemble average of the lattice gas evolution equation. We define the particle probability densities as

[TABLE]

where $\langle\rangle$ implies a non-equilibrium average over an ensemble of microscopic realizations leading to the same macroscopic state. We define the density as

[TABLE]

The evolution equation for the average particle densities is then

[TABLE]

The collision operator is the averaged version of the lattice gas collision operator

[TABLE]

We can define intermediate distribution functions recursively as distribution functions affected by one single collision

[TABLE]

where it is understood that

[TABLE]

The averaged collision operator can then be written as

[TABLE]

From this collision operator we can obtain the equilibrium distributions for the $f_{i}$ by demanding

[TABLE]

Note that it is sufficient for this purpose to consider a single sub-collision of Eq. (55), rather than the full collision operator of Eq. (57). This gives us a quadratic matrix equation in the equilibrium density. One might expect not to obtain a unique solution. However, since we require that $f_{i}^{eq}>0$ we will find below that only one physical solution survives.

IV.1 D1Q3 results

In the one-dimensional case of D1Q3 we find that Eq. (57) gives

[TABLE]

The equilibrium distribution has to be found with respect to the corresponding conserved quantities of mass $N(x,t)$ and the momentum

[TABLE]

This corresponds to an ensemble averaged momentum of

[TABLE]

Furthermore we define a second moment

[TABLE]

We can use a matrix that transforms the $f_{i}$ vector onto the conserved moments and one free moment of Eq. (20).

We then have the moments

[TABLE]

We can then write the collision operator for a single collision of Eq. (55) in this moment space as

[TABLE]

In equilibrium the collision operator does not change the distributions. Reversely, we can find the equilibrium distribution by demanding that the collison operator be zero. Thus we obtain $\pi^{eq}$ in terms of the conserved quantities as

[TABLE]

The solution with the negative sign leads to negative densities, while the other solution stays positive. We obtain the D1Q3 equilibrium distribution

[TABLE]

which can be written as

[TABLE]

where we introduced a weight $w_{v_{ix}}$, anticipating that the weights for higher dimensions can be written as products of one-dimensional weights. For $u=0$ we recover $f_{i}^{eq}=\rho w_{i}$, which was our starting point in selecting the collision probabilities. Once the collision probabilities are fixed, they also imply equilibrium distributions for different conserved velocities. This equilibrium distribution, and its dependence on $u$, is therefore the logical consequence of selecting the collision probabilities corresponding to the imposed zero velocity equilibrium distribution function.

We show a comparison of a measured equilibrium distribution to this prediction in Fig. 4. For the measured equilibrium distribution we show results for $\langle N\rangle=36$, $\lambda=1$, a lattice size of $L_{x}=100$ with periodic boundary conditions. We used $C=10$ collisions per iteration step and the simulation was initially run for $15,000$ steps to equilibrate the simulation, and then another $15,000$ iterations were used for the measurements. We find excellent agreement between this prediction and our measurements of the equilibrium distribution. The agreement between our new equilibrium distribution and that of a standard LB equilibrium distribution function of Eq. (9) is also close for $|u|<0.4$, but diverges thereafter.

This difference implies that the second moment of the equilibrium distribution will not obey the lattice Boltzmann requirement of Eq. (7) for all $u$. This is not surprising since the MCLG equilibrium distribution, unlike it’s LB equivalent, is positive definite. This means that for a velocity of $u=\pm 1$ we require $f_{\pm 1}=\rho$, and therefore the second moment must be zero for this extreme value. However lattice Boltzmann simulations typically require $u\ll 1$, which is usually taken as $u<0.1$, so disagreements between these equilibrium distributions outside the range $|u|<0.1$ have little practical relevance.

We show the second moment of Eq. (7) for our MCLG equilibrium distribution in Fig. 5. As predicted, the effective temperature for the MCLG method goes to zero for large absolute velocities, as is unavoidable for a discrete method with a restricted velocity set that has positive definite occupation numbers. For the relevant region where $|u|<0.1$, however there is excellent agreement between the MCLG and LB results, as shown in the inset.

After identifying the global equilibrium distribution for the D1Q3 MCLG, we should now examine the collision operator more closely. Eq. (64) implies that the collision operator leaves the mass and momentum modes unchanged (as it has to, since those are conserved) and only alters the $\pi$ mode. Let us define its deviation from equilibrium as $\tilde{\pi}=\pi-\pi^{eq}$. If we denote $\tilde{\pi}^{(n)}$ as the value of $\tilde{\pi}$ after $n$ collisions we can write the effect of one collision as

[TABLE]

If we interpret this as a non-linear differential equation we can write the analytical solution after $C$ collisions as

[TABLE]

For the standard BGK approach in lattice Boltzman we would have expected a pure exponential decay, corresponding the the approximation of neglecting $\tilde{\pi}^{2}$ in Eq. (68). The connection between the LG collision operator and the LB BGK collision operator is interesting. The collision MCLG operator is quadratic (and would also involve higher powers if we allowed for multi-particle collisions), but the BGK operator is linear. This apparent mismatch is resolved by writing the collision operator in terms of the moments, and examine only the decay of the non-equilibrium moments. In such a representation, close to local equilibrium, all moments are small, except for the conserved moments, most notably the density. Because of the more complicated dependence of the equilibrium distribution on the momentum, the momentum dependence is more complicated than simply quadratic. So quadratic terms that connect a mode to the local density will dominate the collision operator, and since the density does not change during the collision, the linear relaxation of those moments is recovered.

A BGK collision operator implies that we have the matrix of Eq. (17)

[TABLE]

where $\tau^{\pi}$ is supposed to be a constant independent of the (non-conserved) moments. We can obtain this if we assume $\tilde{\pi}^{(0)}\ll\rho$ and make the usual assumption $u\ll 1$ to get

[TABLE]

The relaxation time can then be deduced from the collision term in moment form

[TABLE]

from which we get

[TABLE]

which becomes a constant, i.e. independent of $\tilde{\pi}^{(0)}$ and $u$, only in the limit mentioned above. This means that the range or reachable relaxation times are given by $\tau^{\pi}\in[1,\infty]$. For simple collisions the important over-relaxation regime for $\tau\in[1/2,1]$ that is often utilized in lattice Boltzmann simulations is unavailable for this lattice gas method. An extension of the collision process that includes over-relaxation is possible for high enough densities, but it derivation is outside the scope of this paper.

For systems far from equilibrium or high mach numbers these relaxation times become functions of the state variables in the lattice gas. The dependence of the relaxation time on $u$ violates Galilean invariance, and can in principle be remedied by setting

[TABLE]

for some constant $\tilde{\lambda}$, as can be seen in Eq. (69).

We test our analytical prediction of the decay of the $\pi$ mode by averaging over 1000 realizations of an initial configuration of $n_{i}(x,0)=\{100,100,100\}$ on a system with $L_{x}=1000$ and $\lambda=1$ and examine the decay of $\langle\pi\rangle$ as a function of the number of collisions. There is a small subtlety in relating the result of the predictions of the Boltzmann averaged lattice gas to actual fluctuating lattice gas simulations. Strictly speaking the analysis we have performed corresponds to relaxations of an ensemble average of systems where the both the local density and the local velocity correspond to the Boltzmann density $\rho$ and Boltzmann velocity $u$ only on average and averaging over the fluctuations is implied. The same is true for the equilibrium distribution: this is an equilibrium average where the realizations averaged over are fluctuating and only the expectation values correspond to the imposed $\rho$ and $u$.

To compare our Boltzmann theory to simulations we therefore have to consider the average relaxation of $\pi$ in a full simulation, including streaming, of a macroscopically homogeneous system. These results are shown in Figure 6, and show excellent agreement between simulation an theory.

The local equilibrium distribution, i.e. without taking the full Boltzmann ensemble average, has a density dependent second moment. This is easily seen as a single particle with mean velocity zero corresponds to the single state $(n_{-1}=0,n_{0}=1,n_{1}=0)$ which has $\pi=-\sqrt{3/2}$. So the local equilibrium distribution in moment representation is

[TABLE]

For two particles we would have two states $(n_{-1}=0,n_{0}=2,n_{1}=0)$ and $(n_{-1}=1,n_{0}=0,n_{1}=1)$ with probability 16/17 and 1/17 respectively according to Eq. (31). With this we obtain

[TABLE]

For larger number of particles and different velocities the detailed calculation of the average occupation numbers becomes more complicated, but the existence of a small negative contribution, which does remain of the order of $-1$ and does not increase with the particle number, remains as a part of the local equilibrium distribution.

Apart from the implementation details, there is an underlying physical reason for this. This second moment is related to a local temperature. It is well known that in a molecular system such a local temperature of a small domain depends on the discretization and at smaller discretization more of the kinetic energy resides in fluctuations of the local velocity and less in the temperature. In the extreme case of a single particle in a discretization cell the temperature would be zero, and all the energy would reside in the momentum of the cell. Therefore the global equilibrium distribution will have a definite momentum and a higher temperature, whereas a local equilibrium will have a fluctuating momentum and a lower temperature. Since the moments here are calculated with respect to the global equilibrium distribution we find that the local equilibrium distributions have small negative deviations of the second moment. This is why the local second moment does not decay towards zero in Figure 6, but will become slightly negative instead. In this Figure we see that the full solution of Eq. (69) gives a very good description of the non-linear decay for situations far from equilibrium.

IV.2 D2Q9 results

The derivation of the D2Q9 collision operator follows the D1Q3 derivation we just showed, but because we have 6 modes that relax, instead of just 1, collision rules become very lengthy to write down here. They are easily derived from the collision classes presented in Figure 3. We now have nine $\lambda_{i}$ which make the equivalent of Eq. (59) again rather lengthy to write. The equivalent of Eq. (64) shows a similar structure: many quadratic terms in the $m_{i}$. Of the $9^{3}=729$ possible combinations 243 actually occur, each typically associated with several $\lambda_{i}$. See Supplemental Material at [URL will be inserted by publisher] for a mathematica notebook that contains all the terms in question.

To calculate the equilibrium distribution we need to set $\Omega_{i}=0$. This initially appeared like a hard problem, and Mathemactia was unable to solve the resulting equations. However, we noticed that using the product of two one-dimensional equilibrium distribution functions

[TABLE]

turns out to be the simple solution (no summations over repeated Greek indices is implied here). We did not initially realize this, but at the 2017 DSFD conference in Erlangen AW had a chance to discuss these results with Santosh Ansumali. To both our surprise it turns out that the equilibrium distribution of entropic lattice Boltzmann derived in Ansumali et al. (2003) and Eq. (78) are identical, even if they are written in a very different form. In retrospect this result is not surprising since both methods will minimize the H-function Wagner (1998)

[TABLE]

To validate that our result for the equilibrium distribution agrees with results of the simulation method we ran a simulation on a lattice with dimensions $L_{x}=100$, $L_{y}=10$, $C=10$, $\lambda_{1}=15/128$, $\lambda_{2-7}=1/4$, $\lambda_{8}=1/8$, $\lambda_{9}=18/144$, $\langle\rho\rangle=36$ and $\langle U_{y}\rangle=0.25$. We varied $\langle U_{x}\rangle$. We then iterated over 20,000 iterations to ensure we are in equilibrium and then measured for another 10,000 iterations. The results of these simulations are shown in Fig. 7. As expected our values for the $f_{i}^{eq}$ agree exactly with the measured values. We also compare the quadratic equilibrium distribution of Eq. (9) to the values of the MCLG equilibrium distribution (78) or, equivalently, the entropic LB. They agree around $u_{x}^{eq}=0$ by design, but for larger velocities the two models begin to diverge. This is expected because as the mean velocity increases, the $f_{i}^{eq}$ of Eq. (78) eventually become negative, which the $n_{i}$ cannot do. Note that we have picked a very large $U_{x}=0.25$ to emphasize the differences between the entropic and polynomial equilibrium distributions.

As an example we show here the relaxation of the stress mode, which is responsible for recovering the shear viscosity. Similar to our analysis for the D1Q3 case we examine the effect of the collision on the deviation of the distribution from the local equilibrium distribution. We write the non-equilibrium part of our moments as

[TABLE]

where we introduced

[TABLE]

so that $dM^{a}=0$ in equilibrium. We write the effect of a single collision on these non-equilibrium moments as

[TABLE]

and we replace the appearances of $f_{i}$ in the collision term with

[TABLE]

The conservation of mass and momentum immediately imply that $\Omega^{0(n)}=\Omega^{1(n)}=\Omega^{2(n)}=0$. For the collision term related to the stress moment $\sigma_{xy}$ we get the rather lengthy

[TABLE]

This expression shows the coupling of the various modes. Unlike in the one-dimensional case we were unable to find an analytical solution for this non-linear relaxation. However, if we restrict our attention to moments that are close to local equilibrium, then the only terms that are not small are the conserved quantities. Here there are three conserved quantities, $\rho$, $j_{x}$, and $j_{y}$. And since $j_{x}=\rho u_{x}$ and $u_{x}<0.1$ for typical applications we can make the further assumption $j_{x}^{2}\ll\rho^{2}$. This means that the leading order terms are contained in the first 6 lines. The next order, small only because $|j|\ll\rho$, are the tems in the last four lines. The terms in lines 9 and 10 are smaller by one factor of $j_{x}/\rho\ll 1$. The remaining terms are quadratic in the non-equilibrium moments (remember those are zero in equilibrium) and can be neglected for systems close to local equilibrium. To leading order, neglecting $|j|/\rho$ terms, we also have $A_{x}=A_{y}=1$. The leading order for the relaxation of the $\sigma_{xy}$ mode is then given by

[TABLE]

Analogously to the one-dimensional case we can derive the effective relaxation time as

[TABLE]

Applying this approximation for all of the modes we obtain

[TABLE]

We see that the relaxation of all moments decouples, except for the $\sigma_{+}$ and $t$ moments. To recover a diagonal collision matrix in moment space $\Lambda^{ab}$ is equivalent to decoupling here. The 8th moment couples to $\sigma_{+}$ with a the same coefficient that the 5th moment couples to the eighths. We can set this coefficient to zero by demanding:

[TABLE]

which then decouples the modes and recovers the diagonal collision matrix for D2Q9. Additionally, to recover an isotropic stress tensor, we demand that the $\sigma_{xy}$ and $\sigma_{-}$ be relaxed at the same rate:

[TABLE]

The moments $q_{x}$ and $q_{y}$ are automatically relaxed at the same rate.

To derive the Boltzmann average of this MCLG method of Eq. (52) we first need to write the effect of $C$ collisions for Eq. (105). Once we ensured that the matrix terms decouple by imposing Eq. (106) we obtain for the collision matrix

[TABLE]

where the relaxation times are given by (95) for the $\sigma_{xy}$ moment and equivalent expressions for the relaxation of the other modes. This means we have derived a multi-relaxation-time BGK collision operator (with the assumption that we are close enough to equilibrium that the density mode is much larger than the non-conserved moments and that $\rho\gg|j|$) from the underlying lattice gas collision rules. We believe that this is the first such derivation. Earlier introductions of a BGK collision operator Higuera et al. (1989); Qian et al. (1992) were unable to derive this from a lattice gas, and instead did so as an ad-hoc assumption. The derivation is quite general: in moment space the dominant terms are those which are proportional to the density, and these leading order terms will contain only one power of one non-equilibrium moment, meaning that the collision operator can be written as a linear operator, i.e. a matrix collision applied to the non-equilibrium part of the particle distribution $f_{i}-f_{i}^{eq}$.

The Boltzmann average of the Monte Carlo Lattice Gas algorithm, with appropriate approximations discussed above, can now be written as

[TABLE]

where the collision matrix is given by Eq. (108) and the local equilibrium by Eq. (78).

Since this is exactly the Entropic lattice Boltzmann method (without over-relaxation for the moment) the standard derivations show that we recover the continuity

[TABLE]

and Navier-Stokes equation

[TABLE]

with

[TABLE]

To validate our prediction for the shear viscosity we set up a sinusoidal shear profile

[TABLE]

for which the Navier-Stokes equation has the analytical solution

[TABLE]

To simulate this system with the MCLG method we set up the lattice densities as 500 particles per lattice site, $L_{x}=100,L_{y}=101$, $\lambda$ values the same as those given below Eq. (79). To initialize the profile we set only the rest density $n_{0}$ and either $n_{1}$ or $n_{2}$ different from zero such that $\sum_{i}n_{i}=500$ and $\sum_{i}n_{i}v_{ix}=500\times 0.1\sin(2\pi y/L_{y})$, where the last equality is rounded to the nearest integer value. The densities are then equilibrated by performing 10,000 collisions on each lattice site (without streaming) to generate the initial configuration for the start of the simulation. The results of our MCLG simulations show a fluctuating version of the decay of the sinusoidal profile of such a shear wave. We average the results by averaging over the 100 $x$ values. To reduce fluctuations further we averaged this over $15$ separate realizations of this simulation. The simulation results show clean sinusoidal profiles with amplitudes that decay exponentially as a function of time. We obtain the amplitude from Eq. (121), given in an appendix. The results are shown in Figure 8, where we see that there is excellent agreement between the decay of the shear wave predicted by theory and the simulations, thus confirming our predictions for the shear viscosity in the MCLG method.

V Fluctuating properties of the MCLG

So far we have focused on the Boltzmann average of our lattice gas. The original reason we were interested in the lattice gas were its fluctuating properties. The fluctuations in an ideal gas have been discussed by Landau (Landau and Lifshitz, 1969, §114), where he shows that for a classical Boltzmann gas the number density in sub-volumes are Poisson distributed. The argument here is trivially extended to lattice gases showing that each density $n_{i}$ should be Poisson distributed:

[TABLE]

As a consequence we obtain

[TABLE]

In the lattice Boltzmann context trying to recover a Poisson distribution is problematic, since the densities $f_{i}$ are continuous, and there is no generally accepted extension of the Poisson statistic to continuous variables. Instead of trying to impose a Poisson distribution, Eq. (116) is used as the starting point for the derivation of fluctuating lattice Boltzmann methods. Using this second moment one finds that the lattice Boltzmann densities approximately represent a Poisson distribution, as seen in recent papers by Kaehler et al. Kaehler and Wagner (2013a) as well as by Wagner et al. Wagner and Strand (2016). However, these solutions are never exact and it is important to be aware of where the assumption of Poisson distributed $f_{i}$, which is used as an input to derive these methods, breaks down in the practical implementation.

It is important to note that this result is independent on the transport parameters, or in our case the number of collisions $C$, the various acceptance prefactors $\lambda_{i}$, the density or an imposed fluid velocity. We therefore tested the prediction of Poisson distributed occupation numbers $n_{i}$ in the MCLG method for a large variety of densities and imposed velocities and throughout found excellent agreement with the prediction from Eq. (115). An example is shown in Figure 9. Similarly we can test the independence of the $n_{i}$ by looking at

[TABLE]

Again we found excellent agreement for all densities and velocities that we tested. For the results shown here we used the same values for $\lambda_{i}$ as above, we initialized the lattice with a density of 360 and an initial velocity $u_{x}$ with $u_{y}=0$. We performed $C=10$ collisions per site per iteration. We then discarded the first 150,000 iterations to ensure we are looking at an equilibrium systems. We averaged over the next 500,000 steps on a lattice with $L_{x}=100$, $L_{y}=10$. We show the worst case result with respect to the results of Fig. 10 in Table 1. We see that we have excellent agreement with our prediction of Eq. (117).

To examine the agreement for different velocities without printing many of these tables we look at a measure of the deviation of the correlation matrix from the expected value

[TABLE]

A graph of this error function as a function of the mean velocity is shown in Figure 10. We see that there is no apparent systematic variation of this error function with $u$ even for extreme values for $u$ and with increasing averaging the error continues to diminish. This suggests that we have found a reference lattice gas implementation that exactly fulfills the predictions for a fluctuating ideal gas, a feat that has not been fully achieved to date for fluctuating lattice Boltzmann methods Kaehler and Wagner (2013a).

As a consequence of the correctly Poisson distributed fluctuations of the $n_{i}$ follows that all equal time corellators of composite quantities like the density or the momentum are correct and obey their respective Poisson or Skelam distributions.

VI Conclusions

In this article we have introduced an integer lattice gas method which employs a Monte Carlo collision operator. We show that the Boltzmann limit of this lattice gas recovers the entropic lattice Boltzmann method, which agrees with the usual lattice Boltzmann distribution with a polynomial distribution function for moderate velocities $|u|<0.1$. Remarkably we were able to derive a BGK collision operator directly from the lattice gas collision rules. Previous the BGK collision operator had been postulated only, but not derived. The equilibrium distribution recovered here is identical with that of entropic lattice Boltzmann Ansumali et al. (2003) and Boltzmann average of the lattice gas obeys the standard H-theorem.

Remarkably the fluctuating properties of this MCLG recover Poisson statistics exactly in a fully Galilean invariant manner. This suggests that an optimized version of this MCLG method may become a promising contender for the simulation of fluctuating fluids. Extending this result to non-ideal fluids will require a better understanding of fluctuations in these systems. We will examine what fluctuations for non-ideal systems with the help of the Molecular Dynamics lattice gas (MDLG) introduced by Parsa et al. Parsa and Wagner (2017). Developing an implementation that can compete with lattice Boltzmann approaches with regard to speed will require an implementation of a lattice gas collision operator that can perform multiple collisions at one time and a way to introducing over-relaxation that has been the way of lowering transport coefficients below the value obtained for full relaxation to local equilibrium. This is the subject of current research, and we anticipate to present solutions to these questions in subsequent publications.

Appendix A Amplitude of sine-wave decay

We can look at the problem of finding the amplitude $A(t)$ that corresponds to the best fit of a profile $u_{x}^{th}=A(t)\sin(2\pi y/L_{y})$ to some data $U_{x}(t,y)$. The mean-square deviation $E(A)$ would be given by

[TABLE]

and the requirement of this being minimal gives us

[TABLE]

so we get

[TABLE]

We can look at $E(A)/A^{2}$ as a measure for the quality of the fit.