# BU-CCS-991001 A Particulate Basis for an Immiscible Lattice-Gas Model

## 1 Introduction

In 1986 it was discovered that certain mass and momentum conserving lattice-gas automata gave rise to the isotropic Navier-Stokes equations in the hydrodynamic limit . In 1988, Rothman and Keller extended this discovery by introducing a hydrodynamic lattice-gas model of immiscible fluids . Their model, and lattice Boltzmann variants thereof, have become an important tool for simulating the hydrodynamics of multiphase flow . In both the original and RK lattice-gas models, the dynamics can be decomposed into two steps: In the first, the particles propagate along the lattice vectors to new sites; in the second, the particles entering each site collide by redistributing mass and momentum.

In the Rothman-Keller (RK) model, the masses of the various immiscible fluid species and the total momentum are conserved locally, but the choice of collision outcome at each site of the RK model depends on the water-minus-oil order parameter, or “color,” of the neighboring sites. See Fig. 1 for one set of possible collision outcomes that might be allowed by the conservation laws. The flux of the color is determined for each such outgoing state, and its local gradient or field is determined by examining the neighboring sites. The negative of the dot product of this color flux and color field is then a measure of the propensity of outgoing particles to move to sites dominated by particles of their own type, and was called the color work by RK; their prescription was then to choose the outcome that minimizes this work in order to create cohesion and interfacial tension. (In case of a tie, the outcome is chosen randomly from among the states with minimal color work.) In later work, Chen, Chen, Doolen and Lee , and Chan and Liang  noted that this minimization of color work is really just the low-temperature limit of a Boltzmann sampling procedure.

This paper proposes a microscopic interpretation of the RK model. In the limit where the ratio of the mean-free path to the interaction range is small, we show that an arbitrary interaction potential can be expanded in terms of dot products of fluxes and fields, of which the RK model is merely the first term. This observation accounts for much of the success and utility of the RK model in describing the hydrodynamics of multiphase flow.

## 2 Hydrodynamic Lattice-Gas Automata

In lattice-gas models of hydrodynamics, incoming particles collide at each site $`𝐱`$ of a lattice $``$, in a manner to be discussed at length shortly, and then propagate to one of $`B`$ neighboring sites $`𝐱+𝐜_i`$, where $`i\{1,\mathrm{},B\}`$. (Note that some of the $`𝐜_i`$’s may be zero in order to accommodate “rest particles” in the model.) We suppose that the occupancy of each of these $`B`$ channels can be represented by $`L`$ bits $`n_i^{\mathrm{}}(𝐱)\{0,1\}`$, where $`i\{1,\mathrm{},B\}`$ and $`\mathrm{}\{1,\mathrm{},L\}`$. Throughout the remainder of this paper, we shall illustrate various concepts by applying them to three concrete examples:

* Example 1: In a lattice gas for a single-species Navier-Stokes fluid , we take one bit ($`L=1`$) in each direction that represents the presence or absence of a particle moving in that direction.
* Example 2: In a lattice gas for two immiscible fluids , on the other hand, we might take two bits per direction ($`L=2`$) so that there can be one bit for water particles $`n_i^W(𝐱)`$ and one bit for oil particles $`n_i^O(𝐱)`$ in each direction.
* Example 3: We consider the model of Chen, Chen, Doolen and Lee , in which there is one bit in each direction ($`L=1`$), one “rest” particle direction $`i=R`$ such that $`𝐜_R=0`$, and only the rest particles feel an interaction potential.

We suppose that there is a charge-like attribute $`q_i(𝐱)`$ associated with the bits in direction $`i`$ at site $`𝐱`$, and we specialize to potential energies of the form

$$V=\frac{1}{2}\underset{𝐱,𝐲}{}\underset{i,j}{}q_i(𝐱)q_j(𝐱+𝐲)\varphi (|𝐲|),$$
(1)   
where the factor of $`1/2`$ prevents double counting. In Example 1, the charge-like attribute might be equal to the channel occupancy $`n_i(𝐱)`$. In Example 2, on the other hand, the charge-like attribute might be the order parameter

$$q_i(𝐱)=n_i^W(𝐱)n_i^O(𝐱)$$
(2)   
which measures the excess of water over oil. In Example 3, we take $`q_i(𝐱)=n_i(𝐱)\delta _{iR}`$, since only the rest particles have a charge-like attribute. In what follows, we shall also have occasion to refer to the total (summed over directions) charge-like attribute of a site $`q(𝐱)_iq_i(𝐱)`$.

