# Three dimensional hydrodynamic lattice-gas simulations of domain growth and self-assembly in binary immiscible and ternary amphiphilic fluids.

###### Abstract

We simulate the dynamics of phase assembly in binary immiscible fluids and ternary microemulsions using a three-dimensional hydrodynamic lattice gas approach. For critical spinodal decomposition we perform the scaling analysis in reduced variables introduced by [1, 2]. We find a late-stage scaling exponent consistent with the inertial regime. However, as observed with the previous lattice-gas model of Appert et al. [3] our data does not fall in the same range of reduced length and time as that of Kendon et al. [1, 2]. For off-critical binary spinodal decomposition we observe a reduction of the effective exponent with the volume fraction of the minority phase. However, the Lifshitz-Slyzov-Wagner droplet coalescence exponent is not observed. Adding a sufficient number of surfactant particles to a critical quench of binary immiscible fluids produces a ternary bicontinuous microemulsion. We observe a change in scaling behaviour from algebraic to logarithmic growth for amphiphilic fluids in which the domain growth is not arrested. For formation of a microemulsion where the domain growth is halted we find a stretched exponential growth law provides the best fit to the data.

## 1 Introduction

The study of phase ordering kinetics has become a testbed for new complex fluid simulation methods. Despite intensive analysis by many methods, it remains a field in which many interesting and fundamental questions go unanswered. Constructing a model which correctly includes hydrodynamics but which is computationally simple enough to reach late times is a major theoretical challenge. New mesoscale models such as lattice gas (LG) [4, 3], lattice Boltzmann (LB) [5, 2], dissipative particle dynamics (DPD) [1], and Boltzmann-Vlasov [6] treatments have been crucial in increasing our understanding of these systems over the past decade. Such methods require significantly smaller computational resources than earlier molecular dynamics (MD) and Cahn-Hilliard approaches, and can therefore access more easily the late time regime in which hydrodynamics plays an important role.

A much studied system which exhibits hydrodynamic influences on phase segregation is a mixture of immiscible fluids quenched into the two phase region of its phase diagram (a so-called critical quench). Such a system undergoes spinodal decomposition, where initially bicontinuous domain structure coarsens by surface-tension driven flows. Spinodal decomposition elicits much interest because of its fundamental and industrial importance. For example, the mechanical properties of alloys depend on the dynamics of the phase separation process. Despite extensive theoretical [7, 8, 9], numerical [10, 11, 12, 13, 14, 15, 16, 17, 4, 3, 1, 5, 2, 6] and experimental investigation [19, 20] some doubt remains about the true asymptotic late time growth behaviour of these systems.

The dynamics of phase segregation in binary immiscible fluids are also dependent on the composition of the mixture. Mixtures which do not have a ratio of the species (so-called off-critical quenches) are much less studied than their critical counterparts. As the quantity of the minority phase decreases, the domain structure ceases to be bicontinuous, so that nucleation and coalescence of droplets dominate the coarsening mechanism.

The addition of surfactant to a system of binary immiscible fluids radically alters the equilibrium properties and growth dynamics of such mixtures [21]. In particular, the adsorption of surfactant (amphiphilic) molecules at oil-water interfaces leads to a dramatic reduction in the interfacial tension. This property is the origin of much industrial interest in surfactant systems. Only recently have computational techniques and sufficiently powerful parallel platforms become available which permit numerical simulation of the hydrodynamic behaviour of amphiphilic fluids in three dimensions. Our hydrodynamic lattice gas model has recently been implemented in three dimensions, and its ability to reproduce many important phenomenological features of amphiphilic systems confirmed [22]. In the present paper we demonstrate this model’s ability to capture quantitative features of the dynamics of binary immiscible and ternary amphiphilic fluids.

Continuum approaches regard spinodal decomposition as a solution of the equations of motion for two phase flow. These are given by the Navier-Stokes equation within each phase, and by boundary conditions at the interfaces. The incompressible Navier-Stokes equations are, for each fluid phase:

(1) |

(2) |

where is the fluid velocity, is the shear viscosity and is the scalar pressure. The first boundary condition is that the velocity of the two phases must match that of the interface:

(3) |

where is the interface normal and , and are respectively the velocities of phase and at the interface, and the velocity of the interface itself. The second boundary condition requires that the stress difference at the interface is balanced by the interfacial tension:

