{'fraction_non_alphanumeric': 0.029503825643403665, 'fraction_numerical': 0.00969935852847979, 'mean_word_length': 4.496707731520815, 'pattern_counts': {'":': 0, '<': 4, '<?xml version=': 0, '>': 3, 'https://': 0, 'lorem ipsum': 0, 'www.': 0, 'xml': 0}, 'pii_count': 0, 'substrings_counts': 0, 'word_list_counts': {'cursed_substrings.json': 2, 'profanity_word_list.json': 0, 'sexual_word_list.json': 0, 'zh_pornsignals.json': 0}}
Sedimentation of a two-dimensional colloidal mixture exhibiting liquid-liquid and gas-liquid phase separation: a dynamical density functional theory study

We present dynamical density functional theory results for the time evolution of the density distribution of a sedimenting model two-dimensional binary mixture of colloids. The interplay between the bulk phase behaviour of the mixture, its interfacial properties at the confining walls, and the gravitational field gives rise to a rich variety of equilibrium and non-equilibrium morphologies. In the fluid state, the system exhibits both liquid-liquid and gas-liquid phase separation. As the system sediments, the phase separation significantly affects the dynamics and we explore situations where the final state is a coexistence of up to three different phases. Solving the dynamical equations in two-dimensions, we find that in certain situations the final density profiles of the two species have a symmetry that is different from that of the external potentials, which is perhaps surprising, given the statistical mechanics origin of the theory. The paper concludes with a discussion on this.


I. INTRODUCTION
Colloid science is a field where physical chemistry and statistical mechanics meet. Colloidal suspensions exhibit a number of interesting phenomena, such as flocculation, adsorption (e.g. of polymers), bridging, osmosis, Tyndall effect, etc. For applications, it is important to understand and control the stability of colloidal suspensions, i.e., to determine the conditions under which the colloids remain dispersed or when the colloids gather to form clusters of a macroscopic size. The interactions between colloids are often dominated by excluded volume effects, and so statistical mechanical theories developed for simple fluids, based on a simple molecular model description, can often be used to understand fundamental properties of colloidal systems. Even crude approximations, where only the repulsive forces between colloids are considered, can give good insight into the phase behaviour of mixtures of colloids and the effects of packing, or 'depletion' interactions, between the different species. 1 Clearly, more complex models are required for charged colloids that may be treated using DLVO (Derjaguin-Landau-Verwey-Overbeek) theory, or for mixtures of colloids and (non-adsorbing) polymers, for which a simple Asakura-Oosawa-Vrij potential provides a surprisingly good description. 1 Colloidal fluids are also interest because of the ability to observe directly the fluid phenomena which occur on length scales accessible to confocal microscopy and which are much harder to observe in molecular fluids. For example, in Refs. 2,3 the detailed dynamics of fluid interfaces and the coalescence of liquid drops could be studied at the particle level.
There is also considerable interest in two-dimensional (2D) model colloidal fluids. These enable the study of the properties of colloidal suspensions which are effectively pinned to a 2D plane; this is the case for par-ticles adhered by capillary forces to a fluid interface 4 (either a liquid-gas or liquid-liquid interface) or those whose motion in one dimension is restricted by an applied external field. For example, 2D models have been studied by means of computer simulations 5,6 and density functional theory 7 to understand the origin of the structures formed by certain colloidal monolayers at the air-water interface 4,8 . In these particular systems one observes complex microphase-separated equilibrium structures of various morphologies that form due to a competition between a short-range attraction and a longerranged repulsion between the colloids. Colloids at the liquid-gas interface experience a capillary attraction between themselves, that can lead to aggregation and phase separation. 4,9 Colloidal monolayers that have sedimented onto the flat bottom of a container can also act as a 2D fluid, and binary colloidal mixtures of such particles can exhibit freezing 10 and phase separation 11 .
A description of the dynamical properties of such systems presents a challenging task. Dynamical properties of fluids are commonly studied by solving the Navier-Stokes equations that treat the fluid as a continuum. However, a more microscopic theory is required, if one wishes to obtain a description on the colloidal lengthscale. Over the last decade or so, a dynamical density functional theory (DDFT), originally formulated by Marconi and Tarazona, 12,13 has been developed to offer such a description for colloidal fluids. This approach stems from the assumption that the colloids can be modelled as an assembly of Brownian particles whose motion is governed by stochastic (Langevin) equations of motion. By making the assumption that two-body correlations in a fluid out of equilibrium can be approximated by those of an equilibrium fluid with the same one-body density profile, a deterministic equation of motion for the onebody density ρ(r, t) can be derived. [12][13][14][15] The resulting equation of motion for the time evolution of ρ(r, t) depends on the functional derivative of the equilibrium free energy functional, F [ρ]. Thus, the theory takes as input the accurate and verified approximations for F [ρ] that have been developed over the last 30 years or so in classical density functional theory (DFT). 16 The advantage of basing the dynamical theory on equilibrium DFT is that the approach guarantees (at the very least) a reliable prediction of the equilibrium fluid density profile in the long time (equilibrium) limit, for external potentials that have no time dependence.
In our previous work in Ref. 17, we studied the interplay between sedimentation under the influence of a lateral driving force (e.g. due to gravity) and the gasliquid phase separation of a 2D one-component colloidal suspension. Note that by colloidal 'gas' we mean a low density suspension and by 'liquid' we mean a high density suspension. This phase separation is driven by the attractive forces between the colloids and is completely analogous to the gas-liquid phase separation of molecular fluids. During this colloidal sedimentation the attractive forces between the colloids can lead to aggregation and the lateral driving also plays an important role in determining the final equilibrium state of the system and the degree to which the colloids are spread over the hard wall at the bottom of the system. The free energy functional used in Ref. 17 is a rather simple approximation, consisting of a local density approximation (LDA) for the hard core repulsive interactions between the colloids and a simple mean-field approximation for the attractive part of the pair interactions. Nevertheless, comparison of the DDFT results with Brownian dynamics simulations revealed surprisingly good agreement. However, this was only achieved when the DDFT was solved as a 2D problem allowing the system to exhibit, in some situations, symmetry breaking in the x-direction, which is the direction parallel to the wall (and perpendicular to the lateral driving force acting in the y-direction). When we also calculate the time evolution of the density profile assuming that it only varies in the y-direction (the direction in which the external potential varies), for cases where symmetry breaking is observed in the 2D DDFT calculation, we find that the agreement between the one-dimensional (1D) DDFT results with simulation is rather poor, with a significant overestimation of the average density in the vicinity of the wall. Much better results are obtained by averaging over the symmetry-broken density distribution obtained from the 2D DDFT -i.e. averaging ρ(x, y) over the x-direction. These 2D density profiles correspond to a series of distinct drops of the colloidal liquid phase. Such states of course only occur due to the fact that the system consists of a fixed number of colloids confined within a finite size system. The reasons for this is discussed in detail in Ref. 17 and we also return to this matter in the last section of this paper.
In the present paper, we generalise our previous work 17 by considering a 2D binary mixture. The good agreement between the simulation and DDFT results for the one-component fluid gives us confidence that the DDFT also provides a reliable description of the mixture. Thus, here we solely present a DDFT treatment of the binary fluid, going to systems with sufficient numbers of particles that simulation is prohibitive. The interplay between the lateral driving (hereafter assumed to be due to a gravitational field and referred to as such), the (bulk) phase separation and "surface" (strictly, line, in 2D) effects, from the interfaces that are produced between the different phases, is shown to produce an extremely rich behaviour. Our goal here is not to construct a global 'phase diagram' of all the possible morphologies that can be observed for this system, which would be an extremely lengthy task, even in the case when only a restricted number of the parameters of the model are varied. Instead, we display some typical results which illustrate the complexity of the confined mixture, giving a detailed discussion of the characteristic behaviour of selected systems. We find that, in some situations, small changes of the system composition can lead to quite dramatic changes in the behaviour of the system. The model differs from the one-component fluid because the binary system has the possibility of exhibiting three phase coexistence. Furthermore, not only is the bulk fluid phase behaviour of the mixture more complex than its one-component counterpart, but the dynamics is also substantially richer.
This paper proceeds as follows: We start in Sec. II by defining our model system, specifying the nature of the interactions between the particles and the external field. The DDFT equations are then briefly described, together with the simple approximation for the free energy functional that we use. In Sec. II C we calculate the equilibrium bulk fluid phase diagram. In Sec. III we present our 2D DDFT results. We first focus on the regime where only liquid-liquid demixing is observed. Following this, we consider the regime where the system exhibits phase separation to form states with three-phase coexistence. We conclude with a summary of our results and a discussion of some of the underlying theoretical issues.