As noted above, any fluid model has some set of quantities that must be conserved in collisions. In Examples 1 and 3, we should demand conservation of the mass $`_in_i(𝐱)`$ and the momentum $`_in_i(𝐱)𝐜_i`$. In Example 2, on the other hand, collisions must conserve water mass $`_in_i^W(𝐱)`$, oil mass $`_in_i^O(𝐱)`$ and total momentum $`_i[n_i^W(𝐱)+n_i^O(𝐱)]𝐜_i`$. Alternatively, we can say that all three examples conserve total mass and the total momentum, while Example 2 also conserves the charge-like attribute, $`q(𝐱)`$. These conserved quantities naturally partition the set of all states of a given site(s) into equivalence classes of states with the same values for all of the conserved quantities. For example, Fig. 1, relevant to Example 2 above, illustrates the equivalence class comprised of the six possible collisional outcomes that may result when one water particle and one oil particle enter a single site on a two-dimensional triangular lattice from opposite directions (and, hence, with zero total momentum).

## 3 Collisional Energetics

We denote the postcollision charge-like attribute with velocity $`𝐜_i`$ at site $`𝐱`$ by $`q_i^{}(𝐱)`$. Upon subsequent propagation, the charge $`q_i^{}(𝐱)`$ will be at position $`𝐱+𝐜_i`$, and the charge $`q_j^{}(𝐱+𝐲)`$ will be at position $`𝐱+𝐲+𝐜_j`$. This is illustrated in Fig. 2. The change in the potential energy due to both collision and propagation is then given by

$`\mathrm{\Delta }V`$ $`=`$ $`{\displaystyle \underset{𝐱,𝐲}{}}{\displaystyle \underset{i,j}{}}q_i^{}(𝐱)q_j^{}(𝐱+𝐲)\varphi \left(\left|𝐲+𝐜_j𝐜_i\right|\right){\displaystyle \underset{𝐱,𝐲}{}}{\displaystyle \underset{i,j}{}}q_i(𝐱)q_j(𝐱+𝐲)\varphi \left(\left|𝐲\right|\right)`$ (3)   
$`=`$ $`{\displaystyle \underset{𝐱,𝐲}{}}{\displaystyle \underset{i,j}{}}q_i^{}(𝐱)q_j^{}(𝐱+𝐲)\left[\varphi \left(\left|𝐲+𝐜_j𝐜_i\right|\right)\varphi \left(\left|𝐲\right|\right)\right]`$   
$`+{\displaystyle \underset{𝐱,𝐲}{}}{\displaystyle \underset{i,j}{}}\left[q_i^{}(𝐱)q_j^{}(𝐱+𝐲)q_i(𝐱)q_j(𝐱+𝐲)\right]\varphi \left(\left|𝐲\right|\right)`$   
$`=`$ $`\mathrm{\Delta }V_c+\mathrm{\Delta }V_n,`$   
where we have defined the contribution to $`\mathrm{\Delta }V`$ due to the movement of the interacting particles,

$$\mathrm{\Delta }V_c\underset{𝐱,𝐲}{}\underset{i,j}{}q_i^{}(𝐱)q_j^{}(𝐱+𝐲)\left[\varphi \left(\left|𝐲+𝐜_j𝐜_i\right|\right)\varphi \left(\left|𝐲\right|\right)\right],$$
(4)   
and that due to nonconservation of the charge-like attribute,

$$\mathrm{\Delta }V_n\underset{𝐱,𝐲}{}\left[q^{}(𝐱)q^{}(𝐱+𝐲)q(𝐱)q(𝐱+𝐲)\right]\varphi \left(\left|𝐲\right|\right).$$
(5)   
Note that $`\mathrm{\Delta }V_c`$ vanishes for systems in which the interacting particles do not move, including our Example 3. Likewise, note that $`\mathrm{\Delta }V_n`$ vanishes for systems with a conserved charge-like attribute, including our Examples 1 and 2.

## 4 The Flux-Field Decomposition

If the potential kernel $`\varphi (y)`$ is analytic, we may expand $`\mathrm{\Delta }V_c`$ in a Taylor series in the ratio of the characteristic lattice spacing<sup>1</sup><sup>1</sup>1Since lattice gases are usually dense fluids for which particles undergo a collision at every step, $`c`$ can also be thought of as a mean-free path. $`c`$ to the characteristic interaction range $`y`$. This is done in Appendix A where it is shown that every term of this series can be expressed as the complete inner product of a tensor flux with a tensor field. More specifically, we find that