(4) |

where and denote the stress tensor in phases and , and are the two principal radii of the interfaces and is the interfacial tension.

Clearly, a numerical integration scheme for the above equations would be of enormous computational complexity. In order to obtain a numerical model which captures the essential physics in a computationally tractable way, an alternative approach is required. Remaining with a phenomenological description of the fluid system, it is convenient to define a scalar order parameter :

(5) |

where and are the densities of oil and water at position . One may then write a Cahn-Hilliard equation for the evolution of the order parameter:

(6) |

and,

(7) |

and

(8) |

This equation is coupled to the Navier-Stokes equations through a forcing term proportional to the gradient of the chemical potential:

(9) |

where and are gaussian noise fields satisfying an appropriate fluctuation-dissipation theorem [14].

The equations (5 - 9) define a model for immiscible fluid dynamics commonly referred to as Model H. This model enables one to obtain the scaling regimes for critical spinodal decomposition by dimensional analysis. The viscous regime is obtained when one may neglect the inertial terms on the left hand side of equation (9); then, balancing the forcing terms by the viscous terms one obtains:

(10) |

The inertial regime is obtained by balancing the forcing terms by the inertial terms in eqn (1):

(11) |

Equating these two regimes implies that the crossover from viscous to inertial scaling occurs at . However, no three-dimensional simulation method so far developed can reach a sufficient range of length and timescales to observe both viscous and inertial regimes in a single simulation. The work of Jury et al. [1] and Kendon et al. [2] overcomes this difficulty by introducing scaling variables and which enable data from separate simulations to be combined.

If the only physical effects involved in critical spinodal decomposition are capillary forces, viscous dissipation and fluid inertia, then the parameters governing domain growth are the surface tension , fluid mass density , and viscosity . As emphasised by Jury et al. [1], only one length, , and one time can be constructed from these parameters. Data sets (L,T) from any model of spinodal decomposition can be expressed in units of reduced length:

(12) |

and time,

(13) |

If all other physics is excluded from late-stage growth, then the dynamical scaling hypothesis [8] states that:

(14) |

where is the domain size expressed in reduced length units and is a non-universal constant which allows for a period of growth in which molecular diffusion leads to the formation of sharp interfaces. The function should then approach a universal form, identical for all incompressible binary fluid mixtures following a deep quench.

The Model H equations have been integrated numerically in three dimensions both with (a quench to finite temperature), and with (a quench to zero temperature) by a number of integration schemes [12, 15]. The free energy defined in eqn (8) is also the basis for the lattice-Boltzmann method of Kendon et al. [2]. Although the lattice-Boltzmann method has a physical motivation which replaces velocity and density fields by single particle distribution functions, it is essentially equivalent to direct numerical integration of Model H such as that performed in [12, 15] (i.e. it should be regarded as a finite-difference scheme for solving these equations in the absence of noise). The work of [2, 12, 15] is therefore a confirmation of the scaling laws derived above. The contribution of Kendon et al. [2], in whose work both growth regimes were clearly seen for the first time, shows that the simple arguments given above are incorrect in the crossover regime, and that the crossover in Model H extends over in reduced time units. It should be noted that this crossover was previously believed from simple scaling arguments to occur at .

Recently there has been further theoretical work by Grant and Elder: those authors suggest [23] that growth in the inertial regime implies a Reynolds number that increases without limit. This would eventually lead to turbulent remixing of the fluids; the requirement that the Reynolds number be self-limiting (i.e. that the critical Reynolds number is not reached and therefore that turbulent remixing does not occur) in the asymptotic regime implies a growth exponent of . There is at present no numerical evidence for the , and recent theoretical challenges to this idea have been made. Novik and Coveney point out in [24] that the relative strength of the interface (characterised by the Weber number) must also be taken into account. A small Weber number could delay the onset of turbulent remixing indefinitely.

The hydrodynamic lattice gas used in the present paper does not rest on the macroscopic free energy functional proposed in eqn (8). The model, which is particulate in nature, is described in detail below, and has a theoretical justification from a ‘bottom up’ perspective. In a bottom-up view the lattice-gas technique may be regarded as a simplification of the molecular dynamics of a binary fluid, abstracting the key microscopic properties in a fictitious microworld. Recent work has derived a microscopic basis for the Rothman-Keller model of binary immiscible fluids [4, 25]. However, no systematic method exists for coarse-graining a real molecular dynamics description of a system to the lattice-gas model used here (although such a systematic method does now exist for the DPD algorithm [26]).