A. Model and equation of motion
We consider a 2D model binary mixture of colloids that interact through the pair potentials where the index i, j = 1, 2 labels the two different species of colloids in the system. The pair potential is composed of two contributions, the purely repulsive hard-disk potentials and we assume the following simple form for the attractive contribution to the pair potentials Here, r = (x a − x b ) 2 + (y a − y b ) 2 is the distance between a pair of particles (labelled with the indices a and b), whose centers are at r a = (x a , y a ) and r b = (x b , y b ). σ ij = (σ ii + σ jj )/2, where σ ii is a hard-core diameter of a species i, and the parameter ε ij > 0 determines a strength of the attractive interactions. In total, there are N = N 1 + N 2 particles in the system, where N 1 and N 2 are the number of particles of species 1 and 2, respectively. We consider the binary fluid when it is confined between two hard walls at y = 0 and y = L and also under the influence of a constant lateral driving force due to gravity, with amplitude |a i |. Thus, the external potential exerted on a particle of species i is: where the hard-wall potential We define a i < 0 so that the colloids are pushed down to the lower wall at y = 0. Note too that in our DDFT calculations we replace the hard wall potentials ϕ hw i by the following slightly softened potentials with γ = 10 and where β = 1/k B T is the inverse temperature; T is the temperature and k B is Boltzmann's constant. Using (6) rather than (5) has almost no noticeable effect on the results, but it does make the numerical computations more stable. We assume that the dynamics of the colloids is governed by the following over-damped stochastic equation of motion:ṙ where X a (t) is a white noise term, Γ i ≡ βD i is the mobility, where D i is the diffusion coefficient for colloids of species i, and with i(a) = 1 for a ≤ N 1 and i(a) = 2 for a > N 1 , is the potential energy of the system. 18 B. DFT and DDFT for the system DFT states that the equilibrium fluid density profiles of a binary mixture are those which minimise the grand potential functional: 16 (9) where F is the Helmholtz free energy functional and µ i is the chemical potential of a species i. The Helmholtz free energy functional can be written as a sum of several contributions: The first term is the ideal-gas contribution, where Λ i is the thermal de Broglie wavelength of species i. The second term F ex is the excess contribution, due to the interactions between the colloids and the final term is the contribution due to the external potentials. The excess free energy can be further split into two distinct contributions: The first from the repulsive part of the interactions between the particles and the second from the attractive part, i.e. from the two different terms in the right hand side of (1), respectively: The second term, stemming from the attractive interactions between the colloids, is approximated in a simple mean-field manner: dr dr ρ i (r)ρ j (r )u att ij (|r − r |) . (12) For the first term in Eq. (11), accounting for the repulsive hard-disk part of the interactions, we adopt a simple local density approximation where f hd is a free energy density for a mixture of harddisks with uniform one-body densities {ρ i }. Here, we use the approximation given by scaled particle theory 19 where the packing fraction η = π 4 i ρ i σ 2 ii ≡ i η i , is the sum of the partial packing fractions, and the total density ρ = ρ 1 + ρ 2 is the sum of the densities of the two species.
Instead of the LDA approximation in Eq. (13), an improvement would be to use the non-local fundamental measure theory for hard disks, developed in Ref. 20. However, for the present system, due to the fact that the confining walls are hard, the LDA is surprisingly good. 17 This is because the density near to the wall is rather low and so the density profiles have very little of the oscillatory structure due to packing effects that fluids normally have close to a wall to which the particles are attracted. The LDA is unable to describe such density oscillations, but since the amplitude of the oscillations is very small in the present system, the LDA actually performs reasonably well. 17 We now consider the statistical mechanics of the system when it is out of equilibrium. The dynamics of the system can be described using DDFT. The quantity of interest here is the non-equilibrium fluid one-body density profiles ρ i (r, t). This is defined as an average over all possible realisations of the white noise in Eq. (7), starting from an ensemble of particle positions consistent with the initial system set-up. In the present case, we consider a system where the particles are initially uniformly distributed between the two walls. On integrating the Smoluchowski equation corresponding to Eq. (7), we obtain the following pair of coupled equations for the time evolution of the one-body densities: [12][13][14][15]18 (15) In deriving Eq. (15) the assumption is made that the instantaneous two-body correlations can be approximated by those corresponding to an equilibrium system with the same one-body density distributions. It is this assumption that leads to the equilibrium Helmholtz free energy functional being input into Eq. (15). Additionally, here we approximate F [{ρ i (r)}] using the functional given in Eqs. (10)- (12).

