Equilibrium and non-equilibrium cluster phases in colloids with competing interactions

Abstract

The phase behavior of colloids that interact via competing interactions – short-range attraction and long-range repulsion – is studied by computer simulation. In particular, for a fixed strength and range of repulsion, the effect of the strength of an attractive interaction (ε) on the phase behavior is investigated at various colloid densities (ρ). A thermodynamically stable equilibrium colloidal cluster phase, consisting of compact crystalline clusters, is found below the fluid–solid coexistence line in the ε–ρ parameter space. The mean cluster size is found to linearly increase with the colloid density. At large ε and low densities, and at small ε and high densities, a non-equilibrium cluster phase, consisting of elongated Bernal spiral-like clusters, is observed. Although gelation can be induced either by increasing ε at constant density or vice versa, the gelation mechanism is different in either route. While in the ρ route gelation occurs via a glass transition of compact clusters, gelation in the ε route is characterized by percolation of elongated clusters. This study both provides the location of equilibrium and non-equilibrium cluster phases with respect to the fluid–solid coexistence, and reveals the dependencies of the gelation mechanism on the preparation route.

---

1. Introduction

In experimental colloid physics, colloid–polymer mixtures have been used as a model system to study various aggregation phenomena such as physical gelation¹ and crystallization.² The advantage of colloid–polymer mixtures is that attractive and repulsive parts of the interaction can be independently tuned. The attraction between the colloids is purely entropic: the so-called depletion interaction.³ The range and strength of the depletion-induced attraction can be tuned by varying the molecular weight of the polymers and their concentration, respectively. The range of attraction is dictated by the radius of gyration of the polymer in the dilute regime, and hence this range can be made very short (<1% of the diameter of the colloid) using a low-molecular weight polymer. An additional repulsive interaction of the colloids can be induced by an electrostatic charge of the colloids. These competing effects of long-range repulsion and short-range attraction dictate the phase behavior of the colloids.

Several experimental, theoretical and numerical studies have been conducted on this type of colloidal model system. An important finding is the so-called cluster phase. In this phase, particles self-assemble into aggregates. As the particles are attractive at short range and repulsive at long range, beyond a certain size, addition of a particle to the cluster contributes more repulsion than attraction, hence growth is energetically not favored. The first experimental evidence for the cluster phase in colloids with competing interactions was reported by Segre et al.⁴ An equilibrium theory developed by Groenewold and Kegel explained the existence of the equilibrium cluster phase and established its necessary conditions.⁵ This theory was confirmed by Stradner et al., who reported that stable clusters were found in a lysozyme protein solution under salt-free or low-salt conditions.⁶ These authors also reported an equilibrium cluster phase in a suspension of charged poly(methyl)methacrylate (PMMA) colloids with polystyrene as the depleting polymer. The size of these clusters was found to increase linearly with the volume fraction of the colloids, as predicted by the theory.⁵ Campbell et al. reported the formation of clusters in a highly charged suspension of PMMA colloids at a dilute volume fraction of colloids.⁷ As the volume fraction was increased, they observed that the clusters aggregated into a network of Bernal spirals, which triggered gelation at sufficiently high attractive strength. Zhang et al. systematically investigated the effect of strength of attraction on the cluster size and cluster morphology at a fixed colloid density.⁸ The authors reported that, in the range of attraction strength from 5.5 to 16 k_BT, the cluster size decreased with increasing attraction strength. They attributed this observation to a higher nucleation rate at higher attraction strength leading to a large number of smaller nuclei, and resulting in non-equilibrium effects in the dynamics of rearrangement.