In our model equations (1-4) are emergent macroscopic properties. The single phase FHP lattice gas is known to satisfy eqn. (1), and far from interfaces this behaviour is reproduced in our model [27, 30]. In Section 3.2 we demonstrate that our model has realistic surface tension behaviour at interfaces.

The purpose of the present paper is to investigate the dynamics of domain growth of both binary immiscible and ternary amphiphilic phases in our model. Section 2 contains a brief description of our model, while Section 3 contains a description of the quantitative measurements of surface tension and viscosity. In Sections 4, 5 and 6 we present our results for self-assembly in critical and off-critical binary immiscible and ternary amphiphilic fluids respectively. We close the paper with discussion and conclusions in Section 7.

## 2 The lattice gas model

Our lattice-gas model is based on a microscopic, bottom-up approach, where dipolar amphipile particles are included alongside the immiscible oil and water species. Lattice-gas particles can have velocities , where , and is the number of velocities per site. We shall measure discrete time in units of one lattice timestep, so that a particle emerging from a collision at site and time with velocity will advect to site where it may undergo the next collision. We let denote the presence () or absence () of a particle of species (, , denoting red (oil), blue (water) and green (amphiphile) species respectively) with velocity , at lattice site and time step . The collection of all for will be called the population state of the site; it is denoted by

(15) |

where we have introduced the notation for the (finite) set of all distinct population states. The amphiphile particles also have an orientation denoted by . This orientation vector, which has fixed magnitude , specifies the orientation of the amphiphile particle at site and time step with velocity . The collection of the vectors at a given site and time step is called the orientation state. We also introduce the colour charge associated with a given site,

(16) |

as well as the total colour charge at a site,

(17) |

The state of the model at site and time is completely specified by the population state and orientation state of all the sites. The time evolution of the system is an alternation between an advection or propagation step and a collision step. In the first of these, the particles move in the direction of their velocity vectors to new lattice sites. This is described mathematically by the replacements

(18) |

(19) |

for all , and . That is, particles with velocity simply move from point to point in one time step. In the collision step, the newly arrived particles interact, resulting in new momenta and surfactant orientations. The collisional change in the state at a lattice site is required to conserve the mass of each species present

(20) |

as well as the -dimensional momentum vector

(21) |

(where we have assumed for simplicity that the particles all carry unit mass). Thus, the set of population states at each site is partitioned into equivalence classes of population states having the same values of these conserved quantities.

We assume that the characteristic time for collisional and orientational relaxation is sufficiently fast in comparison to that of the propagation that we can model this probability density as the Gibbsian equilibrium corresponding to a local, sitewise Hamiltonian function; that is

(22) |

where is an inverse temperature, is the energy associated with collision outcome , and is the equivalence-class partition function. The sitewise Hamiltonian function for our model has been previously derived and described in detail for the two-dimensional version of the model [31], and we use the same one here. It is

(23) |

where we have introduced the colour flux vector of an outgoing state,

(24) |

the total director of a site,

(25) |

the dipolar flux tensor of an outgoing state,

(26) |

the colour field vector,

(27) |

the dipolar field vector,

(28) |

the colour field gradient tensor,

(29) |

the dipolar field gradient tensor,

(30) |

defined in terms of the scalar director field,

(31) |

and the kinetic energy of the particles at a site,

(32) |

where is the average velocity of all particles at a site, the mass of the particles is taken as unity, and , , , and are coupling constants. To maintain consistency with previous work we use the coupling constants as previously defined in [22]. The values of these constants are:

, , ,

These values were chosen in order to maximise the desired behaviour of sending surfactant to oil-water interfaces while retaining the necessary micellar binary water-surfactant phases. It should be noted that the inverse temperature-like parameter (whose numerical value is varied in this paper) is not related in the conventional way to the kinetic energy. For a discussion of the introduction of this parameter into lattice gases we refer the reader to the original work by Chen, Chen, Doolen and Lee [32], and Chan and Liang [33]. Eqs. (23)-(30) were derived by assuming that there is an interaction potential between colour charges, and that the surfactant particles are like “colour dipoles” in this context [31]. The term parameterised by models the interaction of colour charges with surrounding colour charges as in the original Rothman-Keller model [4]; that parameterised by describes the interaction of colour charges with surrounding colour dipoles; that parameterised by accounts for the interaction of colour dipoles with surrounding colour charges (alignment of surfactant molecules across oil-water interfaces); and finally that parameterised by describes the interaction of colour dipoles with surrounding colour dipoles (corresponding to interfacial bending energy or “stiffness”). This model has been extensively studied in two dimensions [31, 34, 35, 36, 37], and the three-dimensional implementation employed in the present paper is described in more detail by Boghosian, Coveney and Love [22].