C. Bulk fluid phase behaviour
In the absence of any external potentials, the fluid densities are uniform and so the free energy density of our model fluid becomes: where V is the area of the 2D system and α ij ≡ u att ij (r)dr = −2πε ij σ 2 ij . For two bulk phases A and B to coexist, the temperature, chemical potentials and pressure in the two phases must be identical, i.e.: In the present system, the chemical potentials are given by: and the pressure is Before proceeding to the description of the phase behaviour of the binary mixture, it is useful to construct the phase diagram for the model involving only one component. The temperature versus density phase diagram for the one-component fluid is displayed in Fig. 1, where both binodal and spinodal are displayed. The binodal is the locus of the coexisting densities and the spinodal is given by the condition ∂ 2 f ∂ρ 2 = 0. It defines the limit at which the uniform density state is linearly unstable. 14 The critical density ρ c σ 2 = 0.274 and critical temperature k B T c /ε = 0.685 follow from the additional condi- 3 ρc,Tc = 0. Thus, when the uniform bulk fluid is cooled below the critical temperature, the system can phase separate, exhibiting gas-liquid phase coexistence. Returning now to the binary fluid, in this case the spinodal limit of linear-stability condition is: This equation defines a temperature surface, as a function of the dimensionless total density ρ * = ρ 1 σ 2 11 +ρ 2 σ 2 and concentration x = ρ 1 σ 2 11 /ρ * , below which a uniform fluid is linearly-unstable and will phase separate, exhibiting either gas-liquid coexistence, which is an extension of the phase transition exhibited by the one-component fluid already discussed above, or the system can also exhibit liquid-liquid phase separation or even three phase gas-liquid-liquid phase coexistence. An example of this spinodal surface is displayed in Fig. 2.
In the remainder of this paper, we consider the case of a symmetric binary mixture, with σ 11 = σ 22 = σ 12 ≡ σ, with ε 11 = ε 22 ≡ ε and with equal mobilities Γ 1 = Γ 2 = Γ. Thus, the like-species interactions for both species are identical. However, the cross-interaction potential between colloids of species 1 with species 2 is not the same, because we choose ε 12 = ε, specifically, we set ε 12 /ε = 0.04. Thus, the present system is rather similar to the 2D binary colloidal fluid studied in Ref. 11. How-ever, in the present work we use the much simpler form for the attractive interactions between the colloids given in Eq. (3). Taking σ and ε as the length and energy units, respectively, we define reduced quantities such as the dimensionless temperature T * = k B T /ε, density ρ * = ρσ 2 and time t * = (Γk B T /σ 2 )t, used throughout this paper. The spinodal surface for the system with ε 12 /ε = 0.04 is displayed in Fig. 2. We do not display the binodal surface obtained from solving Eqs. (17), which lies above the spinodal surface, touching it only along lines of critical points. Slices through such a coexistence surface can be seen for example in Refs. 21-23.