Inspired by these experimental findings, several computer simulation studies were reported on the colloids with competing interactions. For several combinations of strength and range of repulsion, Mossa et al. calculated the ground-state energies of clusters of various sizes to provide evidence of an equilibrium cluster phase.⁹ Fortini et al. reported the phase behavior of charged colloid–polymer mixtures for various ranges of attraction in the limit of high electrostatic screening (where the range of repulsion ≈0.02 times the diameter of the colloid).¹⁰ In this parameter limit, the authors did not find any cluster phase. Charbonneau and Reichman showed, using a potential similar to the one studied in the present work, that in a two-dimensional system microphase separation occurs below the fluid–solid coexistence.¹¹ This microphase separation corresponds to a coexistence of low-density and high-density phases. Sciortino et al. also reported a similar study with a Lennard-Jones (20-10) and Yukawa repulsive potential.¹² At low temperature (or large ε) they observed an extended one-dimensional cluster growth.

Despite the experimental and numerical studies mentioned above, a comprehensive phase diagram showing the region of stability for the cluster phase is not available. Yet, the location of the cluster phase is necessary to understand the pathways of gelation and crystallization. Here we present the phase diagram for a model system with very short ranged attraction and long ranged repulsion. We compute the fluid–solid coexistence of the model and identify the location of equilibrium and non-equilibrium cluster phases. We infer from the phase diagram that the gelation can be reached via an equilibrium cluster route or a non-equilibrium cluster route, with two structurally dissimilar particulate gels.

2. Model & simulation methods

The pair interaction between the colloids, manifested by the combination of depletion attraction and electrostatic repulsion, is modeled as a combination of the 100-50 Lennard-Jones potential and the screened Coulombic potential as below

In this equation ε is the strength of attraction, σ is the diameter of the colloid, A is the prefactor in the screened Coulombic potential which is related to the surface potential of the colloids¹³ and ξ is the Debye screening length. The large exponents in the LJ potential ensure that the range of interaction is very short. We fix A = 2 k_BT and ξ = 1.794σ in all the calculations reported in this paper. For a typical colloid size of σ = 100 nm, a Debye screening length of ξ = 179 nm corresponds to a 1 : 1 salt concentration of 3 μM in water. In the following, energy is expressed in units of k_BT, where k_B is the Boltzmann constant and T is the temperature. Length scales are expressed in units of σ. The unit of time is considered as , where m is the mass of the colloid. Here we take the mass to be unity, which effectively sets the time scale. Fig. 1 plots the overall potential for several values of ε, and it can be seen that the attractive interactions are short ranged (<0.025σ) and the electrostatic part is long ranged (>3σ).

Fig. 1 Plot of inter-colloid interactions for ε = 2 (black), ε = 5 (red) and ε = 10 (blue) from eqn (1).

We perform molecular dynamics simulations in the NVT ensemble for phase-diagram prediction and Monte Carlo simulation in the NPT ensemble in the Gibbs–Duhem integration. The dynamics of the system is integrated in molecular dynamics simulations with constant volume and constant number of particles (NVT) using either N = 1372 or 256 colloids for densities ranging from 0.025 to 1 and ε ranging from 1 to 20. The density is defined as the number of colloids per unit volume (in terms of σ³). We use a cubic periodic box in the simulation. We employ a cut and shifted potential with a cut-off set at 5. The equations of motion are integrated using the velocity Verlet algorithm with a time step of 0.001. The temperature is controlled by using a thermostat with a time constant of 0.2, rendering the dynamics in the high friction limit. In all the simulations, the reduced temperature T is kept unity. Equilibration and production runs are performed with 10⁸ MD steps each. We use the following criterion to define a cluster: if the distance between two colloids is within a cluster cut-off, which is fixed at 1.025, they belong to a cluster. The minimum of the potential is at around r = 1.014 and the distance at FWHM is around 1.03, though the potential is not a Gaussian function. These numbers slightly vary as a function of ε. Therefore, we chose 1.025 as the cut-off for the cluster criterion. As the potential is very steep near the minimum, the choice of cut-off distance does not make substantial difference for the analysis of the results.

