# Lattice Gases and Cellular Automata

## 1 Historical Background

### 1.1 The Ising Model

The use of lattice gases for the study of equilibrium statistical mechanics dates back to a 1920 paper of Lenz  in which he proposed to model a ferromagnet by a regular $`D`$-dimensional lattice $`𝐋`$ of two-state “spins.” Physically, these may be thought of as the magnetization vectors of elemental magnetic domains, and the model constrains them to point in one of two directions, say “up” and “down.” For a two-dimensional lattice, this is illustrated schematically in Fig. 1. Mathematically, the state of the system can be described by the collection of variables $`S(𝐱)`$, indexed by the lattice points $`𝐱𝐋`$, and taking their values from the set $`\{1,+1\}`$; here $`S(𝐱)=+1`$ means that the spin at site $`𝐱`$ is pointing up, and $`S(𝐱)=1`$ means that it is pointing down. If we suppose that the lattice has a total of $`N|𝐋|`$ sites, then the total number of possible states of the system is $`2^N`$

To use these spins as a model of ferromagnetism, it was necessary to assign an energy to each of these $`2^N`$ states, in such a way as to make it energetically favorable for each spin to align with an externally applied magnetic field $`\alpha `$, and for neighboring spins to align with each other. The first of these goals is achieved by including an energy contribution $`\alpha S(𝐱)`$ for each spin present, and the second by including an energy contribution $`JS(𝐱)S(𝐲)`$ for each pair of neighboring sites $`𝐱`$ and $`𝐲`$. Thus, the full energy of the system is

$$H(𝐒)=\alpha \underset{𝐱}{}S(𝐱)\frac{J}{2}\underset{𝐱}{}\underset{𝐲𝐍(𝐱)}{}S(𝐱)S(𝐲),$$

where $`𝐍(𝐱)`$ denotes the set of sites neighboring site $`𝐱`$, and the factor of $`1/2`$ in front of the second term prevents double-counting of the pairs of spins.

To use this to study the equilibrium properties of a ferromagnet, it is necessary to compute the partition function

$$Z(K,h)\underset{N\mathrm{}}{lim}\underset{𝐒}{}\mathrm{exp}\left[\frac{H(𝐒)}{k_BT}\right],$$

where $`T`$ is the temperature, $`KJ/(k_BT)`$, $`h\alpha /(k_BT)`$, the sum over $`𝐒`$ includes all $`2^N`$ possible states of the system, and we have taken the thermodynamic limit by letting the number of spins go to infinity.

Lenz posed the problem of calculating this quantity to his student Ising, who solved it for a one-dimensional lattice of spins in 1925 . While Ising’s $`D=1`$ solution is elementary, Onsager’s $`D=2`$ solution for $`h=0`$ required almost another twenty years  to complete, and is significantly more complicated. The solution for the critical exponents for $`D=2`$ with $`h0`$ is a much more recent development, first published by Zamalodchikov  in 1989. The problem for $`D=3`$ is outstanding, even for $`h=0`$.

### 1.2 Universality and Materials Science

One might wonder why so much effort has been devoted to the Ising model when it is clearly only a crude idealization of a real ferromagnet. Certainly, nobody expects the detailed functional form of, say, the dependence of the Ising model’s magnetization

$$M(K,h)=\frac{_𝐱^NS(𝐱)\mathrm{exp}\left[\frac{H(𝐱)}{k_BT}\right]}{_𝐱^N\mathrm{exp}\left[\frac{H(𝐱)}{k_BT}\right]}=\frac{\mathrm{ln}Z(K,h)}{h}$$

on the temperature $`T`$ to be valid for any real material. There are, however, good reasons to believe that certain features of this functional form are universal – that is, model-independent. This is particularly true near criticality (in the $`D=2`$ and $`D=3`$ Ising models), where the spin-spin correlation length diverges, and fluctuations at all length scales are present. For example, at zero applied field and near criticality, the magnetization varies as