III. DDFT RESULTS
In this section we present typical results for the dynamical properties of the model as determined from DDFT calculations. Eqs. (15) are solved by discretising the 2D system on a Cartesian grid of size 64 σ × 64 σ, with gridspacing ∆x = ∆y = 0.5σ. The two hard walls are located at y = 0 and y = L and separated by a distance L = 45σ. There are periodic boundary conditions in the x-direction, that runs parallel to the walls. We set the initial time t = 0 state to correspond to that where the particles are uniformly mixed and uniformly distributed between the two confining walls. The corresponding density profiles are ρ i (x, y) =ρ i + δ ρ (x, y) for 0 < y < L and ρ i (x, y) = 0 otherwise, whereρ i is the average fluid density of species i in the system and δ ρ is a small random value, which is added as a perturbation to the density profiles in order to allow the system to form states that break the symmetry of the confining potentials; δ ρ (x, y) is a randomly generated number from the uniformly distributed interval (−ρ i /20,ρ i /20). The gravitational force pushes the colloids towards the bottom confining wall at y = 0. This is the same confining potential as was considered in our earlier work in Ref. 17, for the one-component fluid.
For the one-component fluid there are two main factors that determine the final state of the system: 17 i) the strength of the external field pushing the particles towards the bottom of the system and ii) the strength of the interactions between the particles (in particular, whether it is strong enough to lead to phase separation). Since the confining walls are purely repulsive, the colloidal density right near the surface of the walls is low when there is attraction between the particles. This is because any given colloid that is close to the wall can lower its energy by moving away to be closer to the other colloids -i.e. the colloidal fluid does not spread at the wall, unless the gravitational force pushing the particles downwards is strong enough to overcome the cohesive forces between the colloids, or the density of particles in the system is high enough that the weight of the colloids above or the upper wall pushes the lower colloids onto the bottom wall. If this is not the case and when the attraction between the colloids is sufficiently strong, then the colloids gather together to form a series of drops on the bottom wall, with a contact angle that can be close to 180 • .
In the present case of a binary mixture, much more complex behaviour is observed. This is, of course, partly due to a richer bulk fluid phase behaviour, compared to the one-component case -see Sec. II C. In addition, the effect of the external fields, both the gravitational component and the hard-walls, is particularly important and gives rise to a number of different morphologies of the fluid. Furthermore, as the system evolves it can form one morphology that then can become unstable after a certain time period and may then transform to another morphology before the equilibrium is reached. When the gravitational force is relatively weak, at early times, the main effects on the dynamics of the system are the thermodynamic forces driving the demixing of the fluid, coupled with the presence of the confining walls, which together lead to a surface-directed spinodal decomposition. This process is well understood in simple models and at early times in the dynamics, the behaviour of the present system is very similar to that observed in the more simple systems. 24,25 However, the system we consider here also exhibits gas-liquid demixing (in addition to the liquidliquid demixing) and this, together with the gravitational force pushing the particles towards the lower wall, lead to a very rich later stage dynamics.
We start our survey of typical results with the case when both species are subject to an equally strong external field, a 1 = a 2 = 0.01ε, and when the dimensionless temperature T * = 1 and average density in the system ρ * = 0.6, corresponding to a state point where there is no liquid-gas separation but liquid-liquid demixing is possible, depending on the concentration of the mixture. In Fig. 3 we display a time series of plots of the local species 1 concentration, defined as follows: ρ1(r,t)+ρ2(r,t) . The total average concentration in the system x = N1 (N1+N2) is fixed. In Fig. 3 we see in the early stages the strong influence of the walls and the system forms several alternating horizontal layers with the layers closest to the walls rich in the minority species 1. The demixing is due to the fact that the attraction between like species particles is much stronger than that between unlike species. The layering is parallel to the walls because the presence of the walls enhances the growth of demixing-density fluctuations perpendicular to the wall. However, due to the fact that there is a greater amount of species 2 in the system, such a configuration is not energetically favourable because the total amount of interface between the two demixed phases is fairly large in this configuration and so is the interfacial tension contribution to the free energy. The system can lower the free energy by reducing the amount of interface between the two demixed phases. Consequently, the thin layers rich in species 1 become eventually unstable with respect to the density fluctuations parallel to the walls and breakup into a number of droplets. It is interesting to note, however, that this instability hits only the inner layers while the layers at the top and at the bottom nearest to the walls survive long after the first droplets appear. This is because the walls have a stabilising effect on these layers. The distribution of radii of the drops emerging as the layers break-up is quite broad. The larger drops then grow at the expense of the smaller ones until the small ones vanish. Since the gravitational force is relatively weak and equal for both species in this case, the three resulting drops remain in the center of the system, surrounded by species 2 rich fluid. However, as soon as one of the drops (the one on the left in our plot) slightly shifts down, whether due to the external field or due to the randomness in the initial density profiles, a disturbance of both bottom and top layers becomes evident. This is because as the drop moves down, it pushes the fluid beneath which, in turn, pushes on the surface of the bottom layer. At the same time, an inverse process occurs above the drop: here the fluid is dragged in the direction of the moving drop which causes a local thickening of the top layer. Following this, over a rather long time period, both of the remaining stripes become unstable and break up into drops on the walls.
In Fig. 4 we display results for a higher value of the total concentration, x = 0.4. At early times, the system behaves similarly to in the previous case in Fig. 3 (which was for x = 0.3). Again, four layers rich of species 1 are formed, separated by three layers of liquid rich in species 2. However, in this case, the inner layers rich in species 1 are much thicker than previously, because of the increased amount of species 1 particles in the system. We clearly observe the growth of modulations along the borders of the inner stripes. However, there is the competing effect of a transfer of material from the outer stripes (adjacent to the walls) to the inner stripes, which results in the inner ones becoming 'fatter', so that they remain intact over a long time period. Interestingly, the mechanism for the transfer of material from the outer stripes to the inner ones is via diffusion through the other liquid phase, rather than via the coalescence of drops, as was the case in Fig. 3. During this process the top layer continuously shrinks, but the layer at the bottom additionally becomes unstable and forms drops at the wall that then vanish over time in favour of the layer above. The remaining two thick layers seem to be quite stable, although breaks in the layer can and do form.
Next, we increase the amount of species 1, so that the concentration is x = 0.5, i.e. there is exactly the same amount of each species of particles in the system. A time series of local concentration profiles is displayed in Fig. 5. Owing to the fact that this is a symmetric mixture, the effect of the wall in directing the spinodaldecomposition by favouring the minority species is not present in this case. Instead, the system initially phase separates away from the walls in the bulk of the fluid in an isotropic fashion. However, over time, the effect of the walls becomes evident in the orientation of the local concentration profiles and also the interfaces between the two phases sharpen. There is also a striking difference from the previous case: Since in this case the mixture is completely symmetric and therefore neither species can be preferentially adsorbed at the walls, the resulting morphology corresponds to vertical rather than horizontal stripes.
So far, we have considered systems where the temperature and density are both rather high and shown typical microstructure morphologies characteristic for systems exhibiting liquid-liquid demixing. It should be noted that in these cases the effect of gravity is small and the evolution of these systems is very similar to the case with no gravity, a i = 0. For stronger fields, such as a 1 = a 2 = 0.1ε (not shown here), the shift downwards of the center of mass of the entire system becomes substantial, leading to strong gradients of the density along the y-direction. The difference in densities at the top and bottom of the system can be as high as 50%. Nevertheless, apart from that, the sedimentation relaxation remains qualitatively similar to that presented above for the weaker field. We now show results corresponding to the lower temperature T * = 0.5, i.e. below the critical temperature of the gas-liquid phase separation, c.f. Fig.  FIG. 6: A time series of local concentration profiles c(r, t), (upper panels) and the corresponding total number density ρ(r, t) ≡ ρ1(r, t) + ρ2(r, t) (lower panels) for a1 = a2 = 0.01ε,ρ * = 0.3, T * = 0.5, and x = 0.3. The plots correspond to the following times (from left to right): t * = 100, 400, 500, and 600. FIG. 7: A time series of local concentration profiles c(r, t), (upper panels) and the corresponding total number density ρ(r, t) ≡ ρ1(r, t) + ρ2(r, t) (lower panels) for a1 = a2 = 0.1ε,ρ * = 0.3, T * = 0.5, and x = 0.3. The plots correspond to the following times (from left to right): t * = 100, 200, 400, and 600.