To obtain the fluid–solid coexistence lines, we employ the Kofke's Gibbs–Duhem integration method.¹⁴ This method allows one to start from an a priori known coexistence point, and integrate along the coexistence line using the Gibbs–Duhem equation (GD) or alternatively, the Clapeyron equation. Since we only know the hard-sphere coexistence in advance, this procedure is done in two steps. In the first step, starting from hard-sphere initial conditions, the integration is carried out by slowly increasing the repulsive interaction, i.e. the parameter A in the hard-sphere Yukawa potential of the form

The equation to integrate is then derived from the Gibbs–Duhem equation,

dμ = vdP − sdT + λ₁dA (3)

In this equation, μ is the molar chemical potential, v is the molar volume, P is the pressure, s is the molar entropy and λ₁ is a coupling parameter. For two phases – fluid and FCC crystal – to be in equilibrium at constant temperature, it directly follows from eqn (3) that

where Δ denotes the difference in the properties (λ₁ and v) between fluid and solid phases, i.e. Δv = v_l − v_s and Δλ₁ = λ_1,l − λ_1,s.

From eqn (3), it follows that

where the angular brackets denote an ensemble average at constant T, P and ξ. Note that the reduced temperature T is unity. Eqn (4) is integrated from A = 0 to 2 with an increment of 0.2, using the initial conditions of the hard sphere solid–liquid coexistence with a coexistence pressure β_P = 11.7, leading to a density of the fluid, ρ_f, = 0.943 and a density of the FCC crystal, ρ_s = 1.041.¹⁵

Fig. 2 shows the result of the GD-integration. From this figure, we can observe that upon increasing the repulsive interactions the fluid–solid coexistence shifts to lower densities, and the coexistence region becomes narrower. A similar trend was observed in previously reported studies.^16,17

Fig. 2 Gibbs–Duhem integration from the hard sphere (A = 0) to Yukawa repulsive potential (A = 2).

At A = 2, the solid–liquid coexistence conditions are: density of the fluid, ρ_f, = 0.9262, density of the FCC crystal, ρ_s, = 0.9722 and coexistence pressure, β_P, = 38.1. Using these values as the initial conditions, a second integration step is performed from the purely repulsive Yukawa potential to the full potential including the (100-50) LJ and repulsive Yukawa potential (eqn (1)) by slowly increasing ε from 0 to the desired value. Following a similar derivation as in the first step, we obtain the analog of the Clapeyron equation as

andwhere again the reduced temperature T is assumed unity. Eqn (7) is integrated from ε = 0 to 12 with an increment of 1.

Note that the two-stage integration of the Gibbs–Duhem equation described above is not sensitive to the path of the integration. In fact, we investigated several combinations of integration routes. In the first case, we started with hard-sphere colloids and simultaneously increased attractive and repulsive potential to a desired value in a single integration step. In the second case, the integration was performed in two stages from the hard-sphere to LJ and from LJ to Yukawa repulsion. We found a good match of results (fluid–solid coexistence curves) between all these three routes.

3. Results and discussion

3.1 Phase-diagram

We perform straightforward NVT-MD simulations for a range of densities ρ and interaction strengths ε, at fixed A = 2 and ξ = 1.794σ. In all the simulations, the initial coordinates of the particles are obtained from an FCC lattice with the lattice spacing corresponding to the desired density, and the initial velocities were obtained from the Maxwell–Boltzmann distribution. The simulations are analyzed based on their cluster distribution, as well as their dynamics. In addition, we perform a GD integration to compute the solid–liquid coexistence as a function of ε. The results from straightforward NVT-MD simulations and the Gibbs–Duhem integration are summarized in the parameter space of ε–ρ in Fig. 3. At low values of ε and ρ, we observe a fluid phase for which the cluster size distribution (CSD) function decays monotonically with the cluster size. Here, the clusters are transiently formed and broken, suggesting a fluid-like nature of this phase. In the density range of ρ < 0.2, an increase of ε leads to the formation of an equilibrium cluster phase. In this regime, the CSD shows a clear peak at a finite cluster size. Such clusters are in equilibrium with individual colloids (monomers). The monomer fraction, x₁, obeys the condition from the theory of self-assembly as¹⁸
where the depth of the well U_w ≈ ε − A can be inferred from the potential depicted in Fig. 1. Fig. 4 shows the CSD corresponding to an equilibrium cluster phase at a fixed attraction strength, ε = 7 k_BT for various colloid densities. The mean cluster size increases with the density of the colloids, in line with theoretical predictions.⁵ The size distribution broadens as the density is increased. Fig. 5 shows the CSD for an equilibrium cluster phase at a fixed density, ρ = 0.075, for various attraction strengths, ε. Fig. 4 and 5 suggest that the equilibrium cluster phase can be reached via either the ρ-route or the ε-route.