$$\mathrm{\Delta }V_c=\underset{𝐱}{}\underset{n=1}{\overset{\mathrm{}}{}}\underset{r=n/2}{\overset{n}{}}𝒥_r^{}(𝐱)\stackrel{r}{}_{n,r}(𝐱),$$
(6)   
where we have defined the (local) $`r`$th-rank tensor flux of outgoing particles,

$$𝒥_r^{}(𝐱)\underset{i}{}q_i^{}(𝐱)\left(\stackrel{r}{}𝐜_i\right),$$
(7)   
and where the (generally nonlocal) $`r`$th-rank tensor fields are defined in terms of $`(nr)`$th-rank outgoing tensor fluxes of neighboring sites,

$`_{n,r}^{}(𝐱)`$ $``$ $`{\displaystyle \frac{(1)^r}{(1+\delta _{n,2r})n!}}\left({\displaystyle \genfrac{}{}{0pt}{}{n}{r}}\right){\displaystyle \underset{𝐲}{}}𝒦_n(𝐲){\displaystyle \stackrel{nr}{}}𝒥_{nr}^{}(𝐱+𝐲).`$ (8)   
In all the above expressions, primes are used to denote dependence on postcollision values, $`^r`$ denotes an $`r`$-fold outer product, $`^r`$ denotes an $`r`$-fold inner product, and we have defined the $`n`$th-rank completely symmetric kernel

$$𝒦_n(𝐲)\underset{m=n/2}{\overset{n}{}}\frac{\varphi _m(y)}{(nm)!}\text{per}\left[\left(\stackrel{2mn}{}𝐲\right)\left(\stackrel{nm}{}\mathrm{𝟏}\right)\right],$$
(9)   
where “per” indicates a summation over all distinct permutations of indices, and where we have defined the following functions related to the derivatives of $`\varphi (y)`$:

$$\varphi _m(y)\left(\frac{1}{y}\frac{d}{dy}\right)^m\varphi (y).$$
(10)   
The first few such kernels are thus given by

$`\left[𝒦_1(𝐲)\right]_i`$ $`=`$ $`\varphi _1(y)y_i`$ (11)   
$`\left[𝒦_2(𝐲)\right]_{ij}`$ $`=`$ $`\varphi _1(y)\delta _{ij}+\varphi _2(y)y_iy_j`$ (12)   
$`\left[𝒦_3(𝐲)\right]_{ijk}`$ $`=`$ $`\varphi _2(y)\left(y_i\delta _{jk}+y_j\delta _{ik}+y_k\delta _{ij}\right)+\varphi _3(y)y_iy_jy_k`$ (13)   
$`\left[𝒦_4(𝐲)\right]_{ijkl}`$ $`=`$ $`\varphi _2(y)\left(\delta _{ij}\delta _{kl}+\delta _{ik}\delta _{jl}+\delta _{il}\delta _{jk}\right)/2+`$ (14)   
$`\varphi _3(y)\left(y_iy_j\delta _{kl}+y_iy_k\delta _{jl}+y_iy_l\delta _{jk}+y_jy_k\delta _{il}+y_jy_l\delta _{ik}+y_ky_l\delta _{ij}\right)+`$   
$`\varphi _4(y)y_iy_jy_ky_l.`$   
For completeness, we note that, since the zero-rank flux is just the charge, the portion of $`\mathrm{\Delta }V`$ arising from nonconservation of the charge-like attribute may also be written in terms of these fluxes,

$$\mathrm{\Delta }V_n=\frac{1}{2}\underset{𝐱,𝐲}{}\left[𝒥_0^{}(𝐱)𝒥_0^{}(𝐱+𝐲)𝒥_0(𝐱)𝒥_0(𝐱+𝐲)\right]\varphi _0(y)$$
(15)   
as follows immediately from Eqs. (5) and (7).

Because the fields themselves depend on the post-collision fluxes at neighboring sites, it is generally necessary to include $`\mathrm{\Delta }V_n`$ in the collisional energetics. There are, however, some very useful exceptions to this rule. Some of the fluxes $`𝒥_r(𝐱)`$ may be conserved, in which case $`𝒥_r^{}(𝐱)=𝒥_r(𝐱)`$, since their values for incoming and outgoing states must be identical. If all of the fluxes that go into the calculation of a field $`_{n,r}^{}(𝐱)`$ are conserved quantities, then we can write $`_{n,r}^{}(𝐱)=_{n,r}(𝐱)`$ as well.