2.
We also consider cases where the overall density in the system is lower, so that there is the possibility of observing the gas phase in addition to the two liquid phases. In Fig. 6 we display a time series of local concentration profiles c(r, t) and in addition the corresponding total density profiles ρ(r, t) ≡ ρ 1 (r, t) + ρ 2 (r, t) for a system with average overall density ρ * = 0.3, concentration x = 0.3 and again subject to external potentials with weak gravitational driving towards the bottom of the system, a 1 = a 2 = 0.01ε. In the initial stages, the system evolves in a manner quite similar to what was observed in Figs. 3 and 4, i.e. several demixed horizontal layers are formed. However, as can be seen most clearly from the first of the time series of total density profiles along the bottom row of Fig. 6, in contrast with the previous cases where the total density was nearly constant, one can now see that the layers rich in species 2 are of substantially higher density than the ones rich in species 1. Thus, this is initially a gas-liquid type phase separation, where the liquid-like regions are rich in species 2, while the lower density gas is rich in species 1. As time proceeds, at time t * = 400 we see that the horizontal bands across the middle of the system break-up into drops. This, in turn allows the colloids of species 1 from the neighbouring regions to gather together. The strong inhomogeneity in the density in this region leads to the formation of the third phase, which is a liquid rich in colloids of species 1. As time proceeds there is also a coarsening as drops of the same phase coalesce to form larger drops. Finally, the system reaches a state consisting of the following mi- crophases: Along the bottom of the system, we observe a low density (gas-like) layer of variable thickness that is rich in particles of species 1. Above this, there is a liquid layer rich in species 2, above which are several drops composed of the liquid phase rich in species 1. Next, there are further drops of the liquid phase rich in species 2. All of these liquid drops are surrounded by a gas phase of roughly equimolar composition. Note that because the gravitational driving towards the bottom of the system is rather weak, the sedimentation of the liquid drops towards the bottom and the displacement of the gas phase towards the top is rather slow (not displayed).
To illustrate the influence of the gravitational force, we now consider the same system (i.e. with T * = 0.5, ρ * = 0.3 and x = 0.3), but where the external potential driving the colloids towards the bottom of the system is substantially stronger, with a 1 = a 2 = 0.1ε. In Fig. 7 we display a time series of concentration profiles and the corresponding total density profiles. As expected, the sedimentation is much faster now, as can be seen from the total density profiles along the bottom row of Fig.  7. In the early stages of the dynamics the strong gravitational driving leads to a large flux of both species of colloids towards the bottom, leading to a rapid increase in the density in the lower half of the system. This situation may be termed as "gravity-dominated", since the external field has a dominating impact on the structure of the fluid, while the effect of the intermolecular forces is somewhat weaker. However, even during this stage the particles show their tendency to demix, as evident from the time t * = 100 concentration plot in the upper left panel of Fig. 6. As time proceeds, the higher density liquid phase formed at the bottom of the system by the strong gravitational driving then undergoes liquid-liquid demixing. This occurs initially by the system forming alternating horizontal layers of the two different liquid phases. Subsequently, the thinner layers rich in the minority species 1 colloids break-up into drops. This dynamics is driven by the fact that it lowers the amount of interface between the two liquid phases. The late-time configurations consist of a single layer of the species 2 rich liquid in the lower half of the system, within which are drops of the species 1 rich liquid (recall that the gravitational force is the same for both species, so neither liquid phase sediments within the other). The gas phase at the top of the system is of a roughly equimolar concentration.
One further thing to note, comparing Fig. 6 (weak gravity) and Fig. 7 (strong gravity), is that in the former the wall at the bottom of the system becomes covered by a thin gas-like layer. Even if disturbed, this thin layer is very stable, since the wall is hard and is thus preferentially adsorbed by the gas phase. However, in the case in Fig. 7, the strong gravitational force overcomes the interfacial effects and pushes the liquid phase onto the wall.
The impact of the gravitational field can also be seen in Fig. 8 where we display plots of the height of the center-of-mass of both species as they vary as a function of time. For all times t > 0, the center-of-mass of the majority species has a lower value than that of the minority species. Initially, the center-of-mass of both species decreases as one would expect. However, interestingly, in all cases the center-of-mass of the minority species is a non-monotonic function of time, decreasing to a minimum value, before increasing again. In the case of weaker external field (upper panel), the initial "gravitydominated" regime is less dramatic and the minimum is shifted to later stages. However, if the the time axis is rescaled by a factor proportional to the amplitude of the external field, the positions of the minima are comparable. The center-of mass rise is especially pronounced when the composition of the mixture is well away from x = 0.5 and for the weaker field. In this second stage thermodynamic effects become important, so that the liquid-liquid phase separation prevents the minority component from sedimentating to the bottom.
We now consider the case of an equimolar mixture, i.e. with x = 0.5, sedimenting under the influence of the same strong external field, with a 1 = a 2 = 0.1ε. A time series of local concentration and total density profiles for this case are displayed in Fig. 9. In the same way we saw already in Fig. 7, the stronger external field induces a rapid build-up of the liquid above the lower wall, during the initial "gravity-dominated" period of the sedimentation, leading to the formation of a layer of liquid at the bottom of the system, with the gas above. The rapid increase in the density at the bottom of the system then triggers the liquid-liquid phase separation within the liquid phase. In the bulk fluid, this would result in the formation of two separate large slabs of the two liquid phases, separated by a single interface. However, in the present case, where the fluid is confined by the walls and also pushed to the bottom by the gravitational force, the system is restricted in its ability to coarsen to form just two slabs, as it would in bulk spinodal phase-separation. Instead, we see the system break up into vertical stripes alternating between the two different liquid phases, respectively rich in species 1 and species 2. The interfaces between the stripes are vertical (contact angle of 90 • at the lower wall) because the mixture is symmetric and so the wall does not favour either liquid phase over the other. There is some coarsening in the number of stripes, but the system is unable to reach the minimum free energy state consisting of just one drop of each of the two different liquid phases at the bottom of the system. This is because to reach the minimum the system must surmount a free energy barrier by a diffusion of particles of one species through a region rich in the other species. Such a process is quite unlikely, given the low temperature of the system. We thus believe that the final time t * = 600 plot in Fig. 9 corresponds to a configuration which is very stable with an extremely long (maybe even infinite) relaxation time for reaching the true equilibrium configuration. Finally, we now discuss a case where species 1 experiences a stronger driving towards the bottom of the system than species 2. This, for example, corresponds to the case where the species 1 colloids are made of a more dense material than species 2, and so have a much higher buoyancy mass. In Fig. 10 we present results for a system with the same parameters as in the previous system, with the only difference that now a 2 = 0.01ε, i.e. the amplitude of the external field experienced by species 2 is one order of magnitude less than the one experienced by species 1 (a 1 = 0.1ε). In Fig. 10 we see that the much greater sedimentation rate of species 1 leads to a separation of the two species. The center-of-mass of species 2 is almost constant (not displayed). This is because some of the slowly falling species 2 rich liquid drops are subsequently pushed upwards by the faster sedimenting species 1 rich liquid. In this situation the disparity in the strength of the external driving forces is large enough to overcome any free energy barriers preventing a drop of one liquid phase being dragged through the other. In the final stages, the system forms a state with the following morphology: a slab of liquid rich in species 1 is located along the bottom of the system, whilst drops of the liquid phase rich in species 2 remain above this layer with the gas phase which has a roughly equimolar composition intruding between the drops. In contrast with what we observed in Figs. 7 and 9, these liquid drops are not in close contact with the bottom liquid but instead almost 'float' in the gas phase. It is because the external force on them downward is much weaker and also because the interfacial tension between the liquid phases and the gas phase is less than the interfacial tension between the two liquid phases. Consequently, these drops lower the overall free energy by rising a short distance above the lower liquid to allow the gas to intrude between the two liquid phases.