Fig. 3 Phase diagram of the LJ + Yukawa colloid showing fluid, equilibrium cluster, non-equilibrium cluster and gel phases, and fluid–FCC coexistence.

Fig. 4 Cluster size distribution for ε = 7 k_BT at various densities of colloids.

Fig. 5 Cluster size distribution for ρ = 0.075 at various attraction strengths, ε.

A further increase in ε (ε > 10, ρ < 0.2) leads to the formation of a non-equilibrium cluster phase. We identify this phase from a CSD showing multiple maxima (may be discrete) and a monomer fraction that is either zero or

A non-equilibrium cluster phase is also present at moderate ρ > 0.2 and low ε < 4.

At moderate and high densities, ρ > 0.2, the colloids can form a gel phase when the attraction is sufficiently strong. The gel is characterized by a CSD that shows one big cluster with a size close to the total number of particles (a percolated cluster).

The fluid–FCC coexistence lines obtained from the GD-integration are also shown in Fig. 3. Upon increasing ε from zero, the coexistence region widens up as the attractive interaction begins to dominate. The coexisting crystal phase becomes denser, while the fluid phase becomes more dilute. Note that the cluster phase and part of the gel phase fall outside the fluid–solid coexistence region. This means that the system will form a gel before crystallizing, and hence it is very likely that crystallization is kinetically hampered in this kind of systems. In addition, whereas a straightforward NVT simulation would get stuck in a metastable state, the Gibbs–Duhem integration allows for the evaluation of the true phase-diagram.

In the following we compare our results with previous simulations by Toledano et al.¹⁹ and others.^9,11 In their study, parameters were ξ = 2σ and A = 0.2ε and the dimensionless temperature (T*) was defined as , where ε is the strength of attraction. Simulations were done for different sets of T* and ϕ, volume fraction of colloids. With the temperature defined in this way, the strength of repulsion (A) also varies when T is varied. In contrast, we treat the system with variable attraction strength while keeping the strength of repulsion constant in all our simulations. Still, the data corresponding to T* = 0.1 in Toledano et al. correspond to our data points at ε = 10 k_BT and A = 2 k_BT. The value of ξ = 1.794 used in our simulation is only 10% lower than the value used in ref. 19. Table 1 shows comparison of our results with that of ref. 19. Despite the different value of ξ used in our simulation, the comparison of the data reasonably agrees with the data of Toledano et al.¹⁹

Table 1 Comparison between the present simulation and data from ref. 19

Data of Toledano et al.^19			Data from this work
Volume fraction, ϕ (at T* = 0.1)	Density, ρ (ε = 10 kBT)	Phase	Density, ρ (at ε = 10 kBT)	Phase
0.04	0.076	Cluster – slightly elongated	0.075	Elongated clusters
0.1	0.19	Cluster – elongated	0.2	Percolated gel
0.16	0.3	Percolated gel	0.3	Percolated gel

---

For the particular potential used in this work, and for the range of parameters studied, we do not find spontaneous formation of any layered or columnar structures. Previous reports (ref. 9–12 of the manuscript) that employed potentials and parameters similar to ours also did not report such modulated structures. Our model is based on depletion-induced attraction and screened electrostatic repulsion that can be experimentally realized. While de Candia et al.²⁰ reported columnar and lamellar phases, in their model, the range of attraction was 0.2σ, whereas in our simulation it is only 0.02σ, 10 times shorter in range. Wu and Cao²¹ calculated ground state energies of several structures such as lamella, cylinder and spherical domains in a continuum approach, and concluded that the stability of these phases depends on the relative ratio of surface free energy to electrostatic energy.