For example, let us return to consider our Example 2. At zeroth order, Eq. (5) indicates that $`\mathrm{\Delta }V_n`$ vanishes because $`𝒥_0`$ is conserved. At first order ($`n=1`$) the field $`_{1,1}^{}(𝐱)`$ depends only on the zero-order flux – namely the total order parameter (water minus oil) – at a given site, and this is conserved. Hence, the first-order energy change can be computed for the order-parameter flux for each outgoing state, and the field based on the incoming states of the neighbors $`𝐱+𝐲`$. If we restrict this set of neighbors to immediately adjacent sites, and look no further than lowest order in the Taylor expansion, we see that this reduces precisely to the RK model. The RK model is thus an approximation that is valid to first-order in $`c/y`$. While this fact may explain much of the utility of the RK model, in the following section we shall see that it may also indicate that the scaling limit of the RK model is more subtle than previously suspected. In any case, this also means that an alternative model that uses the exact potential energy, $`\mathrm{\Delta }V`$ of Eq. (3), to sample post-collision states at each site will certainly be no worse than the RK model which is known to capture much of the phenomenology of immiscible fluid dynamics.

Note that if both the flux and field in one of the terms of the Taylor expansion can be evaluated based on incoming quantities, then that term of $`\mathrm{\Delta }V`$ is the same for all outgoing quantities, and therefore does not discriminate between them, so that it is necessary to go to higher order (at least until one encounters the first nonconserved fluxes) to obtain any interaction at all. We saw this with the vanishing of $`\mathrm{\Delta }V_n`$ for systems with conserved charge-like attribute, such as Example 2 above. To see this happen at higher order, let us consider the simpler Example 1 which, perhaps for this very reason, has been less studied. Both fluxes $`𝒥_0`$ and $`𝒥_1`$ are conserved, since they are the mass and momentum, respectively. It follows that the entire $`n=1`$ term of $`\mathrm{\Delta }V_c`$ is the same for all outgoing states. Nontrivial interaction between particles therefore does not even begin until the $`n=2`$ term of the expansion, since $`𝒥_2`$ (which is, in fact, related to the pressure tensor) is not a conserved quantity.

## 5 The Hydrodynamic Limit of the RK Model

In order to obtain useful quantitative information from a hydrodynamic lattice gas, one must be careful to work in the correct asymptotic regime. This usually involves scaling the various dimensionless parameters of the problem with the Knudsen number $`\text{Kn}\lambda /L`$, where $`\lambda `$ is the mean-free path and $`L`$ is the characteristic size. In incompressible Navier-Stokes flow, for example, one desires that the Mach number $`\text{M}U/C`$, where $`U`$ is the characteristic hydrodynamic velocity and $`C`$ is the sound speed, scale with the Knudsen number $`\text{M}𝒪(\text{Kn})`$. Since the viscosity $`\nu `$ goes as the product of mean-free path and sound speed, $`\lambda C`$, this implies that the Reynolds number Re scales as $`\text{M}/\text{Kn}𝒪(1)`$<sup>2</sup><sup>2</sup>2Note that this does not mean that the numerical value of the Reynolds number must be near unity. Rather it means that Re approaches a constant value in the scaling limit.. We also demand that the Strouhal number $`\text{St}U\tau /L`$ and the fractional density fluctuation $`\delta \rho /\rho _0`$, where $`\tau `$ is the mean-free time and $`\rho _0`$ is the average background density, both scale as $`𝒪(\text{Kn}^2)`$. This limit is well known  to reduce the compressible Navier-Stokes equations to their incompressible counterparts. Since, for a dense LGA, the mean-free path $`\lambda `$ goes as the grid size $`c`$, this means that in order to approach the continuum limit, every time one doubles the size of the lattice (halves Kn), one must quadruple the number of time steps (since St is quartered), and verify that the fractional density fluctuation (a measured “output” quantity in a lattice-gas simulation) is also quartered. Only when this scaling is verified can one be sure that one is working in the correct asymptotic regime.

The presence of an interparticle potential adds an additional length scale – the range $`y`$ of the force – and therefore a new dimensionless parameter $`\lambda /y`$, or equivalently $`c/y`$. To derive the flux-field decomposition, we demanded that this ratio be small, but of order unity; that is, we did not scale this parameter with the Knudsen number. Operationally, this means that every time the size of the lattice is doubled, the range of the force in lattice units should be kept the same. If it is ten lattice units at one resolution, it should be ten lattice units at all resolutions<sup>3</sup><sup>3</sup>3In some sense, this means that lattice artifacts never completely disappear, as the set of sites within a ten-lattice-unit radius of a given site are not distributed uniformly or isotropically. If this is deemed problematic, it may be possible to scale $`c/y`$ with $`\text{Kn}^{1/2}`$, or in some other way such that both $`c/y`$ and $`y/L`$ vanish in the scaling limit; this has the attraction of completely removing such problems in the scaling limit, but further exploration of such considerations lies outside the scope of this paper..