IV. CONCLUDING REMARKS
In this paper, a two-dimensional binary mixture of partially immiscible colloids dispersed in a continuous solvent medium was studied using DDFT. The mixture is subject to a gravitational field, which pushes the particles down towards the bottom of the system at y = 0, where there is an impenetrable hard-wall. In the absence of the gravitational field, the contact angle of each individual species is 180 • for all temperatures T < T c , as usual for a fluid at a hard wall. However, we have seen that the presence of the other species and the gravity strongly influences the interfacial properties at the wall.
This in turn has a significant effect on the structure and microphase behaviour of the entire system both in equilibrium and out of equilibrium.
A description of all possible behaviours that our model can exhibit during the sedimentation relaxation towards equilibrium is a lengthy and almost impossible task. The vast array of possible configurations found in Ref. 26 using continuation for a mesoscopic model binary fluid on a surface together with the gas phase above, should also be observed in the present model and gives an indication of what drop morphologies can be observed in the present system. Our intention here is not to categorise all possible morphologies but instead to point out some of the most characteristic features of the model and to demonstrate the complexity of its dynamical properties. The systems discussed here differ in temperature, density, composition and/or strength of the gravitational field. Based on these examples, we have illustrated an interplay between the bulk thermodynamics, interfacial (both solid-fluid and fluid-fluid) phenomena and gravity that gives rise to a rich variety of behaviours of the system which can often considerably change over time.
The dynamical properties of colloidal suspensions are often modelled using stochastic Brownian dynamics (BD) computer simulations, based on numerically integrating forwards in time the Langevin equations of motion in Eq. (7). Here, we have instead used an alternative microscopic approach, namely DDFT. This has an advantage that for a given size of the system (and a given discretisation) the numerical requirements are independent of the total fluid density, i.e., systems that correspond to say ∼ 10 2 particles or those corresponding to ∼ 10 4 particles can be handled at the same numerical price. Since BD simulations are too computationally expensive to implement for the systems considered here, we solely focussed on applying the DDFT. The surprisingly good agreement between the BD simulation results and those from the rather simple DDFT treatment of the model found previously for the one-component fluid, 17 gives us confidence that our DDFT approach is also reliable for the binary mixture. Note that we have solely presented results from 2D calculations. Since the external potential is translationally invariant along the horizontal axis, in principle one need only treat the system as effectively in 1D. However, using an approximate DFT in 1D leads to poor results -see e.g. Fig. 9 of our previous paper dealing with a one-component fluid. 17 Below we discuss this issue further.
We should stress that here we have assumed that the hydrodynamic interactions between the colloids can be neglected, both in the Langevin equations of motion (7) and also in the resulting DDFT (15). It would be interesting to investigate their influence on the system behaviour, because in certain situations this can be significant. 27 Hydrodynamic interactions can be included within the DDFT framework -see e.g. Refs. 28 and 29, but we leave this aspect to future work.
Our main observations can be summarised as follows: At high densities and/or high temperatures, the fluid exhibits liquid-liquid demixing but no liquid-gas transition.
In this case, depending on the concentration, we observe three characteristic scenarios: 1. When the concentration of one of the species is low, we observe the formation of liquid droplets rich in the minority species that are dispersed within a slab of the majority species rich liquid. Such a state is preceded by the formation of horizontal layers that break up rather quickly to form the drops. However, the presence of a large number of small drops is thermodynamically unfavourable due to the large free energy cost of all the interfaces. The drops thus have a tendency to coalesce whenever two or more of them approach each other. Alternatively, they may also coarsen via the Ostwald ripening mechanism, but the low miscibility of the mixture can lead to this microemulsion-like density distribution being stable over a very long time.
2. At a somewhat higher concentrations of the minority species (but still for x = 0.5), horizontal stripes also form in the initial stages, but in this situation some of the stripes persist without breaking-up.
3. For the equimolar x = 0.5 case, the morphology of the fluid is striking different. This is due to the symmetry of the mixture and so the confining walls do not direct the spinodal decomposition. The phase separation is initiated in the bulk and leads to vertical rather than horizontal stripes, since the final equilibrium state must be invariant with respect to exchanging the labels on the two different species (i.e. the possibility of a preferential adsorption of either of the two species is excluded by symmetry).
We have also examined the behaviour of the mixture at a state point below the critical point of the liquid-gas transition. In this case, the system separates into three phases: a gas and two liquid phases. To the best of our knowledge, what we present here is the first microscopic study of the dynamics of such a system. We focus in particular on the influence of gravity, displaying results for three examples: 1) the gravitational field is weak and same for both components, 2) the gravitational field is strong and same for both components, and 3) the gravitational field is strong for species 1 but weak for species 2. All cases exhibit a very rich microstructure, not only due to the presence of the three phases but also in the manner the different phases are formed and distributed over different parts of the system, forming different shapes and also in the manner they change over time. Coalescence of regions of one phase is also hindered by the presence of the other phases. The effect of the gravitational driving can be profound, with the center-of-mass of the minority species being a non-monotonic function of time, first decreasing and then, surprisingly, moving up again due to the capillary forces in the system. We conclude with a few general remarks on the nature of the results obtained from a 2D implementation of DDFT -see also the discussion in Ref. 17. In the limit t → ∞, the resulting configuration corresponds to the density profiles for an equilibrium fluid, as obtained from DFT, since from Eq. (15), at equilibrium we have Such equilibrium profiles must, of course, reflect the symmetry of the external potentials since DFT is a statistical mechanical theory, which in principle performs a statistical average over all realisations of the ensemble. This condition should be obeyed no matter what approximation is made for the free energy functional or whether the solution is obtained from a 1D or a 2D calculation.
Our results, where we find density profiles that vary in the x-direction, the direction in which the external potentials are invariant, lead to the obvious question: Why do we observe this symmetry breaking? This is a particularly pertinent question, since the results we find do agree qualitatively with what one would observe in any individual realisation of the ensemble -see e.g. the comparisons for the one-component fluid in Ref. 17.
To address this question, it is important to stress that Eq. (15) is a deterministic theory for a local and instantaneous density distribution, i.e. for the local density averaged over the noise realization. Clearly then, if both the external field and the initial density distribution only vary in the vertical direction, then ρ(x, y, t) can also change only in that direction and the solution of the 2D problem is exactly the same as if the system is treated as a 1D problem. Nevertheless, the perturbation field δ ρ (x, y) added to the otherwise uniform initial density profile may break the symmetry and allows in principle for density profiles that exhibit inhomogeneities even along the x-axis. If the system contains a thermodynamically stable phase, then these density fluctuations quickly disappear. However, inside the spinodal (as in our case), the fluid is thermodynamically unstable and it follows that there exist some density fluctuation modes whose amplitude will grow exponentially with time (at least at early times). The length scale corresponding to the fastest growing mode can be calculated for the bulk fluid from the dispersion relation 7,14 resulting from linearizing Eq. (15). This length scale is typically 10σ-20σ, depending on the state point and is in good agreement with the length scale observed at early times in our density profiles.
Another important point to be emphasized is that DDFT is a theory that is formulated in the canonical ensemble, with a fixed number of particles N = N 1 + N 2 . In contrast, DFT is formulated in the grand canonical ensemble, which is coupled to a reservoir that fixes the chemical potentials µ i and with which the system can exchange particles. In the limit t → ∞ DDFT reduces to the Euler-Lagrange equation (21) subject to the additional condition dr(ρ 1 (r) + ρ 2 (r)) = N . Thus, even in the equilibrium limit, the DFT and DDFT results do not generally coincide, unless a proper thermodynamic limit (i.e., N → ∞) is taken. For any finite values of N (and volume V ), though, the finite size of the system acts as an external constraint, which plays a role somewhat akin to a stabilizing field, for interfacial structures that are metastable in the thermodynamic limit. 30 A well-known example is a nucleating drop of liquid in a supersaturated vapor. 31 All these points hold irrespective of whether the free energy functional used in the theory is exact or approximate. Clearly, for any physically relevant system one always deals with an approximate functional and as shown in our previous work for the one-component fluid, 17 ρ(y) is more accurate when obtained from 2D DFT after averaging ρ(x, y) over the x-direction (the direction parallel to the wall) rather than directly from the 1D DFT. We believe this is due to the fact that within our meanfield treatment some of the fluctuation contributions to the free energy, such as capillary waves, are inevitably ignored. This issue is also nicely discussed by Reguerra and Reiss, 32 who studied what fluctuation effects are underestimated in the standard approximations used in DFT when used to study drops of liquid confined in a cavity. As discussed in Ref. 17, within 2D (D)DFT one relaxes some of the constraints made in a 1D (D)DFT calculation, and thus the resulting averages are made over a larger set of microstates.
This situation is akin to the following simple mechanical example: Consider a cylindrical elastic rod with its tip in contact with the floor and its axis of symmetry perpendicular to the floor. If the rod is gently pressed down on the opposite tip with a force parallel to the axis of symmetry, then the rod will bend slightly outwards in the middle, breaking the symmetry. The direction in which it bends is arbitrary. If the experiment is repeated numerous times and an average over all the positions of the rod is calculated, then the average state will consist of an unbent rod, since all the different realisations compensate one another. Thus, such a series of experiments will give the same "mean" of the rod shape as a "1D" consideration but contains more information about the shape fluctuations of the rod. Accounting for the states corresponding to when the rod is bent in a particular direction is, in some sense, the way the present system is treated by the approximate 2D (D)DFT.