3.2 Structure of the clusters

The shape of the clusters in the equilibrium phase (6 < ε < 8, ρ < 0.2) is found to be rather compact and small (∼10 particles per cluster). As the clusters grow in size, they initially show a quasi-linear shape due to the long repulsive tail, i.e., ξ = 1.794σ. Subsequently, the particles rearrange in a second step to more compact clusters. This rearrangement is possible because the attraction strength corresponding to the equilibrium cluster phase ranges from 6–8 k_BT, which translates in an effective attraction (well-depth, U_w) range from 4–6 k_BT. These barriers can be occasionally overcome by thermal energy in the time scale of the simulations. Fig. 6 shows some representative snapshots of equilibrium clusters.

Fig. 6 Snapshots of representative equilibrium clusters at ε = 7 k_BT at various densities: (a) ρ = 0.1, (b) ρ = 0.125 and (c) ρ = 0.15.

In contrast, clusters in the non-equilibrium cluster phase are more elongated and resemble the Bernal spiral. Such structures were indeed found in experiments close to gelation.⁷ As the non-equilibrium clusters form at much higher attraction strengths, 10–20 k_BT, the particle rearrangement becomes very slow, leading to kinetic trapping. As the cluster sizes in the non-equilibrium regimes are large compared to the equilibrium clusters, it is useful to calculate the fractal dimensions of the non-equilibrium clusters.

The fractal dimension is defined via the radius of gyration (R_g) of the cluster as

where N is the number of particles in the cluster and D_f is the fractal dimension. For solid macroscopic one-dimensional, two-dimensional and three-dimensional structures, the fractal dimensions are simply D_f = 1, 2 and 3, respectively. For branched structures the fractal dimension can be fractional. We compute the R_g values of clusters with various N, and extract D_f from a plot of ln(R_g) vs. ln(N). Clusters consisting of 10 to 70 colloids are considered for the R_g calculation.

From the data given in Table 2, we can infer that the non-equilibrium clusters are of elongated shape, as the D_f values are close to that of a one-dimensional solid object. It should be noted that the R_g analysis was done using finite sized clusters in the size range of 10–70. Fig. 7 shows some snapshots of non-equilibrium clusters.

Fig. 7 Snapshots of representative non-equilibrium clusters at ρ = 0.15 at various strength of attractions: (a) ε = 10 k_BT, (b) ε = 15 k_BT and (c) ε = 20 k_BT.

Table 2 D _f-values for some of the state points in Fig. 3 for ρ = 0.15

ε (kBT)	D f
10	1.12
15	1.28
20	1.38

---

In Fig. 8 we plot the radial distribution function (RDF) of the colloids for two different densities at ε = 7 k_BT. These two state points correspond to the equilibrium cluster phase. For the higher density ρ = 0.2 the RDF indicates that the clusters consist of particles arranged in a crystal-like structure. This is evident from the long-range peaks in the RDF. Clearly these long-range peaks are absent in the RDF corresponding to lower density, implying smaller clusters. Visual inspection of the configurations of the clusters also supports this observation (see Fig. 9). Note that the cluster–cluster interaction is strongly repulsive, leading to a fluid of clusters.

Fig. 8 Radial distribution function at ε = 7 k_BT for densities ρ = 0.075 and 0.2.

Fig. 9 Configurations of the clusters at ε = 7 k_BT for densities ρ = 0.075 (a) and 0.2 (b).

3.3 Dynamics of gelation