## 3 Viscosity and surface tension

In this section we present the results of simulations designed to measure the values of the macroscopic physical parameters which control domain growth in fluid systems in our model, according to the top-down theories described in Section 1. As emphasised above, these parameters are the surface tension , the viscosity and the density .

### 3.1 Viscosity measurements

We measured the viscosity by analysing the decay of shear waves. We performed simulations to observe the decay of shear waves with an initial velocity profile of the form:

(33) |

where is the lattice size in the direction, is the amplitude of the shear wave at time , and is the unit vector in the direction. Solving the Navier-Stokes equations for the time evolution of gives:

(34) |

Where is the initial amplitude of the shear wave and is the kinematic viscosity. We therefore initialise the system with a velocity profile given by (33) and observe the decay of the shear wave by calculating:

(35) |

where indicates summation over all lattice sites in the plane. This gives the mean component of velocity, but includes any net velocity the system may possess. We subtract this to obtain the velocity due only to the shear wave,

(36) |

and calculate as the r.m.s. value of this quantity,

(37) |

We performed simulations at different amplitudes of initial velocity profile and obtained kinematic viscosity in lattice units.

### 3.2 Surface tension analysis

In the present section we analyse the surface tension behaviour in our model for both binary immiscible and ternary amphiphilic systems. The central feature of ternary amphiphilic fluids is the lowering of the interfacial tension between oil and water by the adsorption of surfactant at the interface. The existence of a spinodal point in a binary immiscible fluid is an important feature of the fluid’s thermal behaviour. At temperatures above the spinodal point the fluid will not demix into single phase domains. As one increases the temperature towards the spinodal point the surface tension should be reduced to zero.

We used a direct dynamical method for calculating the surface tension across a flat interface. In the vicinity of an interface the pressure is locally anisotropic, as the pressure in the direction parallel to the interface is reduced by the tension on the interface itself. For a system with flat interface perpendicular to the axis the surface tension is given by the line integral over of the component of pressure normal to the interface minus the component tangential to the interface [38]:

(38) |

Where:

(39) |

(40) |

We begin by showing that our model is capable of reproducing physically realistic interfacial tensions in systems of binary immiscible fluids. We performed simulations for values of , using systems of size . The surface tension was measured every time step for time steps . The simulation with was above the spinodal point. The spinodal point itself was located by a computational steering technique, which was found to be computationally more efficient than traditional taskfarm methods. The value of was modified during a simulation and the phase separation behaviour visualised in real time. The value of at the spinodal point was found to be . This ‘compusteering’ technique represents a powerful and economic simulation technique, and will be the subject of a future paper. Equilibration effects were found to be significant only close to the spinodal point, and equilibration times for , , were taken as , and time steps respectively. All other simulations were allowed to equilibrate for time steps and the surface tension was then time averaged for the remainder of the simulation. The surface tension as a function of is shown in fig. (1).

We next analyse the behaviour of the surface tension as a function of surfactant concentration adsorbed onto an oil-water interface. We use two types of system, both initialised with a flat oil-water interface perpendicular to the direction. The first system is initialised with a fixed amount of surfactant at the interface. The second is initialised with the surfactant dispersed throughout the system. We compute the surface tension as above from eqn. (38).

The equilibrium state for such a system in our model is quite complex. The simplistic macroscopic picture of interfaces in ternary systems is fairly static, with a monolayer of surfactant coating the interface and lowering the surface tension. However, in reality as well as in our model the surfactant may also exist in bulk solution far from the interface, either as monomers or micelles or both. Equilibration of the system is achieved when the rates for adsorption and desorption of monomers and micelles are equal. These times are typically long, varying from to time steps in our model.