While Eq. (4) for $`\mathrm{\Delta }V_c`$ is exact, the flux-field decomposition of Eq. (6) is usually used to approximate $`\mathrm{\Delta }V_c`$ only to some specified order in $`c/y`$. Having determined that the RK model is just such an approximation, we are now in a position to examine how the error incurred by this approximation scales in the continuum limit. Let us compare the RK model to a variant of our model, in which we adopt the strategy of permitting the number of terms $`n_{\mathrm{max}}`$ that we retain in the Taylor expansion for $`\mathrm{\Delta }V_c`$ to increase in the hydrodynamic limit. We shall justify this strategy a posteriori. Specifically, let us take $`\xi `$ more terms each time the lattice size $`N`$ is doubled. If we take the system size in physical units $`L`$ to be fixed, then the lattice spacing is $`cL/N`$. For a dense lattice gas, the mean-free path is of order $`c`$, so the Knudsen number Kn scales as $`c/L`$. It follows that

$$n_{\mathrm{max}}=\xi \mathrm{log}_2\left(\frac{N}{N_0}\right)=n_0\xi \mathrm{log}_2\left(\text{Kn}\right),$$
(16)   
where $`N_0`$ and $`n_0`$ are constants. If we then take the interaction range in lattice units $`y/c`$ to be fixed, then the error term $`\epsilon `$ in the Taylor expansion goes as

$$\epsilon \left(\frac{c}{y}\right)^{n_{\mathrm{max}}}\left(\frac{c}{y}\right)^{n_0}\left(\frac{y}{c}\right)^{\xi \mathrm{log}_2(\text{Kn})}\left(\frac{c}{y}\right)^{n_0}\text{Kn}^{\xi \mathrm{log}_2(\frac{y}{c})}.$$
(17)   
For $`\xi =1`$ (adding one more term to the series at each lattice refinement), it follows that by keeping the lattice spacing less than half of the characteristic interaction range, the error will scale as the Knudsen number; likewise, by making the lattice spacing less than a fourth of the characteristic interaction range, the error will scale as the square of the Knudsen number; and so on. One can also increase $`\xi `$ to raise the power of the Knudsen number to which the error scales. Thus, for small but finite values of $`c/y`$ (order unity in the scaling limit), the error can be made to scale subdominantly to terms that are usually neglected in a Chapman-Enskog expansion anyway, providing the a posteriori justification promised above.

A potential problem with the RK model is that it does not refine the definition of the energy in this way at each level of the scaling limit ($`\xi =0`$), and so the corrections that it neglects may indeed matter in that limit. Of course, in most situations, one will choose a particular fixed value of $`n_{\mathrm{max}}`$; in fact, all studies of the RK model to date have used $`n_{\mathrm{max}}=1`$, simply because it has not been previously recognized that the RK model is only the first-order approximation to a more exact model. The RK model is known to exhibit certain anomalous phenomena; for example spurious currents are known to develop near interfaces, even if there is no bulk flow, and the surface tension is known to be slightly anisotropic . It is possible that such anomalies would be eliminated by the more exact treatment of the scaling limit advocated here, but a numerical test of this conjecture is outside the scope of the present work.

Of course, such anomalies may be regarded as tolerable as long as one appreciates that one is working only to lowest order of an asymptotic series. Indeed, because the series is asymptotic, there is no point in being overzealous about the value of $`n_{\mathrm{max}}`$, since the series may begin to diverge at some point. There is at least one situation, however, in which the observation that it is necessary to let $`n_{\mathrm{max}}`$ scale with Knudsen number may be critically important, and that is when one is studying the scaling of (possibly divergent) quantities with system size. For example, Rothman and Flekkøy  recently studied the scaling properties of fluctuating interfaces using the RK model, measuring, among other things, the saturated width of the interface as a function of system size. Superimposed upon the usual power-law behavior of the saturated width, they found a logarithmic correction that resisted theoretical explanation. We suggest that this anomaly may be due to scaling the system size with fixed $`n_{\mathrm{max}}`$ (since the RK model effectively fixes $`n_{\mathrm{max}}`$ at unity), rather than letting $`n_{\mathrm{max}}`$ increase linearly with system size as advocated here. Again, it would be interesting to verify this conjecture by redoing the numerical experiment using our model, but that is outside the scope of the present work.