$$M=\{\begin{array}{cc}0\hfill & \text{for }T>T_c\hfill \\ M_0(\frac{T_cT}{T_c})^\beta \hfill & \text{for }TT_c,\hfill \end{array}$$

where $`T_c`$ is the critical temperature, $`M_0`$ is a proportionality constant, and $`\beta `$ is an example of what is called a critical exponent. The scale invariance of the fluctuations at the critical point allow a renormalization group treatment which indicates that the critical exponent should be rather insensitive to the particular model Hamiltonian used.

In fact, critical exponents should depend on only the dimensionality of the space and the symmetries of the underlying Hamiltonian function. For example, the unmagnetized Ising-model Hamiltonian is invariant under the symmetry group $`Z_2`$ – that is, multiplication in the set $`\{1,+1\}`$ – because the energy is invariant under the inversion of all the spins in the system. Systems with $`Z_2`$ symmetry are expected to have $`\beta =1/8`$ in $`D=2`$, and $`\beta 0.33`$ in $`D=3`$. A related lattice spin model, called the Heisenberg model, endows each spin with a vector orientation in three dimensions and has an interaction Hamiltonian that depends only on dot products of these vectors at neighboring sites. Since these are invariant under the continuous group of SO(3) rotations, we might expect a different critical exponent for $`\beta `$, and in fact this is the case: $`\beta 0.36`$ for the $`D=3`$ Heisenberg model.

Thus, universality teaches us that it is possible to learn some “real physics” by studying highly idealized models such as the Ising model. This realization led to a flurry of variants of lattice spin models, appropriate to various real materials. As an example, we consider an ingenious model developed by Widom  to describe microemulsions. A microemulsion is created by the addition of a surfactant or amphiphile to a mixture of two immiscible fluids, such as oil and water. An amphiphile is a chemical that typically has an ionic hydrophyllic end that likes to sit in water, covalently bonded to a hydrocarbon chain which is hydrophobic in that it prefers to sit in oil. This situation gives rise to two crucial properties: First, the free energy of the amphiphile is lowest when it lives on the interface between the two fluids. Second, the presence of the amphiphiles on the interface gives it a rigidity, or bending energy, that is proportional to the square of the local mean curvature <sup>1</sup><sup>1</sup>1More generally, it is proportional to the square of the difference of the local mean curvature and some spontaneous value thereof. Note that this rigidity is in addition to the surface tension that is always present between immiscible fluids.. Together, these two properties can give rise to some spectacular behavior. For example, thanks to the first property, when there is insufficient interface to accomodate the amphiphile, it becomes energetically profitable for the amphiphile to create new interface to inhabit. It does this by breaking up the separated mixture of water and oil into emulsion droplets. The droplets become smaller as more amphiphile is added, until they become so small (on the order of 50 nm) that the curvature energy mentioned in the second property removes the incentive for them to get any smaller. If the amount of surfactant continues to increase past this point a sponge phase can result, or, at lower temperatures, ordered lyotropic phases consisting of alternating sheets or tubes of oil and water. The self-organization of these structures as a result of relatively simple chemical properties is a methodology for nanofabrication, and scientific and industrial uses of these materials abound.

To model microemulsions using a lattice gas, Widom  introduced, in the mid 1980’s, a model very similar to that of Ising, in that each site on a Cartesian grid could be in one of two states, $`S(𝐱)\{1,+1\}`$. One clever innovation of this model is that it situates the particles on the links between the lattice vertices, rather than on the vertices themselves. Links can be characterized by the values of $`S`$ at the vertices they connect. Thus, each link can be in any one of four possible states: Links that connect two $`S=0`$ vertices are said to contain water. Links that connect two $`S=1`$ vertices are said to contain oil. Finally, links that begin at an $`S=1`$ ($`S=0`$) vertex and end at an $`S=0`$ ($`S=1`$) vertex are said to contain amphiphile whose hydrophyllic end points toward the end (beginning) of the link. The great advantage of this model is that the representation itself literally forces the amphiphile to live on the interface between the oil and water, and orients it correctly. This is illustrated in Fig. 2.

The Hamiltonian for Widom’s model includes two-spin interactions, as were present in the Ising model Hamiltonian, but it also includes three-spin interactions, necessary to capture the curvature energy described above and required for the rigidity of the interface. Since the number of particles of each type may change when a spin is flipped, it is also necessary to append terms of the form $`\mu _Wn_W+\mu _On_O+\mu _An_A`$, where the $`\mu `$’s are the (fixed) chemical potential per particle and the $`n`$’s are the particle numbers of each species, and run the Monte Carlo simulation in the grand canonical ensemble.

Widom’s model has been used to reproduce much of the above-described phenomenology associated with amphiphilic fluids. In particular, droplet, sponge and lyotropic phases are all seen, as are coexistence regions between these phases.

### 1.3 From Statics to Dynamics

The Ising model is discussed in many (perhaps most) statistical physics textbooks. Our purpose for bringing it up in this context is merely to note that it has two “false identities.” First, it looks a bit like a Hamiltonian dynamical system because it has an energy function (which we have obligingly denoted by $`H`$) that is exponentiated and integrated to get a partition function – just as is done with the Hamiltonian of a system of classical particles when deriving, e.g., the Mayer cluster expansion. Second, it looks a bit like a cellular automaton (CA) because its state can be represented by a single bit at each site. In fact, as defined above, the Ising model is neither of these things. Indeed, it is not a dynamical system at all. It is posed only as an equilibrium model. You can use it to find a partition function, and from that derive all of its thermodynamic properties, but thermodynamics (its name notwithstanding) is not dynamics at all. When it is used to study the difference between two states of a given system, it pays little or no attention to the particulars of the dynamical process that takes the system between those states. By contrast, Hamiltonian dynamical systems require a continuous phase space, endowed with a Poisson bracket that is antisymmetric and obeys the Jacobi identity; given such a bracket structure, any phase function $`H`$ defines a deterministic dynamics on the phase space. Likewise, a CA is a Markovian dynamical system on a discrete phase space, with transition rules that are applied to each site simultaneously and depend only upon the states of neighboring sites.

When the Ising model is simulated on a computer to study its equilibrium properties, it begins to look more like a dynamical system. The usual method is to perform a Monte Carlo simulation by flipping spins randomly, and then accepting or rejecting these flips according to the Metropolis algorithm , while sampling desired equilibrium properties. Because the Metropolis algorithm is Markovian in nature, this is a dynamical system in the mathematical sense of the word. However it is still neither Hamiltonian nor a CA. The discrete phase space would seem to distinguish it from most Hamiltonian systems <sup>2</sup><sup>2</sup>2There have been interesting attempts to construct fully discrete Hamiltonian systems, but such ideas are still in the embryonic stage of development., while the sequential nature of the algorithm (changes are considered to only one spin at a time) distinguishes it from the latter. More pertinently, it should be noted that the dynamics of the Metropolis algorithm do not necessarily have any relation to physical dynamics – other than the property of converging to the equilibrium state. The physical dynamical properties of a system cannot be studied by (ordinary) Monte Carlo methods; molecular dynamics, or some other such microcanonical algorithm must be used.

Given the utility of the Ising model in studying the equilibrium properties of materials, it should not be surprising that people have tried very hard to extend such models to include real dynamics. For a ferromagnet, this would mean being able to study the approach of the magnetization to equilibrium. One way to approach this problem is to invent a microcanonical dynamics for the Ising model. Such a model would evolve the system in such a way as to conserve energy globally and maintain the condition of detailed balance. If the probability that the system is in state $`𝐒`$ is denoted by $`P(𝐒)`$, and if the probability of transition of the system from state $`𝐒`$ to state $`𝐒^{}`$ is denoted by $`A(𝐒𝐒^{})`$, then the latter condition, which ensures that the dynamical process will converge to a Boltzmann-Gibbs equilibrium, may be written

$$P(𝐒)A(𝐒𝐒^{})=P(𝐒^{})A(𝐒^{}𝐒).$$

Such a microcanonical dynamics for the Ising model was developed by Creutz in 1983 . To localize the energy change incurred by the flip of a spin, it was necessary to update the lattice in a “checkerboard” pattern – first updating the black squares and then the red ones. This has the effect of “freezing” a site’s neighbors while it is flipped, so that the energy change incurred by the flip can be calculated in advance. Because the algorithm is to be microcanonical, this energy must be stored somewhere. Creutz solved this problem by allowing for each site to have a “bank” where it can make local deposits and withdrawals of energy. The bank must have a finite capacity though, and if flipping the spin at a site would result in an overflow or underflow of the local bank, the flip is not done; otherwise, it is done.

Creutz demonstrated that his dynamics equilibrate the $`D=2`$ Ising model. The lattice is initialized with a total energy (particles plus banks) that never changes in the course of the simulation. The portion of the energy in the spins themselves, however, may change, and must therefore be monitored by the simulation. The important point, however, is that Creutz’ model is able to make a CA out of the Ising model. The checkerboard updating is not a real problem, since it can be incorporated in the CA framework by including a “time parity” bit at each site, initialized in a checkerboard pattern, and mandated to toggle at each time step. The value of this bit might then be used to determine which sites will attempt to flip, and the rest is naught but fully deterministic and reversible nearest-neighbor communication and simultaneous updating at each site – in other words, a reversible CA.

### 1.4 From Ferromagnetism to Hydrodynamics

The approach to equilibrium of a ferromagnet is certainly interesting, but, because the order parameter (magnetization) is not a conserved quantity, it happens relatively fast. The path to equilibrium in systems with conserved order parameters is generally more tortuous and difficult – and hence interesting – because such systems are more constrained. Hence, for example, a viscous fluid at high Reynolds number achieves equilibrium only after turbulent relaxation in which intricate structures may be spawned across a wide range of scales in length and time.

To be sure, it is possible to conserve the ferromagnetic order parameter. The so-called Kawasaki dynamics  of the Ising model flips only pairs of oppositely directed spins, so that the total magnetization is conserved. If this idea were combined with Creutz’ microcanonical dynamics, it would be possible to perform a microcanonical simulation of a system with a conserved order parameter, but I am not aware that this has been tried. The more natural setting for systems with conserved order parameters is hydrodynamics. Fluids generally have a conserved mass and momentum, and if compressibility effects are considered, the conservation of energy will also play an important role.

There was, however, a conceptual problem with the CA simulation of fluids: Using lattice gases to model ferromagnets is intuitive enough because it is not difficult to picture a regular array of magnetic domains. It seemed much less clear how to model a fluid on a lattice, however, and almost a half a century elapsed between the development of the Ising model and the first lattice-gas models of hydrodynamics.

## 2 Hydrodynamic Lattice Gases

### 2.1 The Kadanoff-Swift Model

The goal of a hydrodynamic lattice gas is to take the same “minimalist” approach to fluids that the Ising model takes to ferromagnets. The object is not a precise model of the dynamics at the finest scales, but rather to invent a fictitious microdynamics whose coarse-grained behavior – in the thermodynamic limit – lies in the same universality class as the phenomenon under study. Along these lines, the approach that seems the most promising is to model the fluid at the level of fictitious “molecules” that can move about and collide, as they do in a real fluid, conserving mass, momentum and (for compressible fluids) energy as they do so.

The first attempt (known to me) along these lines was undertaken by Kadanoff and Swift in 1968 . They considered a $`D=2`$ Cartesian lattice, each site of which may be occupied by a particle, as shown in Fig. 3. Each particle is tagged with one of four momenta, oriented along the diagonals,

$`𝐜_1`$ $`=`$ $`+\widehat{𝐱}+\widehat{𝐲}`$   
$`𝐜_2`$ $`=`$ $`\widehat{𝐱}+\widehat{𝐲}`$   
$`𝐜_3`$ $`=`$ $`\widehat{𝐱}\widehat{𝐲}`$   
$`𝐜_4`$ $`=`$ $`+\widehat{𝐱}\widehat{𝐲},`$   
where $`\widehat{𝐱}`$ and $`\widehat{𝐲}`$ are the unit vectors in the $`x`$ and $`y`$ directions. An Ising-like Hamiltonian is defined. Then, at each step, a particle is randomly selected and one of three things happen to it:

* It moves in the direction of its velocity vector to the next site where it could land without violating conservation of energy.
* It diffuses, to any empty neighboring site, carrying its momentum with it.
* It exchanges momentum vectors with another neighboring particle.

The dynamics thus defined conserves mass, momentum and energy, and obeys the principle of detailed balance. Note that it is not technically a cellular automaton, because of the sequential nature of the particle updates, but it might be made into one by using some generalization of checkerboard updating and/or Cruetz’ “banks.” To my knowledge, however, this has never been tried.

The Kadanoff-Swift (KS) model exhibits many features of real fluids, such as sound-wave propagation, and long-time tails in velocity autocorrelation functions. As the authors noted, however, it does not faithfully reproduce the correct equations of motion of a viscous (or, for that matter, any other kind of) fluid. In particular, the model exhibits a strong lattice anisotropy; the decay of sound waves, for example, depends on their direction of propagation with respect to the underlying lattice.

### 2.2 The HPP Model and Kinetic Theory

The next advance in the lattice modelling of fluids came in the mid 1970’s, when Hardy, de Pazzis and Pomeau introduced a new lattice model with a number of innovations that warrant discussion here. The HPP model, named for its authors, also resides on a $`D=2`$ Cartesian lattice. Particle velocities are taken from the set

$`𝐜_1`$ $`=`$ $`+\widehat{𝐱}`$   
$`𝐜_2`$ $`=`$ $`+\widehat{𝐲}`$   
$`𝐜_3`$ $`=`$ $`\widehat{𝐱}`$   
$`𝐜_4`$ $`=`$ $`\widehat{𝐲},`$   
and there may be anywhere from zero to four particles at each site. The only restriction is that there may not be more than one particle with a particular velocity, at a particular site, at a particular time. This “exclusion principle” makes it possible to represent the state of any site $`𝐱`$ by four bits; bit $`n_j(𝐱,t)`$, where $`j\{1,2,3,4\}`$, encodes the presence ($`1`$) or absence ($`0`$) of a particle with velocity $`𝐜_j`$ at site $`𝐱`$ and time $`t`$.

Given this representation, the dynamics is defined as follows: At each time step, each site may experience a purely local collision, in which its particles rearrange their velocity vectors in such a way as to conserve mass and momentum. A moment’s thought indicates that the only nontrivial way for this to happen is if exactly two particles enter a site from opposite directions, and exit in the other two opposite directions; for any other configuration, the particles simply retain their incoming velocities, “passing through” each other without interacting. After the collisions, all particles “stream” to the site in the direction of their velocity vector. Note that all sites experience collisions simultaneously, and all particles stream simultaneously as well. The HPP model can thus be regarded as a CA.

The above described dynamical rule can actually be expressed algebraically as follows:

$$n_j(𝐱+𝐜_j,t+\mathrm{\Delta }t)=n_j(𝐱,t)+\omega _j,$$
(1)   
where $`\mathrm{\Delta }t`$ is the time associated with one step of the dynamical process (equal to unity if natural “lattice units” are adopted), and where $`\omega _j`$ is called the collision operator. If the collision operator were not present, this equation would simply state that a particle with velocity $`𝐜_j`$ will exist at site $`𝐱+𝐜_j`$ at time $`t+\mathrm{\Delta }t`$ if it existed at site $`𝐱`$ one time step earlier. This captures the streaming process. To include the collisions, $`\omega _j`$ must subtract (add) a particle from direction $`j`$ if the incoming state will undergo a nontrivial collision that will deplete (augment) that direction. For the HPP rule described above, we have

$`\omega _1`$ $`=`$ $`n_1(1n_2)n_3(1n_4)+(1n_1)n_2(1n_3)n_4`$   
$`\omega _2`$ $`=`$ $`n_2(1n_3)n_4(1n_1)+(1n_2)n_3(1n_4)n_1`$   
$`\omega _3`$ $`=`$ $`n_3(1n_4)n_1(1n_2)+(1n_3)n_4(1n_1)n_2`$   
$`\omega _4`$ $`=`$ $`n_4(1n_1)n_2(1n_3)+(1n_4)n_1(1n_2)n_3.`$   
Here we have made use of the fact that multiplication of bits is equivalent to the logical “and” operation, substraction from one is equivalent to the logical “not” operation, and addition of mutually exclusive bits is equivalent to the logical “or” operation. Thus, for example, the quantity $`(1n_1)n_2(1n_3)n_4`$ is equal to one if directions $`2`$ and $`4`$ are occupied, and $`1`$ and $`3`$ are not. Thus, this term appears with a plus sign in $`\omega _1`$ and $`\omega _3`$, since the resulting collision augments those directions, and with a minus sign in $`\omega _2`$ and $`\omega _4`$ since those directions are depleted.

This algebraic description of the exact microscopic motion of the particles is somewhat akin to the Klimontovich description of continuum kinetic theory . It is made somewhat simpler by the discreteness of the spatial lattice, and the finite number of allowed velocities. To use this microscopic description to find the fluid equations obeyed by the coarse-grained density and hydrodynamic velocity requires all of the tricks of kinetic theory. The first step is to determine a kinetic equation for the single-particle distribution function. We define this by an ensemble average, supposing that we have a large number of such lattices, with initial conditions sampled from some (unspecified) distribution, and writing

$$N_j(𝐱,t)n_j(𝐱,t),$$

where the angle brackets denote the ensemble average. Note that, whereas $`n_j`$ is bit-valued, $`N_j`$ is real-valued. By taking the ensemble average of Eq. (1), we arrive at

$`N_j(𝐱+𝐜_j,t+\mathrm{\Delta }t)`$ $`=`$ $`N_j(𝐱,t)+\omega _j`$   
$`=`$ $`N_j(𝐱,t)n_j(1n_{j+1})n_{j+2}(1n_{j+3})`$   
$`+(1n_j)n_{j+1}(1n_{j+2})n_{j+3},`$   
where all subscripts are evaluated modulo $`4`$. At this point, we see that we have a problem. The collision operator is nonlinear in the $`n_j`$’s, and the average of a product is not generally equal to the product of the averages – not if the quantities involved are correlated. Thus, the dynamical equation of the $`N_j`$’s will involve averages of products, such as $`n_jn_{j+2}`$. To know these, it would be necessary to write kinetic equations for these two-point correlations, but these will involve still higher correlations, etc. This infinite series of equations is the lattice-gas analog of the BBGKY hierarchy of kinetic theory.

To truncate this hierarchy and obtain a closed equation for the $`N_j`$’s, it is necessary to make a physical approximation: We shall assume that the particles entering a collision are uncorrelated. This approximation is tantamount to Boltzmann’s famous molecular chaos assumption. It is unlikely to be true, especially for high densities and low-dimensional lattices. For now we note that, under this assumption, it is possible to replace the average of products above by the product of averages, resulting in a closed equation for the single-particle distribution function,

$`N_j(𝐱+𝐜_j,t+\mathrm{\Delta }t)`$ $`=`$ $`N_j(𝐱,t)+\mathrm{\Omega }_j`$ (3)   
$`=`$ $`N_j(𝐱,t)N_j(1N_{j+1})N_{j+2}(1N_{j+3})`$   
$`+(1N_j)N_{j+1}(1N_{j+2})N_{j+3}.`$   
This is called the lattice-Boltzmann equation. From this, it is possible, using the Chapman-Enskog analysis of classical kinetic theory, to derive the hydrodynamic equations obeyed by the mass density,

$$\rho (𝐱,t)\underset{i}{}N_i(𝐱,t)$$
(4)   
and the momentum density,

$$\rho (𝐱,t)𝐮(𝐱,t)\underset{i}{}𝐜_iN_i(𝐱,t).$$
(5)   
Note that Eqs. (4) and (5) are the discrete-velocity analog of the usual integration over velocity space to obtain the hydrodynamic densities. In this way, the fully compressible hydrodynamic equations for the HPP lattice gas were worked out in the original papers by HPP . For our purposes, we note that the result of this exercise in the incompressible limit is $`\mathbf{}𝐮=0`$ and

$`{\displaystyle \frac{u_x}{t}}+g(\rho ){\displaystyle \frac{}{x}}u_y^2`$ $`=`$ $`{\displaystyle \frac{1}{\rho }}{\displaystyle \frac{P}{x}}+\nu (\rho ){\displaystyle \frac{^2u_x}{x^2}}\left(\nu (\rho )+{\displaystyle \frac{1}{2}}\right){\displaystyle \frac{^2u_y}{xy}}`$   
$`{\displaystyle \frac{u_y}{t}}+g(\rho ){\displaystyle \frac{}{y}}u_x^2`$ $`=`$ $`{\displaystyle \frac{1}{\rho }}{\displaystyle \frac{P}{y}}+\nu (\rho ){\displaystyle \frac{^2u_y}{y^2}}\left(\nu (\rho )+{\displaystyle \frac{1}{2}}\right){\displaystyle \frac{^2u_x}{xy}},`$ (6)   
where $`P`$ is the pressure, and we have defined the functions

$$g(\rho )\frac{1\rho /2}{1\rho /4}$$

and

$$\nu (\rho )\frac{1}{2\rho (1\rho /4)}\frac{1}{2}.$$

Eqs. (6) has some superficial resemblence to the Navier-Stokes equations of viscous fluid dynamics, but closer inspection reveals some important differences. Like the KS model studied earlier, the HPP model gives rise to anisotropic hydrodynamic equations that are not invariant under a global spatial rotation. They involve derivatives with respect to $`x`$ and $`y`$ in combinations that cannot be expressed in terms of the familiar $`\mathbf{}`$ operator of vector calculus. Rather, the grid coordinates $`x`$ and $`y`$ have a preferred status. Hence, for example, the drag of a KS or HPP fluid as it flows past a generic obstacle will depend on that obstacle’s angle of orientation with respect to the underlying lattice. At the time, this was not considered a problem, since the real purpose of the KS and HPP models was to study the statistical physics of fluids, and both models could do this well. Traditional computational fluid dynamicists, however, were not inclined to take notice of this as a serious numerical method unless and until a way was found to remove the unphysical anisotropy.

### 2.3 The FHP Model

Another thirteen years passed from the introduction of the HPP model to the solution of the anisotropy problem in 1986 by Frisch, Hasslacher and Pomeau , and simultaneously by Wolfram . The FHP lattice gas, named after the authors of the first reference given above, is very similar to that of HPP in that the evolution proceeds by alternating collision and streaming steps – hence, it is again a cellular automaton. The only real difference is that it is based on a triangular lattice instead of a Cartesian one, as shown in Fig. 4. Now one would expect that a six-fold symmetric lattice would give rise to a more isotropic model than a four-fold symmetric one. The surprising result of the 1986 studies, however, is that the six-fold version does not only improve the isotropy – it yields perfect isotropy!

To see why isotropy is recovered on a triangular lattice, we generalize the results of the Chapman-Enskog analysis mentioned above. It turns out that the most general form for the viscous term in the hydrodynamic equations above is

$$\frac{u_i}{t}+\mathrm{}=c_{ijkl}_j_ku_l,$$

where $`c_{ijkl}`$ is a fourth-rank tensor that is constructed from the lattice vectors, and hence shares all of their symmetries. If the lattice is invariant under rotation by $`60^{}`$ or $`90^{}`$, so then will be the components of this tensor. Note that if this tensor is isotropic – so that it is expressible in terms of the Kronecker delta – then we have most generally

$$c_{ijkl}=\nu \delta _{jk}\delta _{il}+\mu \delta _{ij}\delta _{kl}+\mu \delta _{ik}\delta _{jl},$$

whence our hydrodynamic equation becomes

$$\frac{𝐮}{t}+\mathrm{}=\nu ^2𝐮+2\mu \mathbf{}\left(\mathbf{}𝐮\right).$$

This form, expressible with vector notation, is indeed isotropic and the two terms on the right may be identified with the shear and bulk viscosity terms of the Navier-Stokes equation, respectively. Thus, a sufficient condition for the attainment of isotropy is to ensure that the $`c_{ijkl}`$ tensor is isotropic. It turns out that the only rank-four tensors in two dimensions that are invariant under $`60^{}`$ rotation are isotropic; by contrast, there exist rank-four tensors in two dimensions that are invariant under $`90^{}`$ rotation but which are not isotropic (cannot be expressed in terms of the Kronecker delta). Thus, generically, isotropy will be recovered on a triangular grid and not on a Cartesian one.

In the incompressible limit, the result of the Chapman-Enskog analysis for the FHP fluid is then $`\mathbf{}𝐮=0`$ and

$$\frac{𝐮}{t}+g(\rho )𝐮\mathbf{}𝐮=\mathbf{}P+\nu (\rho )^2𝐮,$$
(7)   
where $`g(\rho )`$ is defined as it was for the HPP model, and $`\nu (\rho )`$ is a well-behaved function of density whose precise form depends on the particulars of the collision rules (e.g., whether or not three- or four-body collisions are considered). Thus, isotropy is recovered, and there is no problem writing Eq. (7) in the usual notation of vector calculus, as we have done.

There is, however, one last lingering problem: The factor of $`g(\rho )`$ in front of the inertial term (second term on the left-hand side) is equal to unity in the real Navier-Stokes equations. This is a consequence of Galilean invariance (GI) which, like isotropy, is yet another Lie-group symmetry possessed by the real Navier-Stokes equations, but which we may not take for granted in our lattice model. To see that the convective derivative operator

$$\frac{D}{Dt}\frac{}{t}+𝐮\mathbf{}$$
(8)   
has GI, we consider the transformation of coordinates

$`𝐱^{}`$ $`=`$ $`𝐱𝐕t`$   
$`t^{}`$ $`=`$ $`t,`$   
where $`𝐕`$ is a constant vector. This has inverse

$`𝐱`$ $`=`$ $`𝐱^{}+𝐕t^{}`$   
$`t`$ $`=`$ $`t^{}.`$   
Then the derivatives transform as follows:

$`{\displaystyle \frac{}{t}}`$ $`=`$ $`{\displaystyle \frac{t^{}}{t}}{\displaystyle \frac{}{t^{}}}+{\displaystyle \frac{𝐱^{}}{t}}\mathbf{}^{}={\displaystyle \frac{}{t^{}}}𝐕\mathbf{}^{}`$   
$`\mathbf{}`$ $`=`$ $`\left(\mathbf{}t^{}\right){\displaystyle \frac{}{t^{}}}+\left(\mathbf{}𝐱^{}\right)\mathbf{}^{}=\mathbf{}^{},`$   
and the hydrodynamic velocity transforms as follows

$$𝐮=\frac{d𝐱}{dt}=\frac{d}{dt^{}}\left(𝐱^{}+𝐕t^{}\right)=\frac{d𝐱^{}}{dt^{}}+𝐕=𝐮^{}+𝐕.$$

Under this Galilean transformation, we see that

$$\frac{}{t}+𝐮\mathbf{}=\frac{}{t^{}}𝐕\mathbf{}^{}+\left(𝐮^{}+𝐕\right)\mathbf{}^{}=\frac{}{t^{}}+𝐮^{}\mathbf{}^{},$$

so the form of the convective derivative operator is indeed Galilean invariant. The presence of the $`g(\rho )`$ factor invalidates this argument. The lack of GI in the FHP fluid is due to the fact that the lattice itself constitutes a preferred reference frame.

We can now see the isotropy and GI problems in a unified light: Lie-group symmetries are often responsible for particular features of hydrodynamic equations; for example, isotropy implies that the equation can be written in the notation of vector calculus, and GI implies that the convective derivative will have the form of Eq. (8). Such symmetries may be broken by the presence of the fixed lattice, and when this happens the corresponding features of the hydrodynamic equations may be destroyed. We managed to recover isotropy by noting that the maximum rank of the tensors in our hydrodynamic equation is four, and that rank-four tensors with $`60^{}`$ rotational invariance are always isotropic. How do we fix the GI problem? There are several ways to do so. For incompressible flow, for which $`\rho `$ and hence $`g(\rho )`$ are constant, the easiest solution is to scale $`𝐮`$ and the pressure $`P`$ by $`g(\rho )`$; since the term in which $`g(\rho )`$ appears is the only one that is quadratic in $`𝐮`$, this works handily. Another approach introduces additional particle velocities to force $`g(\rho )`$ to unity.

Another consideration is the functional form of the shear viscosity, $`\nu (\rho )`$. For lattice gases that satisfy detailed balance, this quantity will have a minimum as a function of $`\rho `$. This minimum value of the viscosity sets the maximum value of the Reynolds number that may be simulated on a given size lattice.

The result of this 1986 work was a cellular automaton model whose coarse-grained behavior was that of a Navier-Stokes fluid. Now the computational fluid dynamicists took notice, and there was a five-year-long flurry of activity in the field. Some of the accomplishments during this time were

1. the extension of the model to three dimensions.
2. simulations of flow past various obstacles, and comparisons with other more conventional methods of computational fluid dynamics (CFD).
3. clever variations on the collision rules of the basic FHP model intended to achieve lower viscosity minima, and experimental tests of the functional form of $`\nu (\rho )`$.
4. careful measurements of long-time tails in velocity autocorrelation functions, and finite-size corrections to the viscosity.
5. extensions of the model to simulate complex fluid hydrodynamics, including interfaces and surface tension.
6. a host of algorithmic tricks, and even special-purpose hardware, to simulate such hydrodynamic lattice gases on parallel computers.

The number of papers produced in this period is far too great to review here; the interested reader is referred to the excellent text by Rothman and Zaleski  and the secondary references therein, for a review of this work. Items 1 and 5 on the above list will be discussed briefly below in Subsection 2.4 and Section 3, respectively.

One reason that hydrodynamic lattice gases captured peoples’ imagination at this time was that they represented an altogether new way of doing CFD – and, indeed, computational physics in general. Conventional CFD methods began with the Navier-Stokes equations, and discretized them in one of a variety of ways. Lattice gases, by contrast, defined a kind of particle kinetics from which the Navier-Stokes equations were emergent – just as they are emergent for a real fluid. There are definite advantages to such a “physical” approach, aside from its undeniable aesthetic appeal. For example, one often overlooked advantage of lattice-gas models is their unconditional stability. The Navier-Stokes equations have a basis in kinetic theory, as the behavior of a system of particles whose collisions conserve mass and momentum. The fact that these underlying collisions obey a detailed-balance condition ensures the validity of the $`H`$ theorem, the fluctuation-dissipation theorem, Onsager reciprocity, and a host of other critically important properties with macroscopic consequences. When its kinetic origins are cavalierly ignored, and the Navier-Stokes equations are “chopped up” into finite-difference schemes, these important properties can be lost. The discretized evolution equations need no longer satisfy an $`H`$ theorem, and the notion of underlying fluctuations may lose meaning altogether. As the first practioneers of finite-difference simulations on digital computers found in the 1940’s and 1950’s, the result can be the development and growth of high-wavenumber numerical instabilities, and indeed these have plagued essentially all CFD methodologies in all of the decades since. Such instabilities are entirely unphysical because they represent a clear violation of the $`H`$ theorem; indeed, the Second Law of Thermodynamics would preclude their occurrence. Numerical analysts have responded to this problem with textbooks full of ways to “patch up” these anomalies – including upwind differencing, artificial viscosity, and a host of other very clever tricks – but from a physicist’s point of view it would have been much better if the original discretization process had retained more of the underlying physics, so that these problems had not occurred in the first place. Lattice gases represented an important first step in this direction. As was shown shortly after their first applications to hydrodynamics, they can be constructed with an $`H`$ theorem that rigorously precludes any kind of numerical instability. More glibly stated, lattice gases avoid numerical instabilities in precisely the same way that Nature herself does so.

Even the computer implementation of a hydrodynamic lattice gas was novel. All other CFD methodologies make use of floating-point numbers to represent real quantities. In such floating-point numbers, some bits represent the mantissas, others the exponents and others the signs. A lattice gas, by contrast, can use a representation of one bit per particle. All bits thus play an equal role in some sense; this idea is sometimes colorfully referred to as “bit democracy.”

In spite of these successes, the use of lattice gases for the simulation of simple (single-phase) Navier-Stokes fluids declined substantially in the early 1990’s. Ironically, it was largely supplanted by the direct floating-point computer simulation of the lattice Boltzmann equation, Eq. (3). The lattice Boltzmann equation also has kinetic underpinnings in a sense, though its representation is that of a single-particle distribution function, rather than a full particle-level description. Most importantly, the lattice Boltzmann framework allows for greater accuracy than LGA by effectively eliminating kinetic fluctuations. Measured quantities are therefore much less noisy, and require less ensemble averaging to compute accurately. For certain complex fluids such as microemulsions, however, lattice gases remain a very effective simulation methodology.

### 2.4 The FCHC Model

To find a lattice-gas model of isotropic three dimensional hydrodynamics, it was necessary to find a $`D=3`$ lattice, under whose symmetry group the only isotropic rank-four tensors are invariant. The trouble is that no such regular lattice exists in three dimensions. In 1987, however, Frisch et al.  noticed that a lattice with the required symmetry does exist in four dimensions. It is called the face-centered hypercubic (FCHC) lattice, and is self-dual with 24 lattice vectors per site. The 24 lattice vectors can be arranged in three groups of eight, such that any two groups of eight comprise the 16 vertices of a regular four-dimensional hypercube, and the vectors of the third group of eight point in the direction of the centers of the faces of that hypercube (hence, the name FCHC).

The 24 lattice vectors are most easily enumerated by taking all integer quadruples $`(i,j,k,l)`$ that lie at a distance $`\sqrt{2}`$ from the origin. Clearly, two of the integers components must be zero, and the other two must be $`\pm 1`$. There are $`\left(\genfrac{}{}{0pt}{}{4}{2}\right)=6`$ ways to choose which two are zero, and $`2^2=4`$ ways of assigning $`\pm 1`$ to the other two, for a total of 24.

It was noticed  that by projecting this lattice back to three dimensions – by, say, ignoring the fourth coordinate – a simple set of $`D=3`$ lattice vectors was obtained which worked. The three dimensional lattice thus obtained is not a regular lattice, in that not all the lattice vectors have the same length, and some of them have multiplicity two (because they correspond to two different vectors on the four-dimensional lattice), but it works nevertheless. The projected lattice vectors are thus

| $`j`$ | $`c_{jx}`$ | $`c_{jy}`$ | $`c_{jz}`$ |
| --- | --- | --- | --- |
| 0,1 | +1 | 0 | 0 |
| 2,3 | -1 | 0 | 0 |
| 4 | 0 | +1 | +1 |
| 5 | 0 | +1 | -1 |
| 6 | 0 | -1 | +1 |
| 7 | 0 | -1 | -1 |
| 8,9 | 0 | +1 | 0 |
| 10,11 | 0 | -1 | 0 |
| 12 | +1 | 0 | +1 |
| 13 | +1 | 0 | -1 |
| 14 | -1 | 0 | +1 |
| 15 | -1 | 0 | -1 |
| 16,17 | 0 | 0 | +1 |
| 18,19 | 0 | 0 | -1 |
| 20 | +1 | +1 | 0 |
| 21 | +1 | -1 | 0 |
| 22 | -1 | +1 | 0 |
| 23 | -1 | -1 | 0 |

Note that four of them have multiplicity two. Also note that we have listed them in the three groups of eight, described above. A more geometric description of the FCHC lattice and its use for lattice-gas simulations is given in the reference by Adler et al. .

## 3 Applications to Complex Fluids

The simulation of the hydrodynamics of complex fluids, such as immiscible flow, coexisting phases, emulsions, colloids, liquid crystals, gels and foams is one of the principal outstanding challenges of computational condensed matter physics. Hydrodynamic equations for such materials are often not known or ill-posed, so that finite-difference discretizations are not even an option. Molecular dynamics can, of course, be employed, but it advances in time steps that are typically between $`10^2`$ and $`10^3`$ of a mean-free path. In this context, lattice gases offer the possibility of an inexpensive molecular dynamics – one for which particles are still discrete, but can advance in steps on the order of a mean-free time (since the particles typically suffer a collision at each time step). Earlier in this review, we mentioned that lattice models can be useful for studying the equilibrium properties of complex fluids, and we considered the model of Widom in some detail. We finish this survey by demonstrating how lattice gas hydrodynamics have begun to explore the hydrodynamics of such materials.

In 1988, Rothman and Keller  introduced a lattice-gas model of immiscible flow, such as that of oil and water. They accomplished this by tagging the lattice gas particles with two “colors,” to distinguish oil and water. Their collisions were required to conserve the total mass of each color separately, as well as the total momentum. They then skewed the collision outcomes to favor those that send particles towards sites dominated by other particles of the same color. This affinity of particles for other particles of the same color gives the two phases cohesion, and the interface surface tension. Let us call the two colors “red” and “blue,” and use $`n_i^R(𝐱,t)`$ and $`n_i^B(𝐱,t)`$, respectively, to denote the occupation number of each color in velocity direction $`i`$, at site $`𝐱`$ and at time $`t`$. Then Rothman and Keller defined a color field

$$𝐄\underset{i}{}𝐜_i\left[\underset{j}{}\left(n_j^Rn_j^B\right)\right]$$

at each site that points in the direction of increasing (red minus blue) color, and a color flux

$$𝐉\underset{i}{}𝐜_i\left(n_i^Rn_i^B\right)$$

for each possible outgoing state. They then chose the outcome that minimized the color work,

$$H=𝐉𝐄.$$
(9)   
More generally, Chen et al.  pointed out that one should assign probabilities to each of the possible outcomes based on the Boltzmann weights $`\mathrm{exp}\left(\frac{H}{k_BT}\right)`$, where $`T`$ is a temperature. This model has been studied extensively (see ), and it exhibits phase separation and surface tension.

Another rather different model of two coexisting phases, such as water and water vapor, was worked out by Appert and Zaleski in 1990 . The model gives particles an attraction by allowing mass- and momentum-conserving collisions between certain configurations of particles at different sites. This model has also been extensively studied (again, see ). Both this model and that of Rothman and Keller have been cast in the lattice Boltzmann framework as well ; indeed, this remains a very active area of study.

Neither the Rothman-Keller and Appert-Zaleski models, as originally posed, obeyed the principle of detailed balance. Since this is important, for all of the reasons described above, some effort has been devoted to restoring detailed balance to hydrodynamic lattice-gas models with interaction between particles at different sites. This turns out not to be easy. The best attempts to date have taken as a starting point a 1988 method due to Colvin et al. , called Maximally Discretized Molecular Dynamics or “$`(MD)^2`$.” This is essentially a return to the Kadanoff-Swift representation with a maximum of one particle per site, and sequential propagation of particles. (Thus, this model is no longer, strictly speaking, a cellular automaton.) Like the KS model, particles are allowed to exchange momentum with their neighbors, and they can step in the direction of their momentum to the nearest empty site. A hard-sphere interaction is introduced that may extend over more than one site. Mass, momentum and energy are exactly conserved, and detailed balance is maintained.

Arbitrary interaction potentials – possibly beyond just the nearest neighbor in range – were added to the $`(MD)^2`$ model by Gunn et al.  in 1993, and this model remains near the state of the art <sup>3</sup><sup>3</sup>3The algorithm in the paper by Gunn et al. used a continuum-valued velocity, but this is not an essential feature.. The collisions in this model do not change the energy, and may therefore be unconditionally accepted. The propagation step, however, may change the energy, and the algorithm can be implemented in either a microcanonical version (in which the propagation is accepted only if the energy change is zero), or a canonical version (in which it is accepted according to the Metropolis algorithm).

Finally, we note that many variants of these algorithms exist and are useful for various purposes. Boghosian et al.  developed a variation of the Rothman-Keller model that allows for the inclusion of a surfactant phase, allowing the simulation of microemulsions. To see how this works, first note that the quantity

$$\underset{i}{}\left(n_i^Rn_i^B\right)$$

can be thought of as a color charge in the Rothman-Keller model, in that its current is the color flux, and its vector-weighted sum over neighbors is the color field. In this context, surfactant particles are introduced as color dipoles, and numerous terms are added to the Hamiltonian, Eq. (9), to account for the color-dipole interaction that makes the surfactant prefer to live on the interface, and the dipole-dipole interaction that gives rise to the curvature energy. Details are given in the reference .

The $`D=2`$ model has been studied in some detail, and preliminary $`D=3`$ results have been obtained as of this writing. The model is able to track the formation and saturation of droplet (Fig. 5, $`D=2`$), wormlike-micelle (Fig. 6, $`D=3`$), sponge (Fig. 7, $`D=2`$) and lamellar (Fig. 8, $`D=3`$) phases, and the time dependence of this saturation has been studied . Interfacial fluctuations in the presence of surfactant have been studied with this model , and it has also been used for the first simulations of the shear-induced sponge-to-lamellar phase transition  (Fig. 9, $`D=2`$). Note that all of these applications involve nonequilibrium or dynamical processes that have previously been difficult to address. The most lengthy MD simulations to date are barely able to see the self-assembly of a single emulsion droplet. The lattice-gas method, by contrast, is able to study the growth and saturation of many such droplets and larger structures, as can be seen in these figures.

## 4 Conclusion

We have reviewed the use of lattice gases for physical problems, with emphasis on their applications to hydrodynamics. We have traced the evolution of hydrodynamic lattice gases in the 1960’s and 1970’s, leading up to the development of the FHP model in 1986. We have then seen some of the attempts to add interactions – possibly between particles at different sites – to such models, in order to simulate complex-fluid hydrodynamics; we have also discussed the loss of detailed balance in such models, and its recovery by the $`(MD)^2`$ algorithm and its variants. Though much of the CFD-related activity in this field has migrated to lattice Boltzmann methods, it is this author’s belief that lattice gases remain very well suited for the role of “inexpensive molecular dynamics” for the hydrodynamic simulation of complex fluids.

In addition, it seems likely that this kind of simulation can still benefit from the development of special-purpose hardware. Interest in this area has been somewhat stifled to date because the ongoing meteoric rise in workstation performance (per unit cost) has made it dangerous to try to develop any kind of hardware outside of the mainstream. Sooner or later, however, workstation performance will begin to saturate, and the very real differences between hardware optimized for lattice gases and for general-purpose floating-point calculations will once again become a target on which hardware designers might well focus some activity.

## Acknowledgements

The author would like to thank Raissa d’Souza, Norman Margolus and Peter Coveney for their review of the manuscript and helpful comments. The three-dimensional images were made by Andrew Emerton, Peter Love and David Bailey of Oxford University.