The data points given in the phase diagram (Fig. 3) are obtained from independent NVT simulations that are performed at different sets of (ρ, ε). This phase-diagram can be useful in predicting the mechanism of gelation for different experimental conditions. The gelation can be reached either via the ρ-route or the ε-route. In the ρ-route, for a fixed polymer concentration (fixed ε), more colloids are added to the suspension to increase the colloid number density (volume fraction). In this preparation, the gel is likely to be formed from densely packed compact crystalline clusters. This mechanism is commonly denoted as a cluster glass transition.²² In contrast, in the ε-route, for a fixed colloid density (fixed ρ), more polymers are added to the suspension to increase the strength of attraction (epsilon). In this method, the gel is likely to be comprised of percolated networks of elongated non-equilibrium clusters. This mechanism is known as percolation. Thus, depending on the route we may expect different mechanisms of gelation in the same colloid–polymer model system.

The dynamics of gelation can also be captured by the intermediate scattering function, F(k, t), with wavevector k as a function of time t. The function F(k, t) is shown in Fig. 10a and b for different attraction strengths at a constant density of 0.05. Fig. 10a is obtained for an attraction strength of 5 k_BT for various kσ values (0.36 to 1.54), corresponding to the fluid phase according to the phase diagram in Fig. 3. Fig. 10b shows F(k, t) for an attraction strength of 10 k_BT, and dynamical arrest on the time scale of the simulation is observed suggesting non-equilibrium nature of clusters. The function F(k, t) for a higher density of 0.2 and at different attraction strengths, ε = 5 and 10 k_BT, is depicted in Fig. 10c and d. Fig. 10c shows a rapid decay for all k values (0.58 to 1.63), suggesting a fluid of aggregates. Whereas Fig. 10d shows kinetic arrest on the time scale of the MD simulation for an attraction strength of 10 k_BT, where gelation due to percolation of linear-clusters occurs.

Fig. 10 F(k, t) from the MD simulation for different kσ values. (a) ρ = 0.05, ε = 5 k_BT, (b) ρ = 0.05, ε = 10 k_BT, (c) ρ = 0.2, ε = 5 k_BT and (d) ρ = 0.2, ε = 10 k_BT.

Our work clearly shows that for the dynamics of gelation, the structure and other physical properties of the gel can depend on the route of gelation. These conclusions are in agreement with the recent experimental study of Zhang et al.⁸ In the experimental report of Zhang et al.,⁸ it has been shown that at low attraction strengths, crowding of crystalline clusters leads to ‘microcrystalline’ gels, as shown in Fig. 4a of ref. 8. The experimental conditions correspond to ε = 5.5 k_BT and ρ = 0.69. Similarly, for these parameters, we find in our simulation (Fig. 3) gelation due to such crowding of compact clusters. Segre et al.⁴ reported a gelation mechanism that involved kinetic arrest due to crowding of clusters, analogous to the glass transition. In another set of experiments, Zhang et al. have shown that for ε = 8.8 k_BT and ρ = 0.38, percolation of elongated clusters led to gelation. State points close to this experimental condition in our computed phase diagram also show such percolation driven gelation. Of course, while the gelation might show a history dependence, we realize that the underlying equilibrium state will be independent of the route.

Zhang et al.⁸ also reported that the mean cluster size decreased with increasing attraction strengths in their experiments involving colloids with competing interactions. This conclusion was based on experimental data for ε values of 5.5, 8.8 and 16 k_BT at a density of 0.114 (or ϕ = 0.06). Zhang et al.⁸ argued that the equilibrium theory cannot be used to explain this trend because of non-equilibrium effects such as kinetic trapping. Our simulation data for these parameters represent fluid, equilibrium cluster and non-equilibrium cluster phases, respectively. In Fig. 5 we report that the equilibrium cluster size increases with the strength of attraction, according to the equilibrium cluster theory reported in ref. 5. The range of attraction strengths in Fig. 5 is chosen between 7 and 8.5 k_BT, which corresponds to the relevant range where equilibrium clusters are found. Within this parameter range, we observe an increase in the cluster size with the strength of attraction. We therefore expect that one would observe the predicted trend if experimental conditions are restricted to the equilibrium cluster phase.

4. Conclusion