The main point of this section has been to demonstrate that the flux-field decomposition allows one to understand in what sense the RK model is an approximation to an exact interaction potential, why it may be used for trial state sampling, how it might be corrected at higher order in $`c/y`$, and why its scaling limit may be more subtle than previously suspected. We emphasize that the strategy of scaling $`n_{\mathrm{max}}`$ with Kn was invoked only to facilitate this discussion of the scaling limit, and we are certainly not suggesting that it be used in practical simulations. If one is planning to work at any order of $`c/y`$ at which $`\mathrm{\Delta }V_c`$ involves nonlocal interactions, it makes much more sense to use Eq. (4) directly than to try to deal with higher order terms of the flux-field decomposition. For Monte Carlo sampling of the outgoing states, it is possible that the lowest-order local term of the flux-field decomposition could be used for sampling, while the exact expression for $`\mathrm{\Delta }V_c`$ could be used for the acceptance criterion.

## 6 Conclusions

We have shown that the lattice gas model developed by Rothman and Keller for immiscible-fluid hydrodynamics can be derived from an underlying model of particle interactions. From the enhanced understanding provided by our observation, we elucidated the nature of the hydrodynamic limit of the Rothman-Keller model, demonstrating that it is more subtle than previously suspected. Though practical simulations of the particulate model are likely to be significantly more compute-intensive than the original version of the Rothman-Keller model, this work is offerred in the spirit that it is always useful to know the exact model corresponding to any given approximation. We hope that this work helps to provide some theoretical basis for these models’ success, and perhaps for their ultimate improvement.

## Acknowlegements

BMB would like to thank the International Centre for Theoretical Physics (ICTP) for their hospitality during a portion of this work. BMB was supported in part by the United States Air Force Office of Scientific Research under grant number F49620-95-1-0285. The collaboration of BMB and PVC was facilitated by NATO grant number CRG 950356.

## Appendix A Derivation of Flux-Field Decomposition

To derive Eq. (6) for $`\mathrm{\Delta }V_c`$, we begin with the Taylor expansion of a function of the magnitude of a displaced vector,

$$\varphi (|𝐲+𝐜|)=\underset{n=0}{\overset{\mathrm{}}{}}\frac{ϵ^n}{n!}\underset{m=n/2}{\overset{n}{}}\frac{n!\varphi _m(y)}{2^{nm}(2mn)!(nm)!}\left(𝐲𝐜\right)^{2mn}\left(𝐜𝐜\right)^{nm},$$
(18)   
where $`ϵ`$ has been introduced as an expansion parameter to keep track of the order in $`c/y`$ (it is numerically equal to one), $`y|𝐲|`$, and we have defined the following functions related to the derivatives of $`\varphi (y)`$:

$$\varphi _m(y)\left(\frac{1}{y}\frac{d}{dy}\right)^m\varphi (y).$$
(19)   
We let $`𝐜𝐜_j𝐜_i`$ and use the binomial theorem to write

$`\left(𝐲𝐜\right)^{2mn}`$ $`=`$ $`\left(𝐲𝐜_j𝐲𝐜_i\right)^{2mn}`$ (20)   
$`=`$ $`{\displaystyle \underset{l=0}{\overset{2mn}{}}}{\displaystyle \frac{(2mn)!}{l!(2mnl)!}}\left(𝐲𝐜_j\right)^{2mnl}\left(𝐲𝐜_i\right)^l,`$   
and

$`\left(𝐜𝐜\right)^{nm}`$ $`=`$ $`\left(|𝐜_j|^2+|𝐜_i|^22𝐜_j𝐜_i\right)^{nm}`$ (21)   
$`=`$ $`{\displaystyle \underset{k=0}{\overset{nm}{}}}{\displaystyle \underset{p=0}{\overset{nmk}{}}}{\displaystyle \frac{(nm)!|𝐜_i|^{2k}|𝐜_j|^{2p}}{k!p!(nmkp)!}}\left(2𝐜_i𝐜_j\right)^{nmkp}.`$   
Inserting these into Eq. (4) for $`\mathrm{\Delta }V_c`$, we get