For the first type of system we simulate, where all surfactant particles initially reside at the interface, the surface tension prior to equilibration calculated by (38) is an increasing function of time. The surfactant density at the interface decreases with time as monomers and micelles desorb into the bulk. For the second type of system, where all surfactant particles reside in the bulk, the surfactant density at the interface increases as the surfactant adsorbs to the interface. The surface tension calculated by (38) is a decreasing function of time.

We calculate the surface tension as a function of surfactant density in two ways. For low surfactant concentrations where we can reach timescales on which the system is equilibrated we use the surface tension and interfacial density for the equilibrated system. For higher surfactant densities where the equilibration times become prohibitive we plot the surface tension at time against the density at time . The results are shown in fig. (2).

## 4 Critical Spinodal Decomposition

### 4.1 Defining the characteristic size

To analyse the domain growth quantitatively in the following simulations we obtain the first zero crossing of the coordinate-space pair correlation function. This defines the characteristic domain size . The pair correlation function is defined by:

(41) |

where is the colour charge at site , and the integral is taken over the whole system. We compute for each of our simulations, and perform an average over an ensemble of initial conditions. Taking the spherical average of yields , the first zero of which gives the characteristic domain size. We obtain the first zero by performing a linear interpolation between the last point greater than zero and the first point less than zero.

### 4.2 Scaling analysis

The first test of scaling is applied to the correlation functions of Section 4.1. If the domain structures are self-similar at different times during the coarsening, the scaling function defined as

(42) |

should be independent of . Alternatively the scaling criteria may be applied to the Fourier transform of , the structure factor:

(43) |

Spherically averaging this quantity yields S(k,t), which has a similar scaling form:

(44) |

where is the fourier transform of . The functions and for simulations with are plotted as a function of and respectively in figs. 3 and 4. These datasets scale quite well, whereas datasets from simulations with do not. The specific form of g(q) has been the subject of much theoretical attention. The three features which receive the most attention are the small, intermediate and large regimes. We show in figure 4 the behaviour of our model in these three regimes. In the large region we see a clear Porod tail with . This region ends when probes the scale of the interface width. Furukawa speculated that there would be a regime : we see a behaviour closer to , in agreement with Jury et. al. [1]. In the small-q limit we see , with . In [28] the constraint was derived for a Cahn-Hilliard theory without hydrodynamics. However, this derivation assumed a dynamical scaling exponent of and did not include fluctuations. Furukawa speculated that fluctuations would cause the crossover at small .

The above analysis of the correlation functions and structure factors is independent of the form of . We now introduce the reduced length and time variables and as defined in equations (12) and (13). To exclude finite size effects we use data for which the characteristic domain size is less than of the system size. Data from our simulations satisfying this constraint spans a range of and . A study using the DPD algorithm reached a range of and , and a lattice-Boltzmann algorithm has been used to reach a range of and [1, 2]. The analysis method used to plot the data displayed in figures 5 and 6 from simulations with different values of and is identical to that of Kendon et al. [2].

We find that data with show growth with an inertial exponent, and collapses well onto a single scaling function. Simulations with show slow growth where . Visualisation of the order parameter for this data shows that sharp interfaces do not form during the simulation, and so data for these low values of surface tension do not satisfy the criteria for the postulated universal scaling regime. This is confirmed by our scaling analysis of the correlation functions above, where data for did not collapse onto the same curve shown in figure 3. The time evolution of the domain size in reduced units for is shown in figure 5.

This inertial exponent is not consistent with the previous results of Kendon et al. and Jury et al. [1, 2], although it is similar to the results of Appert et al. [3]. In [2], the scaling of Appert et al. was ascribed to ‘excessive diffusion’ in the lattice gas algorithm. The particulate noise characteristic of the lattice gas is not present in their lattice-Boltzmann algorithm, and so it is perfectly plausible that these fluctuations continue to inhibit phase separation, reducing the growth exponent in the viscous regime to some effective exponent close to the usually associated with the inertial regime. An identical effect is well known in two dimensions, where models with fluctuations yield an exponent of in the early time regime, and models without yield an exponent of [40, 34].

This issue deserves some closer examination, however. If we interpret our structure factor and correlation function data as exhibiting goood scaling collapse (i.e. for all times considered the morphology is characterised by a single lengthscale), then the crossover at small q is consistent with speculations by Furukawa that fluctuations would cause the behaviour. As noted above these fluctuations may act to inhibit the phase separation and give rise to a lower effective exponent.