The computational study reported in this paper demonstrates that for colloids interacting via competing short ranged attraction and long ranged repulsion, both equilibrium and non-equilibrium cluster phases with well defined internal structures are present outside the fluid–FCC coexistence region in the ε–ρ parameter space. While the equilibrium cluster phase contains small and compact clusters, non-equilibrium clusters are large and elongated. The very presence of equilibrium and non-equilibrium clusters influences the mechanisms of gelation. The gelation mechanism can follow either a cluster glass transition or a percolation route, leading to gels of different physical properties such as rheological properties. The study also explains why it is difficult to crystallize colloids under competing interaction regimes, mainly because of intervening cluster phases and gel formation.

Acknowledgements

This work is part of the research programme VICI 700.58.442, which is financed by the Netherlands Organization for Scientific Research (NWO).

References

1. P. J. Lu, J. C. Conrad, H. M. Wyss, A. B. Schofield and D. A. Weitz, Phys. Rev. Lett., 2006, 96, 028306.
2. E. H. A. de Hoog, W. K. Kegel, A. van Blaaderen and H. Lekkerkerker, Phys. Rev. E: Stat., Nonlinear, Soft Matter Phys., 2001, 64, 021407.
3. S. Asakura and F. Oosawa, J. Chem. Phys., 1954, 22, 1255.
4. P. N. Segre, V. Prasad, A. B. Schofiled and D. A. Weitz, Phys. Rev. Lett., 2001, 86, 6042.
5. J. Groenewold and W. K. Kegel, J. Phys. Chem. B, 2001, 105, 11702.
6. A. Stradner, H. Sedgwick, F. Cardinaux, W. C. K. Poon, S. U. Egelhaaf and P. Schurtenberger, Nature, 2004, 432, 492.
7. A. I. Campbell, V. J. Anderson, J. S. van Duijneveldt and P. Bartlett, Phys. Rev. Lett., 2005, 94, 208301.
8. T. H. Zhang, J. Klok, R. H. Tromp, J. Groenewold and W. K. Kegel, Soft Matter, 2012, 8, 667.
9. S. Mossa, F. Sciortino, P. Tartaglia and E. Zaccarelli, Langmuir, 2004, 20, 10756.
10. A. Fortini, M. Dijkstra and R. Tuinier, J. Phys.: Condens. Matter, 2005, 17, 7783.
11. P. Charbonneau and D. R. Reichman, Phys. Rev. E: Stat., Nonlinear, Soft Matter Phys., 2007, 75, 050401.
12. F. Sciortino, P. Tartaglia and E. Zaccarelli, J. Phys. Chem. B, 2005, 109, 21942.
13. E. J. W. Verwey and J. T. G. Overbeek. Theory of Stability of Lyophobic Colloids, Dover Publications, Mineola, NY, 1999.
14. A. Z. Panagiotopoulos, Mol. Phys., 1987, 61, 813.
15. W. G. Hoover and F. H. Ree, J. Chem. Phys., 1968, 49, 3609.
16. F. El Azhar, M. Baus, J. -P. Ryckaert and E. J. Meijer, J. Chem. Phys., 2000, 112, 5121.
17. A.-P. Hynninen and M. Dijkstra, Phys. Rev. E: Stat., Nonlinear, Soft Matter Phys., 2003, 68, 021407.
18. J. N. Israelachvili, Intermolecular and Surface Forces, Academic Press, London, 1992.
19. J. C. F. Toledano, F. Sciortino and E. Zaccarelli, Soft Matter, 2009, 5, 2390.
20. A. de Candia, E. del Gado, A. Fierro, N. Sator, M. Tarzia and A. Coniglio, Phys. Rev. E: Stat., Nonlinear, Soft Matter Phys., 2006, 010403.
21. J. Wu and J. Cao, Phys. A, 2006, 371, 249–255.
22. A. Ikeda, L. Berthier and P. Sollich, Phys. Rev. Lett., 2012, 109, 018301.

---