$`\mathrm{\Delta }V_c=`$ (22)   
$`{\displaystyle \frac{1}{2}}{\displaystyle \underset{𝐱,𝐲}{}}{\displaystyle \underset{i,j}{}}q_i^{}(𝐱)q_j^{}(𝐱+𝐲){\displaystyle \underset{n=1}{\overset{\mathrm{}}{}}}{\displaystyle \frac{ϵ^n}{n!}}{\displaystyle \underset{m=n/2}{\overset{n}{}}}`$   
$`{\displaystyle \underset{l=0}{\overset{2mn}{}}}{\displaystyle \underset{k=0}{\overset{nm}{}}}{\displaystyle \underset{p=0}{\overset{nmk}{}}}{\displaystyle \frac{(1)^{nm+lkp}n!}{2^{k+p}l!k!p!(2mnl)!(nmkp)!}}`$   
$`\varphi _m(y)\left(𝐲𝐜_j\right)^{2mnl}\left(𝐲𝐜_i\right)^l|𝐜_i|^{2k}|𝐜_j|^{2p}\left(𝐜_i𝐜_j\right)^{nmkp}.`$   
Eliminating $`p`$ in favor of the new summation index $`rnm+l+kp`$, and reordering the summations, this becomes

$`\mathrm{\Delta }V_c=`$ (23)   
$`{\displaystyle \frac{1}{2}}{\displaystyle \underset{n=1}{\overset{\mathrm{}}{}}}{\displaystyle \frac{ϵ^n}{n!}}{\displaystyle \underset{r=0}{\overset{n}{}}}{\displaystyle \underset{𝐱,𝐲}{}}{\displaystyle \underset{i,j}{}}q_i^{}(𝐱)q_j^{}(𝐱+𝐲){\displaystyle \underset{m=n/2}{\overset{n}{}}}{\displaystyle \underset{l=0}{\overset{2mn}{}}}{\displaystyle \underset{k=\mathrm{max}(0,mn+rl)}{\overset{\mathrm{min}(nm,(rl)/2)}{}}}`$   
$`{\displaystyle \frac{(1)^rn!\varphi _m(y)|𝐜_i|^{2k}|𝐜_j|^{2n2r2m+2l+2k}}{2^{nrm+l+2k}l!k!(nrm+l+k)!(2mnl)!(rl2k)!}}`$   
$`\left(𝐲𝐜_j\right)^{2mnl}\left(𝐲𝐜_i\right)^l\left(𝐜_i𝐜_j\right)^{rl2k},`$   
where we have adopted the convention that a sum is zero if its upper limit is less than its lower limit. By reinterpreting dot products raised to the $`s`$ power as the $`s`$-fold inner product of two $`s`$-fold outer products, we can rewrite this as follows

$`\mathrm{\Delta }V_c`$ $`=`$ $`{\displaystyle \frac{1}{2}}{\displaystyle \underset{n=1}{\overset{\mathrm{}}{}}}{\displaystyle \frac{ϵ^n}{n!}}{\displaystyle \underset{r=0}{\overset{n}{}}}(1)^r\left({\displaystyle \genfrac{}{}{0pt}{}{n}{r}}\right){\displaystyle \underset{𝐱,𝐲}{}}{\displaystyle \underset{i,j}{}}q_i^{}(𝐱)q_j^{}(𝐱+𝐲)`$ (24)   
$`\left[\left({\displaystyle \stackrel{r}{}}𝐜_i\right){\displaystyle \stackrel{r}{}}𝒦_n(𝐲){\displaystyle \stackrel{nr}{}}\left({\displaystyle \stackrel{nr}{}}𝐜_j\right)\right],`$   
where $`^r`$ denotes an $`r`$-fold outer product and $`^r`$ denotes an $`r`$-fold inner product, and where we have defined the kernel

$`𝒦_n(𝐲)`$ $``$ $`{\displaystyle \underset{m=n/2}{\overset{n}{}}}{\displaystyle \underset{l=0}{\overset{2mn}{}}}{\displaystyle \underset{k=\mathrm{max}(0,mn+rl)}{\overset{\mathrm{min}(nm,(rl)/2)}{}}}`$ (25)   
$`{\displaystyle \frac{r!(nr)!\varphi _m(y)}{2^{nrm+l+2k}l!k!(nrm+l+k)!(2mnl)!(rl2k)!}}`$   
$`\left[\left({\displaystyle \stackrel{2mn}{}}𝐲\right)\left({\displaystyle \stackrel{nm}{}}\mathrm{𝟏}\right)\right]`$   
$`=`$ $`{\displaystyle \underset{m=n/2}{\overset{n}{}}}{\displaystyle \frac{n!\varphi _m(y)}{2^{nm}(2mn)!(nm)!}}\left[\left({\displaystyle \stackrel{2mn}{}}𝐲\right)\left({\displaystyle \stackrel{nm}{}}\mathrm{𝟏}\right)\right].`$   
In the very last step above we performed the sums over $`k`$ and $`l`$.