Our immiscible lattice-gas model reduces to a simulation of a convection-diffusion equation for . The diffusion constant in this equation is a function only of the collision rules. It may be possible by a judicious choice of collision rules or addition of rest particles to vary this diffusion constant independent of surface tension to fully investigate its effect.

If, however, we interpret the small behaviour as being indicative of poor scaling collapse, the exponent seen may be interpreted as a crossover between an early time and viscous regime. This issue could only be resolved by larger and longer simulations, and present hardware limitations mean that this must remain a matter for future investigation.

There is at present considerable uncertainty about both the universality of the asymptotic scaling and what the true asymptotic regime is. Theoretical work concerning the dispersion relation on fluid interfaces casts doubt on the universality of hydrodynamic phase separation in three dimensions [29]. Recent work by Kendon [39] proposes a new set of macroscopically deduced scaling laws, based on the energy balance in the system. This analysis suggests that there may be more than one length scale of importance in spinodal decomposition. Another derivation of the growth laws proposes that the is transient and true asymptotic scaling is [18]. However, it would be difficult to resolve such a subtle difference in exponent with any current numerical method. In addition, the derivation of [18] assumes a droplet morphology, even in the symmetric case. No such morphology is seen in any simulation of the symmetric case of which we are aware. The existence of multiple length scales and breakdown of scaling in two dimensional spinodal decomposition is well established [40, 41], and the existence of a single underlying scaling function in three dimensions remains an open question.

## 5 Kinetics of phase separation for off-critical binary fluids

In this section we analyse the dynamics of domain growth in systems where the composition is asymmetrical. In such systems the domain structure is not bicontinuous. The minority oil (or water) phase exists as droplets, and the domains coarsen by diffusion of oil (or water) from the bulk onto the droplet, and by droplet coalescence. The composition for a system where oil is the minority phase is defined as:

(45) |

where and are the reduced densities of oil and water respectively. We performed simulations for values of , and . These simulations were performed on systems for time steps. The correlation function was averaged over five independent simulations for , and three independent simulations for . The results are shown in figures 8 and 9.

The data in our simulations gives effective exponents of and for compositions and respectively. We do not observe the exponent expected from simple theories of droplet coalescence and nucleation. However, as fig. 7 shows, the morphology for both these volume fractions is far from being a collection of isolated drops. It is likely that the exponent we measure is therefore intermediate between the critical case and droplet coalescence and nucleation. The existence of such intermediate exponents is well established. Previous work by Appert et al. [3] on off-critical decomposition for also saw growth more rapid than the expected for simple nucleation and coalescence. The work of [42] for the two-dimensional implementation of our model similarly saw a continuous variation of exponent with composition, as did [24] for equal viscosity DPD fluids.

## 6 Self assembly in ternary amphiphilic fluids

In this section we turn to the analysis of ternary amphiphilic fluids. We concentrate on systems with equal amounts of oil and water, varying the amount of surfactant for each simulation. The presence of surfactant dramatically alters the interfacial energetics and structure, and consequently the dynamics of domain growth. Specifically, the adsorption of surfactant at oil-water interfaces and its concomitant lowering of the surface tension weaken the forces which drive binary immiscible phase separation. For sufficiently large surfactant concentrations the final equilibrium state is a bicontinuous or sponge microemulsion phase, where the domain growth is arrested at some final characteristic domain size . Such a system is depicted in fig. (10). All three components are bicontinuously connected, with the surfactant particles lining the interface in a monolayer between the percolating oil and water regions.

We work with and use the values for the coupling coefficients in eqn (23) as defined earlier (in Sec. 2). The results that follow for the ternary emulsion system have been obtained from a system with periodic boundary conditions in all directions, the system having been intialised from a random quench with the particles placed randomly on the lattice. The total reduced density of the system was kept constant at . The characteristic domain size was measured from the spherically averaged pair correlation function as described in Section 4.1. The growth of sponge as opposed to droplet phases in these systems means that for the lowest values of the surfactant concentration the growth behaviour should represent a perturbation of the critical quenches investigated in Section 4. The results for reduced surfactant concentrations , , and support this conclusion. In all cases the late-time behaviour is consistent with an algebraic growth law of the form . Prior to this there is an an early-time regime in which the growth is consistent with a diffusive algebraic exponent of . Visualisation of the surfactant densities during this period shows that surfactant adsorbs at the oil-water interfaces.