The expression for the potential energy, Eq. (24), has a remarkable symmetry, inherited from Eq. (3). By making the substitutions

$`r`$ $``$ $`nr`$   
$`i`$ $``$ $`j`$   
$`j`$ $``$ $`i`$   
$`𝐱`$ $``$ $`𝐱+𝐲`$   
$`𝐲`$ $``$ $`𝐲`$ (26)   
and noting that $`𝒦_n(𝐲)=(1)^n𝒦_n(𝐲)`$, we can see that the $`r`$th term of Eq. (24) is equal to the $`(nr)`$th term. It also follows that the kernel $`𝒦_n`$ can be chosen to be completely symmetric under interchange of any two of its $`n`$ indices. If we notice that the combinatorial factor

$$\frac{n!}{2^{nm}(2mn)!}$$
(27)   
is precisely equal to the number of distinct ways to assign $`n`$ indices to the tensor $`(^{2mn}𝐲)(^{nm}\mathrm{𝟏})`$, and recalling that only the symmetric part of $`𝒦_n`$ matters, we can rewrite the kernel, Eq. (25), in the remarkably compact form of Eq. (9). The simplicity of this result suggests that there may be an easier way to derive it.

If we now introduce the completely symmetric $`r`$th-rank outgoing tensor fluxes, as defined in Eq. (7), then Eq. (24) may be written as

$`\mathrm{\Delta }V_c=`$ (28)   
$`{\displaystyle \frac{1}{2}}{\displaystyle \underset{𝐱,𝐲}{}}{\displaystyle \underset{n=1}{\overset{\mathrm{}}{}}}{\displaystyle \frac{ϵ^n}{n!}}{\displaystyle \underset{r=0}{\overset{n}{}}}(1)^r\left({\displaystyle \genfrac{}{}{0pt}{}{n}{r}}\right)𝒥_r^{}(𝐱){\displaystyle \stackrel{r}{}}𝒦_n(𝐲){\displaystyle \stackrel{nr}{}}𝒥_{nr}^{}(𝐱+𝐲).`$   
Because of the symmetry, Eq. (26), we can remove the factor of $`1/2`$ and sum $`r`$ from $`n/2`$ to $`n`$, instead of from $`0`$ to $`n`$. The exception to this arises when $`n`$ is even and $`r=n/2`$; in that case, the factor of $`1/2`$ must be retained. We accommodate this case by dividing by $`1+\delta _{n,2r}`$, where $`\delta `$ is the Kronecker delta. If we then define the $`r`$th-rank tensor fields, as in Eq. (8), the expression of Eq. (6) for the potential energy change follows immediately.

As noted in the text, the primes on the fluxes and fields indicate that these are evaluated using the outgoing (post-collision) states. Thus, at each order $`n`$, Eq. (6) expresses the change in potential energy due to the propagation step as the sum of $`n`$ terms, the $`r`$th of which is an $`r`$-fold inner product of an $`r`$th-rank outgoing tensor flux, Eq. (7), with an $`r`$th-rank tensor field, Eq. (8).

Note that we can rearrange the order of summation of $`n`$ and $`r`$ to write Eq. (6) in the following alternative format

$$\mathrm{\Delta }V_c=\underset{𝐱}{}\underset{r=1}{\overset{\mathrm{}}{}}ϵ^r𝒥_r^{}(𝐱)\stackrel{r}{}\left(\underset{n=0}{\overset{r}{}}ϵ^n_{n+r,r}^{}(𝐱)\right).$$
(29)   
This makes it clear that at each order in $`r`$ a new $`r`$th-rank tensor flux must be introduced, and the corresponding $`r`$th-rank field is the sum of $`r`$ terms. The disadvantage of writing the result in this way is that each term of the outer sum over $`r`$ contains terms of differing order in $`ϵ`$.

Finally, we note in passing that tensor fluxes and fields have been considered in the context of a LGA model of microemulsions . While that model did indeed derive its update rule from considerations of single-particle interactions, it did not employ the method used here. In particular, where $`𝐜_j𝐜_i`$ appears in Eq. (3), that study took only $`𝐜_i`$, and the expansion was not carried out to all orders. On the other hand, that work also included vector-valued charge attributes, in order to properly model the orientation of the surfactant molecules. It is likely that a more exact formulation, analogous to our Eq. (3), exists for this microemulsion model, and may well connect this work with the static lattice models of microemulsions due to Matsen and Sullivan .