Once the system reaches the biconinuous oil and water state the inertial hydrodynamic regime begins. The length of the diffusive period increases with increasing surfactant concentration ( time steps for reduced surfactant concentration , time steps for , time steps for and time steps for ). As one increases the surfactant concentration the fluctuations in the domain size at late times increases. This is consistent with the known dynamical nature of the adsorption and desorption of surfactant to and from interfaces in ternary systems as one approaches emulsification. With a reduced surfactant density of the time evolution of the domain size ceases to follow the algebraic growth law. We see a slower-than-algebraic growth law, as previously observed in the two dimensional implementation of our model where a growth function was proposed, based on a comparison with slow growth in systems with quenched impurities [34, 43]. Consequently we look at a plot of against domain size in order to determine whether we have logarithmically slow growth in this regime. The characteristic domain size for surfactant concentrations are plotted against on logarithmic scales in fig. (12). One can clearly see a transition through a regime where logarithmically slow growth dominates throughout the timescale of the simulation. Such a growth law is inconsistent with arrest of the domain growth at late times for reduced surfactant density below . As we increase the reduced surfactant density beyond we see complete arrest of the domain growth at late times. We performed a fit of a stretched exponential function to datasets from simulations with . This function is defined identically with [34]:

(46) |

These fits are shown in fig. (13). To quantify the effect that the surfactant has on the domain size at late times in these simulations, following Emerton et al. [34] and Gyere et al. [43], we define as the domain size at time step , where data which is unaffected by finite size effects is available for all surfactant concentrations. We plot against the inverse of the reduced surfactant density at the interface (having subtracted the micellar and monomeric concentration) in fig. 14 (for systems with quenched impurities the amplitude of the disorder is analogous to the surfactant concentration). We find a linear dependence of on , consistent with the results found in [34] and [43].

## 7 Conclusions

We have studied the dynamics of domain growth in binary immiscible fluids, in both critical and off-critical cases, and in ternary (oil-water-surfactant) emulsions and microemulsions, using a three-dimensional hydrodynamic lattice gas model. In the critical binary case we performed a scaling analysis similar to that of Jury et al. [1] and Kendon et al. [2] and find late-time growth behaviour which is consistent with an inertial hydrodynamic exponent. However the position of the late-time domain growth on a plot of reduced time variables is not consistent with that of [1, 2], in the same way that the results of Appert et al. [3] showed an inertial exponent in a crossover region of the reduced scaling plot. In the off-critical case we find effective exponents of and for compositions and respectively. In the microemulsion case we observe a slowed algebraic growth with exponent , followed by a regime in which we see evidence for logarithmically slow growth. For concentrations of surfactant high enough to completely arrest domain growth we observe good agreement with a stretched exponential fit to our data. Overall, all our results are fully consistent with behaviour we observed previously in our two-dimensional lattice gas model.

## Acknowledgements

We are indebted to numerous people and organisations for their support of and contributions to this work. They include Jean-Bernard Maillet, Keir Novik, Nelido Gonzalez, Jeremy Martin, Mike Cates and Julia Yeomans. Simulations were performed on the Oscar Origin 2000 service at Oxford Supercomputing Centre and the Cray T3E at the Manchester CSAR service. Resources for the T3E were allocated under EPSRC grant number GR/M56234. PJL would like to thank EPSRC and Schlumberger Cambridge Research for funding his CASE studentship award. The collaboration between PVC and BMB was supported by NATO grant number CRG950356. BMB was supported in part by the United States Air Force Office of Scientific Research under grant number F49620-99-1-0070.

## References

- [1] S. Jury, P. Bladon, S. Krishna, and M. E. Cates, Phys. Rev. E 59, R2535 (1999).
- [2] P. Bladon, V. Kendon, J.-C. Desplat and M. E. Cates, Phys. Rev. Lett 83, 579 (1999).
- [3] C. Appert, J. Olson, D. H. Olson, and S. Zaleski., J. Stat. Phys. 81, 181 (1995).
- [4] D. Rothman and J. Keller, J. Stat. Phys. 52, 1119 (1988).
- [5] F. Alexander, S. Chen, and D. Grunau, Phys. Rev. B. 48, 634 (1993).
- [6] S. Bastea and J. L. Lebowitz, Phys. Rev. Lett. 78, 3499 (1997).
- [7] A. Bray, Adv. in Physics 43, 357 (1994).
- [8] E. Siggia, Phys. Rev. A. 20, 595 (1979).
- [9] V. Kumaran, J. Chem. Phys 108, 3038 (1998).
- [10] A. Chakrabarti, R. Toral, and J. Gunton, Phys. Rev. B 39, 4836 (1989).
- [11] A. Shinozaki and Y. Oono, Phys. Rev. Lett. 66, 173 (1991).
- [12] T. Koga and K. Kawasaki, Phys. Rev. A. 44, 817 (1991).
- [13] S. Puri and B. Dunweg, Phys. Rev. A. 45, R6977 (1992).
- [14] O. Valls and J. E. Farrell, Phys. Rev. E. 47, R36 (1993).
- [15] T. Lookman, Y. Wu, F. Alexander, and S. Chen, Phys. Rev. E 53, 5513 (1996).
- [16] W. Ma, A. Maritan, J. Banavar, and J. Koplik, Phys. Rev. A. 45, R5347 (1992).
- [17] M. Laradji, S. Toexvaerd, and O. Mouritsen, Phys. Rev. Lett. 77, 2253 (1996).
- [18] F. J. Solis and M. O. de la Cruz, Phys. Rev. Lett. 84, 3350 (2000).
- [19] T. Hashimoto, H. Jinnai, H. Hasegawa, and C. C. Han., Physica A. 204, 261 (1994).
- [20] P. Guenoun, R. Gastaud, F. Perrot, and D. Beysens, Phys. Rev. A. 36, 4876 (1987).
- [21] G. Gompper and M. Schick, Phase Transitions and Critical Phenomena 16, 1 (1994).
- [22] B. M. Boghosian, P. V. Coveney, and P. J. Love, Proc. Roy. Soc. London A. 456 1431 (2000).
- [23] M. Grant and K. Elder, Phys. Rev. Lett. 82, 14 (1999).
- [24] K. E. Novik and P. V. Coveney, Phys. Rev. E. 61, 435 (2000).
- [25] B. Boghosian and P. V. Coveney, Computer Physics Communications 129, 46 (2000).
- [26] E. Flekkøy and P. V. Coveney, Phys. Rev. Lett. 83, 1775 (1999).
- [27] S. Wolfram, J. Stat. Phys. 45, 471 (1986).
- [28] C. Yeung, Phys. Rev. Lett. 61, 1135 (1988).
- [29] C. Yeung, Phys. Rev. E. 48, 1984 (1993).
- [30] U. Frisch, D. d’Humières, B. Hasslacher, P. Lallemand, Y. Pomeau and J. P. Rivet, Complex Systems 1, 1 (1987).
- [31] B. Boghosian, P. V. Coveney, and A. Emerton, Proc. Roy. Soc. A 452, 1221 (1996).
- [32] H. Chen, S. Chen, G. Doolen, and Y. Lee, Phys. Rev. A. 40, 2850 (1989).
- [33] C. Chan and N. Liang, Europhys. Lett. 13, 495 (1990).
- [34] A. N. Emerton, P. V. Coveney, and B. M. Boghosian, Phys. Rev. E 55, 708 (1997).
- [35] P. V. Coveney, A. N. Emerton and B. M. Boghosian, J. Amer. Chem. Soc. 118, 10719 (1996).
- [36] A. N. Emerton, P. V. Coveney, and B. M. Boghosian, Physica A. 239, 373 (1997).
- [37] A. N. Emerton, P. V. Coveney, and B. M. Boghosian, Phys. Rev. E. 56, 1286 (1997).
- [38] J. Rowlinson and B. Widom, Molecular Theory of Capillarity (Clarendon, Oxford, 1982).
- [39] V. Kendon, Ph.D. thesis, Edinburgh University, (1999).
- [40] A. J. Wagner and J. M. Yeomans, Phys. Rev. Lett. 80, 1429, (1998).
- [41] A. J. Wagner and C. Scott, cond-mat/9911067 (1999).
- [42] F. Weig, P. V. Coveney, and B. M. Boghosian, Phys. Rev. E. 6877 (1997).
- [43] M. Gyure, S. Harrington, R. Strika, and H. E. Stanley, Phys. Rev. E. 52, 4632 (1995).