The effects of polydispersity and metastability on crystal growth kinetics

The effects of polydispersity and metastability on crystal growth kinetics


We investigate the effect of metastable gas-liquid (G-L) separation on crystal growth in a system of either monodisperse or slightly size-polydisperse square well particles, using a simulation setup that allows us to focus on the growth of a single crystal. Consistent with experiment and theoretical free energy considerations, we find that, inside the metastable binodal, a layer of the gas phase ‘coats’ the crystal as it grows. Crucially, the effect of this metastable G-L separation on the crystal growth rate is qualitatively altered by a very small degree of polydispersity as compared to the monodisperse case, suggesting a highly local fractionation process which is facilitated by the gas layer. Our results show that polydispersity and metastability, both ubiquitous in soft matter, must be considered in tandem if their dynamical effects are to be understood.

82.70.-y, 64.75.Gh, 64.60.-i

I Introduction

Substances in the category of ‘soft condensed matter’ (colloidal suspensions, polymers, proteins etc.), as well as having widespread industrial and medical importance, are interesting in their own right because they exhibit physics analogous to that of simpler molecular or atomic systems, but on much longer time- and length-scales Prasad et al. (2007). The analogy stems from the importance of thermal fluctuations in the microscopic dynamics of the constituents, which means that the same thermodynamic and statistical mechanical approaches can be applied in both cases Onsager (1949). Nonetheless, the study and exploitation of soft matter presents a number of challenges, two of which are the focus of the present work.

Firstly, a collection of mesoscopic particles which are nominally the same will almost certainly be polydisperse, i.e. will exhibit a distribution in properties such as size, charge, shape or chemical makeup. This is in contrast to e.g. a molecular fluid of pure water in which, in a strict sense, every molecule is identical. Efforts by various workers to describe the thermodynamics of polydisperse systems Evans et al. (1998); Evans (2001); Sollich and Wilding (2010, 2011); Fasolo and Sollich (2004); Wilding and Sollich (2004); Fantoni et al. (2006); Bartlett and Warren (1999) have resulted so far in a reasonably firm understanding of the effects of a mild degree of polydispersity on the equilibrium phase diagrams of hard-spheres and related systems. Beyond this, polydispersity remains a pervasive but poorly understood complicating factor in soft matter physics – in particular, its effects on phase ordering kinetics as actually enacted in real systems are unclear Zaccarelli et al. (2009); Williams et al. (2008); Martin et al. (2003, 2005); Schöpe et al. (2007); Liddle et al. (2011), and quantitative theoretical work thereon is almost non-existent B. Warren (1999); Evans and Holmes (2001).

Secondly, and in common with simpler systems, soft matter’s path towards equilibrium may be influenced by the presence of metastable states. These are the result of local minima in the free energy landscape which, although they do not fully minimise the free energy of the system (and therefore do not appear in the equilibrium phase diagram), may be encountered as an intermediate stage. They may be long-lived, especially if the system must overcome some free energy barrier in order to reach its global minimum. A prime example is the metastable gas-liquid (G-L) separation of attractive particles, which may subsequently progress to equilibrium crystal-fluid coexistence Liddle et al. (2011). The influence of metastability on phase ordering is a problem of general importance, with particular application in e.g. colloid-polymer mixtures and protein crystallisation Anderson and Lekkerkerker (2002); Dixit and Zukoski (2000); Haas and Drenth (2000, 1999).

In this paper, we perform simulations which focus on crystal growth (as distinct from nucleation) in a model colloidal system exhibiting the metastable G-L coexistence described above. By varying the interaction strength and size-polydispersity, we study the effects of metastability and polydispersity on the crystal growth dynamics. For our parameters, we find that the metastable G-L separation results in a gaseous layer coating the growing crystal, an observation which we relate to theoretical free energy curves and existing experimental work. When polydispersity is present, this layer can act to speed up crystal growth, in contrast to the monodisperse case in which crystal growth is slowed by the presence of the gas layer. Interpreted as a dynamical phenomenon, this suggests that fractionation at the crystal interface, relying on self diffusion, is facilitated by the gas layer. Our results shed light on the influence of polydispersity and metastability on the crystal growth process, demonstrating the complex way in which these factors can interact.

The structure of the paper is as follows. In Section II we describe the simulation apparatus and the parameters used, with reference to the equilibrium monodisperse phase diagram of our system. In Section III we outline two pertinent theoretical aspects of the work, viz. the free energy landscape of the system, and the diffusive growth of split interfaces of the kind exhibited in our simulations. In Section IV we present the simulation results and discuss the effects of metastability, polydispersity, and the two combined on the crystal growth process. Finally, in Section V, we outline and briefly test a possible explanation for the complex interaction between metastability and polydispersity, in the course of which we measure fractionation (de-mixing) associated with crystal formation, and present novel findings relating to local size correlations in polydisperse systems. We conclude and motivate future work in Section VI.

Ii Simulation

ii.1 Kinetic Monte Carlo algorithm

The simulation is built on a Kinetic Monte Carlo (KMC) algorithm, in which the available trial moves are limited to small stochastic ‘hops’, representing the inertia-free Brownian motion of colloidal particles dispersed in a solvent. Therefore, in contrast to equilibrium MC methods, KMC is suitable for studying the dynamics of phase transitions since unphysical moves such as cluster rearrangements, particle resizing etc., which are helpful in speeding up equilibration, are not used. In common with comparable Molecular or Brownian Dynamics methods, hydrodynamic interactions are neglected. Further details of our implementation are available in Ref. Williamson and Evans (2012).

ii.2 Parameters, crystal template

The system studied consists of spherical particles with diameters drawn from a Bates (pseudo-Gaussian) distribution of polydispersity (defined as the ratio of the standard deviation to the mean) (‘monodisperse’), or 1. The mean hard core diameter sets the length unit, and the time unit is the mean time taken for a free particle of diameter to diffuse a distance equal to its own diameter. The hard particle cores are surrounded by pairwise additive square wells of depth and range . The interaction range between specific particles and depends on the particle diameters in a scalable fashion Williamson and Evans (2012); Evans and Holmes (2001):


with . While the range parameter is the same for all particles, Equation 1 shows that the actual range between a given particle pair depends multiplicatively on the size of their hard cores.

The parent volume fraction (hereafter also referred to as ‘concentration’) is . In the monodisperse limit, the values of and are such that the equilibrium state is a coexistence of crystal and gas Liu et al. (2005). However, by choosing (and therefore the effective temperature ) appropriately, the system’s state point can be made to lie either inside or outside the metastable gas-liquid binodal, as shown in FIG. 1 (reproduced from Ref. Liu et al. (2005)). Given that the critical effective temperature corresponds to , we use a well depth of inside the G-L binodal and outside, where the energy unit is . In soft matter, one is often free to vary the interaction strength in this way, for instance by adding more or less polymer into a colloid-polymer mixture.

Figure 1: Phase diagrams reproduced with permission from Ref. Liu et al. (2005) for (a) , (b) and (c) in the density-temperature plane. With our choice of units, and . As decreases, the gas-liquid coexistence becomes metastable with respect to crystal-gas separation. The relevant diagram for our system is (c). The points corresponding to the coordinates used in our study are indicated with hashed circles: point 1 is ‘inside the G-L binodal’ or ‘gas-mediated’ and point 2 is ‘outside the binodal’ or ‘gas-free’. The crystal, the (stable) gas, and the liquid side of the metastable G-L binodal are marked.

Since our aim is to study the process of crystal growth, not nucleation, we introduce crystallisation artificially by placing a template of immobilised particles, arranged in the (100) face of the FCC lattice, at one end of the simulation cell. The simulation cell is long in the axis (), and relatively short in the (periodic) and dimensions (), so that the growth of the templated crystal along the axis can be studied for a long period of time. The template’s lattice parameter is set to match the volume fraction of the equilibrium crystal Liu et al. (2005). Increasing the and dimensions of the lattice had no measurable effect on the results presented. Crystalline particles are identified by constructing bond order vectors from the spherical harmonics ten Wolde et al. (1996); Williams et al. (2008): Particle is flagged crystalline if , where the sum is over the neighbours within of particle , being the average separation between particle and its closest 6 neighbours. The result of the crystal identification algorithm can be seen in FIG. 11.

After initialising the system as an amorphous fluid, the square well attraction is turned on and the template positioned, defining . An illustrative simulation snapshot shortly thereafter is shown in FIG. 2. In principle, the nucleation of independent crystals in the bulk fluid is possible, but we observe no such nucleation in the simulations presented here.

Figure 2: Illustrative snapshot of a monodisperse system with at , showing the crystal templating method. The FCC template is positioned at , and in this case 1 or 2 further crystalline layers have so far been deposited from the bulk fluid. Approximately half the length of the simulation cell is shown. Simulation snapshots were produced using OVITO Stukowski (2010).

Iii Theory

iii.1 Free energy landscape

Although this work is concerned with the dynamics of phase transitions, it is possible to gain a surprising level of insight into the expected behaviour from thermodynamic considerations alone Poon et al. (1999, 2000); Renth et al. (2001), by examining the free energy landscape of the system. To that end, we now calculate approximate fluid and crystal free energy curves for our system, in the monodisperse case, by perturbation to a hard-sphere reference system.

The free energy density is given by:


where and is the number density. The integral is over the square well range. For the hard-sphere contribution , we use the Carnahan-Starling free energy density Carnahan and Starling (1969) for the fluid, and Hall’s expression for the crystal Hall (1972). For the hard-sphere radial distribution function , we use the Percus-Yevick expression in the fluid Stell (1963) (using the analytical form for the first neighbour shell given in Ref. Trokhymchuk et al. (2005)) and that of Choi et al. for the crystal Choi et al. (1991).

Figure 3: Free energy density as a function of volume fraction for a monodisperse system with and . The solid curve shows the fluid branch of the free energy while the circles show the crystal branch. The dotted line indicates the common tangent for metastable gas-liquid coexistence, and the dashed line shows the tangent for the equilibrium crystal-gas coexistence. There is no common tangent between the liquid and crystal phases.

FIG. 3 shows the resulting free energy densities plotted as functions of volume fraction for the case (inside the G-L binodal, for our value of , at point 1 in FIG. 1). Allowed coexistences are given by the common tangent construction. Although quantitative agreement with the simulation data of Ref. Liu et al. (2005) (in terms of the coexistence volume fractions of the gas and liquid phases) is relatively poor, some important qualitative features are present. The common tangent linking the minimum of the crystal branch with the gas is lower than that linking the liquid and gas – the crystal-gas coexistence therefore has a lower overall free energy and is hence the equilibrium state, while the gas-liquid coexistence is metastable. Furthermore, there is no common tangent between the crystal and liquid, which means that on the basis of free energy considerations, the crystal cannot locally coexist with the metastable liquid.

The lack of a crystal-liquid common tangent implies a growth scenario like that proposed for the experimental colloid-polymer mixture observations in Refs. Evans et al. (2001); Renth et al. (2001), in which growing crystallites are coated by a layer of colloidal gas because of their inability to coexist with the metastable liquid. This ‘boiled-egg crystal’ should deplete the surrounding bulk fluid until the required (equilibrium) crystal-gas coexistence is achieved overall. Our simulation setup allows us to model the effects of a metastable gas layer on the growth of a single crystal in a system for which we have a calculation of the corresponding free energy landscape, allowing new insight into some of the physics described in Refs. Evans et al. (2001); Renth et al. (2001).

iii.2 Split interfaces

The expected growth scenario of the ‘boiled-egg crystal’ described in Section III.1 is that of a ‘split interface’. That is, the crystal should form an interface with a gaseous region, which in turn forms an interface with the metastable liquid, since the crystal and metastable liquid cannot coexist on the free energy grounds discussed above. Such interfaces and a method for predicting their growth rates are described in Refs. Evans and Poon (1997); Evans et al. (1997), wherein the example used was of a crystal-liquid-gas interface as opposed to the crystal-gas-liquid interface we expect here. We now make use of that theory to predict the evolution of the split interface that should be formed when our system is inside the G-L binodal.

As shown in FIG. 4, we assume the growth scenario to be such that the phases on each side of the interfaces are at their ‘correct’ densities according to the equilibrium phase diagram (FIG. 1), the densities corresponding respectively to: the metastable liquid, the metastable gas, the stable gas, and the stable crystal. The assumption is therefore one of local equilibrium at the interfaces.

Let be the position of the gas-liquid interface, and that of the crystal-gas interface. The respective fluxes on to and away from the interfaces are , , and , as shown in FIG. 4. We assume , i.e. that the liquid exists uniformly at its metastable density, and within the crystal region. The remaining fluxes can then be related to the interface speeds and the densities :


We now make the simplifying approximation that the gas region supports a uniform gradient between the densities and , as is shown in FIG. 4, and further that the diffusion constant therein is equal to the ideal Stokes-Einstein diffusion constant due to the low density of the gas. Therefore:


where . Substituting into Equations 3 and 4 and separating variables, with the necessary constants of integration given by the initial conditions , we find expressions for the interface positions:


where we have defined and .

To make contact with our simulations, we read off the values of from FIG. 1 at the appropriate effective temperature for point 1, inside the G-L binodal, and use the Stokes-Einstein diffusion constant defined by our choice of time unit, . The resulting interface positions through time are shown in FIG. 5, and in the next section are compared with simulation data. The characteristic diffusive growth is apparent, with the gas-liquid interface leading the crystal-gas interface as expected.

We note at this stage that the crystal-gas-liquid split interface scenario is closely analogous to the crystal-liquid-gas interfaces described in Refs. Evans and Poon (1997); Evans et al. (1997). The only difference is that in the present work, the ‘pivot’ phase which coexists locally with both others (and therefore coats the crystal) is the gas, i.e. is the equilibrium phase, whereas in the aforementioned work it is the metastable, i.e. nonequilibrium liquid phase.

Figure 4: Schematic density profile of the expected crystal-gas-liquid split interface inside the G-L binodal, where the crystal cannot coexist with the metastable liquid. The phases from left to right are the crystal, gas, and metastable liquid. The gradient in the gas region is drawn uniform, as is assumed in our calculations.
Figure 5: The predicted evolution of the crystal-gas (solid line) and gas-liquid (dashed line) interfaces for a split crystal-gas-liquid interface when our system is inside the G-L binodal (, ).

Iv Results and Discussion

In this section, we examine the simulation results in terms of their time-dependent concentration profiles and in terms of their crystal growth rates, elucidating and discussing the separate and combined effects of metastability and polydispersity on the growth process. Then, in Section V, we outline and briefly test a possible explanation for our findings in terms of fractionation and rearrangement at the crystal interface.

We present results for polydispersities of (monodisperse), and . The and coordinates used are marked as points ‘1’ and ‘2’ on the monodisperse phase diagram, FIG. 1. We use square well depths of (point 1, inside the G-L binodal, referred to as ‘gas-mediated’) or (point 2, outside the G-L binodal, referred to as ‘gas-free’). The combinations of and give 6 state points in total. At each state point, 6 independent initial configurations were used.

iv.1 Concentration profiles – monodisperse

The clearest way to visualise crystal growth over the whole simulation time is by plotting the time-dependent concentration profile. To achieve this, we plot the ‘area fraction’ of particles intersecting planes at a given coordinate, corresponding to a local volume fraction in an infinitely narrow interval . Since the simulation geometry is such that the crystal interface advances along the axis, these plots are a powerful way of observing the time evolution of concentration differences along the direction of growth while averaging over the axes perpendicular to it.

In this section, we qualitatively compare representative concentration profiles in the monodisperse and cases. Except where noted, the phenomena described occurred similarly in all trajectories of the state point in question.

Let us consider the monodisperse case first. FIGs. 6 and 7 show example trajectories for simulations at the gas-free (point 2 on FIG. 1) and gas-mediated (point 1 on FIG. 1) state points respectively. In the gas-free case, the growth scenario is, as expected, relatively simple. The bulk fluid remains homogeneous, since we are above the critical for G-L separation. The crystal template causes the growth of a single crystalline region with an average volume fraction of , which is the expected equilibrium volume fraction for the crystal Liu et al. (2005). The crystal, being at a higher volume fraction than the bulk fluid, depletes its surroundings of particles, resulting in a concentration gradient between the local fluid and the bulk; it is this gradient which transports particles toward the crystal so that it can continue to grow.

Figure 6: Time-dependent concentration profile along the axis for one of the monodisperse gas-free simulations. The greyscale indicates the local volume fraction, ranging from (white) to (black). As the simulation progresses, the crystal advances along the axis, depleting the fluid in front of the interface.
Figure 7: As FIG. 6, for a monodisperse gas-mediated simulation. Gas-liquid separation takes place in the bulk, and the advancing crystal is coated by a distinct gas layer which advances ahead of it. The theoretical crystal-gas and gas-liquid interface curves from FIG. 5 are rotated and superimposed, showing approximate agreement with the data.
Figure 8: Snapshot of the crystal interface for one particular monodisperse gas-mediated (i.e. inside the G-L binodal) trajectory in which the crystal happened to exhibit dendritic growth around an impinging gas bubble. Such dendritic growth was not generally observed.

Inside the G-L binodal (FIG. 7), the scenario is quite different. Firstly, there is clearly some G-L separation taking place in the bulk. Focusing next on the crystal template at , we can see that a region of low density gas forms in front of the crystal almost immediately, shielding it from the liquid with which it cannot locally coexist according to the free energy considerations in Section III.1. This is in contrast to FIG. 6, in which the fluid immediately next to the crystal template retains its density for quite some time until the growing crystal depletes it of particles. Note also that whereas the depleted region in FIG. 6 fades smoothly into the higher density bulk fluid, the gas next to the crystal in FIG. 7 forms a substantially sharper interface with the fluid next to it, indicating a distinct phase boundary between it and the metastable liquid. The formation of this well-defined, macroscopic gas layer shielding the crystal is consistent with the free energy considerations above, and with the experimental observations of Renth et al. (2001).

The speed of advance of the crystal-gas and gas-liquid interfaces is in approximate agreement (to within 1 particle diameter) with the theoretical prediction in Section III.2 (FIG. 5), which we have rotated and superimposed on FIG. 7 for comparison. This agreement is particularly satisfying given that the input densities required for the theoretical calculation in Section III.2 were taken from the equilibrium phase diagram, so that there are no free fitting parameters. The accuracy of the prediction thus gives credence to the assumptions of local equilibration of the interfaces and of essentially ideal gas-like Stokes-Einstein diffusion in the gas region. We note also that the theoretical prediction does not take account of the ‘start-up’ stage when the gas layer is just forming, since it assumes a gas layer (albeit one of zero size) to be present from – this may explain the slight discrepancy at early times.

Finally, it is interesting to note that for some trajectories, the crystal can grow incomplete layers, as shown in FIG. 8. This is the result of dendritic growth of the crystal around an impinging gas bubble. Such growth was not found to be generally present in the other trajectories at this state point – the crystal typically formed complete layers covered by a layer of gas spanning the whole width of the simulation cell.

iv.2 Concentration profiles – polydisperse

Figure 9: As FIG. 6, for a gas-free system of polydispersity . As in the monodisperse case, the bulk fluid remains homogeneous. However, the crystal template does not induce any growth on the timescale simulated.
Figure 10: As FIG. 6, for a gas-mediated system of polydispersity . The bulk fluid shows G-L separation, and the gas layer is present in front of the crystal. In contrast to FIG. 9, the crystal is able to grow, although more slowly than in the corresponding monodisperse case (FIG. 7).

We now discuss the polydisperse case. FIG. 9 shows the gas-free scenario, corresponding again to point 2 on FIG. 1. On the simulated timescale, essentially no crystal growth is seen to take place, except perhaps for a slight ‘wetting’ of the template by ordered particles. As in the monodisperse case, the bulk fluid remains homogeneous, indicating that the polydispersity has not altered the location of the critical in such a way as might bring this state point inside the G-L binodal.

Moving to the gas-mediated case inside the binodal (at point 1 on FIG. 1), for which a trajectory is shown in FIG. 10, we can see that the crystal is able to grow substantially despite the polydispersity. As in the corresponding monodisperse case (FIG. 7), the bulk fluid separates into gas and liquid regions, and a gas layer covers the crystal as it grows. An illustrative snapshot of the gas-coated polydisperse crystal is shown in FIG. 11.

To summarise the qualitative observations in this section: the crystal template successfully induces a physically realistic growth process, in which the crystal depletes its surroundings and forms ordered, high density layers. When the system is inside the G-L binodal (at point 1 in FIG. 1), a gas layer forms in front of the crystal, shielding the crystal from the metastable liquid. This is exactly the scenario predicted for our system in Section III.1 and proposed in Refs. Evans et al. (2001); Renth et al. (2001). The presence of polydispersity substantially slowed the crystal growth overall. In the gas-free case (point 2 on FIG. 1) the crystal showed barely any growth on the simulated timescale, whereas the gas-mediated case (point 1) showed substantial crystal growth despite the polydispersity.

Figure 11: Illustrative snapshot of a polydisperse gas-mediated trajectory at , showing the gas-coated crystal region. Particles flagged crystalline by the bond order parameter algorithm described in Section II.2 are shown in black, all others are light grey.

iv.3 Monodisperse crystal growth rates – the effect of metastability

Having summarised the qualitative features of the monodisperse and simulations, we next quantitatively examine the effect of those features on the crystal growth rate. This is done by plotting the crystallinity (the proportion of particles flagged as crystalline, according to the bond order criterion described in Section II.2) through time. FIG. 12 shows the crystal growth at each of the 6 state points studied, each averaged over the 6 independent realisations of the state point.

In the monodisperse case, the effect of metastability alone is apparent. The gas layer in front of the crystal, which forms when we are at point 1 in FIG. 1, significantly slows crystal growth compared to the gas-free case, corresponding to point 2 in FIG. 1.

Since the crystal exists at a higher density than the rest of the system, as soon as the immediate surroundings are depleted of particles, particles must be transported toward the crystal from the bulk in order for it to grow. This transport takes place via collective diffusion down a concentration gradient between the low concentration near the crystal and the relatively higher concentration away from the crystal. The diffusion is described by Fick’s law,


in which is the concentration flux, is the local concentration and is the (concentration-dependent) collective diffusion coefficient. Note that the collective diffusion coefficient, in contrast to the self diffusion coefficient, increases with Tirado-Miranda et al. (2003).

We now consider the phase diagram for our system, as shown in FIG. 1 (c). In the absence of the gas layer, and assuming that the interface is locally equilibrated, a concentration gradient exists between at the interface and in the bulk. However, inside the G-L binodal, when the crystal is covered by a distinct gas layer, the relevant concentration gradient is that which exists across the gas layer itself – the gas-liquid interface beyond truncates the concentration gradient. The gradient in the gas layer is between the equilibrium gas () and the slightly higher concentration metastable gas at . The gas layer therefore results in a significantly smaller concentration difference being spread across a comparable distance – from FIGs. 6 and 7, the gas layer is similar in size to the depletion region in the gas-free case.

Therefore, the low-concentration gas layer affects particle transport in two ways, relative to the gas-free case: (a) the collective diffusion coefficient is lower, due to reduced concentration; (b) the gas’s low concentration and macroscopic size mean it can only support a relatively weak concentration gradient in comparison to the gas-free case. As evidenced in the monodisperse crystal growth rates (FIG. 12), the net result is that the flux onto the crystal, and therefore the interface speed, is lower when the gas is present. Although, on thermodynamic grounds, the driving force for crystallisation appears higher inside the G-L binodal, the growth of a given crystal is slowed due to the gas layer’s effects on particle transport to the interface.

The effect of the gas layer in the monodisperse case is comparable to that of the liquid layer in Refs. Evans and Poon (1997); Evans et al. (1997), which was also found to inhibit crystal growth. There, as here, the crystal is coated by a phase which advances ahead of the crystal and slows down its growth. It is interesting that this comparison holds, given that in our case the crystal is being coated by the equilibrium gas, whereas Refs. Evans and Poon (1997); Evans et al. (1997) concern a crystal coated by a nonequilibrium liquid.

Figure 12: (colour online) Crystallinity through time for the 6 state points studied, in each case averaged over 6 independent realisations. Filled symbols: gas-mediated (point 1 on FIG. 1); open symbols: gas-free (point 2 on FIG. 1). Triangles: monodisperse; circles: ; squares: . The standard errors are approximately the size of the symbols. The filled triangles (monodisperse gas-mediated) are shown in grey (red online) to distinguish them from the filled circles ( gas-mediated).

iv.4 Polydisperse crystal growth rates

We consider first the highest polydispersity studied here, . As shown in FIG. 12 and in agreement with the observations in Section IV.2, the gas-free state point now shows only very slight growth on the timescale simulated, whereas the gas-mediated state point shows significant crystal growth (although less than in either of the monodisperse systems). This is a qualitative difference compared to the monodisperse case, in which the gas layer instead strongly slows the crystal growth. That is, the presence of polydispersity qualitatively alters the effect of the metastable gas layer on crystal growth.

In the case, the gas-mediated system is initially slightly faster, before being overtaken by the gas-free system by around . This is suggestive of some kind of crossover phenomenon between whatever factors are dominant in the fully monodisperse and cases.

Given the relative proximity of state point 1 to the gas-liquid critical point (FIG. 1), the possible influence of critical fluctuations was considered. We therefore tested a stronger well parameter, , , at , to see any effect of moving further away from the gas-liquid critical point. The crystal growth rate was unaffected within error, suggesting that the critical point is not close enough to state point 1 for critical fluctuations to have an effect. In addition, the observed gas-liquid interface on the far side of the gas layer appears sharp – we do not observe the large fluctuations characteristic of critical phenomena.

Also for , moving the (gas-free) state point 2 to (to increase the driving force for crystallisation) slightly speeded up growth, but it was still much slower than in the gas-mediated case. Further study in this direction is left to future work.

V Crystallisation mechanism

v.1 Fractionation

The results presented so far show that metastable gas-liquid separation can affect the rate of crystal growth, because it results in a gaseous layer coating the crystal as it grows. However, the qualitative nature of that effect depends strongly upon polydispersity. In the monodisperse case, the resultant gas layer impedes crystal growth by reducing the efficiency of particle transport to the crystal. At , the metastable separation instead strongly enhances crystal growth. It is this effect of polydispersity on the crystal growth mechanism that we now discuss.

We propose an explanation in terms of a local fractionation process at the interface of the polydisperse crystal. In size polydisperse systems, it is well known that phase separation is typically associated with some degree of thermodynamically-driven fractionation in which one phase e.g. contains on average larger or smaller particles than another, changes its overall polydispersity, etc. Such fractionation has been predicted theoretically Evans (2001); Fasolo and Sollich (2004), and observed experimentally Evans and Fairhurst (2004) and in simulation Williamson and Evans (2012). Of course, this kind of fractionation pertains to the equilibrium phase composition and may not be fully achieved in real systems on accessible timescales due to kinetic factors. This is a particularly important consideration for the crystal phase, where particles are essentially stationary once incorporated and, thereafter, the long-range particle transport required for fractionation is facilitated only by the presence of defects. In general, phase composition is expected to relax slowly in comparison with overall density so that fractionation may lag behind phase separation somewhat B. Warren (1999), leading in the crystal phase to the ‘freezing-in’ of a nonequilibrium composition Evans and Holmes (2001). (Note, however, that recent work by us has shown that some fractionation is possible in the very early stages of gas-liquid separation Williamson and Evans (2012).)

Nevertheless, we would expect that fractionation at the interface of the growing crystal, where particles are still mobile and may be easily exchanged with the fluid, is quite feasible. Indeed, the slowing effect of polydispersity on crystallisation is generally taken as evidence that fractionation is involved to some extent Martin et al. (2003); Schöpe et al. (2007). We propose that such interfacial fractionation takes place not just in nucleation but during the subsequent growth of the crystal, so that it is facilitated by the low-density gas layer that forms in front of the crystal when inside the G-L binodal. In contrast to the collective diffusion discussed in Section IV.3, the self diffusion required for fractionation is hindered by particle density, so that fractionation would be frustrated outside the G-L binodal, where the fluid side of the interface is at a relatively high density compared to that in the gas-mediated case.

In addition to an overall preference for a narrower size distribution, it seems reasonable to suppose that the crystal may be selective in some way as to which particles are incorporated where, on a local basis. For instance, it may be frustrated by the presence of regions occupied solely by very large or small particles. Whatever the nature of any local size ordering, we expect that the crystal should have some kind of preference as to how particle size is distributed within it, given that the largest and smallest particles in the system differ in size by around of the mean. This local ordering is conceived of as being in addition to any fractionation in terms of the overall balance of particle sizes between the phases.

From the outset, we stress that a full explanation of the results in the previous section requires further work to elucidate and test other possible contributing factors, such as the influence of polydispersity on the equilibrium phase diagram. The proposals here are motivated by the dynamical nature of our simulation approach and by the intuition that the increased diffusivity of particles in the gaseous layer would help any dynamical sorting processes at the crystal interface, a consideration which we expect to be a significant aspect of any full explanation of the results. In the following, we perform further analysis on the simulation data in an attempt to detect the presence of fractionation processes.

v.2 Polydispersity of the crystal

The results we now present concern the gas-mediated (point 1 in FIG. 1) and gas-free (point 2 in FIG. 1) simulations. To ensure good sampling in the crystal, we allowed these to run up to . By this time, the crystal interface in the gas-mediated simulations reached , while the gas-free simulations had only grown to . We check for fractionation by measuring the mean diameter and polydispersity in the crystal. For the purposes of the present analysis, we identify the crystal interface at a given time as the centre of an averaging window of length in in which, moving through the system from , the average local volume fraction first drops below the parent value , due to the depleted region in front of the crystal. The identified interface positions are checked against the concentration profiles (e.g. FIG. 7) and against direct visual observations, and the region lower in than the interface is considered to be the crystalline region. For the case of FIG. 8, in which the crystal interface is not flat, our definition means that only full or nearly-full crystal layers fall into the ‘crystalline region’.

To within statistical error, was equal to the parental mean, . However, changes in were detected. In FIG. 13, the evolution of in the gas-mediated and gas-free crystals is shown, as a function of time and of , the mean number of particles in the crystal. In the gas-mediated case, there is a clear reduction in the polydispersity of the crystal compared to the parental , demonstrating that the crystal is selecting a narrower subset of the parent size distribution, eventually attaining a value of at the end of the simulated time. This is qualitatively consistent with Fasolo and Sollich’s equilibrium calculations on hard-spheres, in which even a small parental polydispersity was found to lead to reduced polydispersity in a crystal coexisting with a fluid Fasolo and Sollich (2004). This measurement in itself is significant: to our knowledge, such fractionation has not previously been measured in the dynamical crystal growth of polydisperse spherical colloids (an analogous example for colloidal platelets exists in Ref. Byelov et al. (2010)), although it is proposed as an explanation of slow crystallisation in such systems Martin et al. (2003); Schöpe et al. (2007).

The data for the gas-free crystal are subject to significantly larger error, because very little crystal has actually formed, but appear to be consistent with a small reduction in polydispersity corresponding to the frustrated crystal growth. In any case, it is clear that the successful growth of a crystal involves, as expected from equilibrium work Fasolo and Sollich (2004); Sollich and Wilding (2011), a measurable reduction in polydispersity from the bulk fluid, indicating a fractionation process which would be enhanced by the gas layer.

Figure 13: (colour online) Crystal polydispersity as a function of (upper pane) and (lower pane) for the gas-mediated (black) and gas-free (grey, blue online, larger error bars) crystals. The error on the axis is approximately 30.
Figure 14: (colour online) Upper pane: the local size correlation function as defined in the text, measured at late times for systems of polydispersity . Black crosses indicate the gas-mediated crystal, red circles the gas-mediated outside-crystal region (comprising gas and liquid coexisting), and blue squares the gas-free outside-crystal region (a homogeneous fluid). Lower pane: the corresponding radial distribution functions measured in the same systems. Colours and symbols as for the upper pane. Both the outside-crystal datasets have been multiplied by a factor of 2 for ease of visualisation. Lines between the data points are guides to the eye. Vertical lines linking the panes are provided to show the relationship between the functions – dashed lines for the crystal, and dotted lines for the two fluid datasets, whose dependence on is essentially identical. The data points are placed in the centre of each bin in the axis.

v.3 Local size correlations

To check for local size ordering, in addition to the overall preference for a narrower particle distribution, we define a local diameter-diameter correlation function :


where the averaging in the first term is over neighbours separated by a distance . The function is therefore analogous to a ‘radial distribution function,’ but for size correlations, as opposed to density correlations, and tends to zero in the limit of an ideal system (i.e. one in which particles do not interact and ‘size’ just becomes a label which does not affect the behaviour of the particles).

Measurements are taken from the simulations by binning in , and can be made within either the crystalline or amorphous region, to observe any differences induced by crystallisation. Then, the averaging over the second term in Equation 9 is over all particles in that phase. We compare in the gas-mediated crystal with that outside the crystal in the same system (remember that in the gas-mediated case, ‘outside the crystal’ refers to a metastable coexistence of gas and liquid in the bulk). We also measure outside the crystal in the gas-free case, i.e. a homogeneous fluid of .

The results of this analysis are shown in FIG. 14, in which are also shown the measured radial distribution functions corresponding to the data. The data have been averaged within each simulation over , during which time both and appeared essentially static due to the very slow growth of the crystal (whose size scales as ), and then averaged over the independent simulations.

The comparison reveals a fascinating ’-like’ appearance to the local size correlations both inside and outside the crystal. However, closer inspection reveals that whereas the oscillations in for the crystal appear approximately in phase with those in its , those in the fluid – both the homogeneous and G-L separating fluid cases – are approximately in phase quadrature with .

Our proposed explanation for these findings is as follows. We will consider first the gas-mediated crystal, then the two fluid datasets. Before starting, we note that, in a polydisperse system, larger particles will tend to contribute a stronger signal to structural correlation functions. See, for instance, the measurements of partial static structure factors in Ref. Weysser et al. (2010). One may therefore be guided by this in interpreting : e.g. a positive value indicates a correlation between either ‘big-big’ or ‘small-small’ particle pairs, but will contain a stronger contribution from ‘big-big’ pairs due to the stronger overall structuring of large particles. The terms ‘big’ and ‘small’ are here shorthand for particles greater or lesser in size than the mean .


In the crystal, the features in more or less mimic those of , and the two are approximately in phase, aside from a slight offset for the first peak which decreases through the subsequent peaks. Where there is a peak in , i.e. a pair of neighbours is likely to be found separated by this distance, those particles on average show a slight positive correlation in terms of their size, indicating that the crystal structure has a preference for particles in the same region to be of similar size, presumably to minimise distortion of the lattice. Conversely, particle pairs separated by the unusual distances corresponding to minima in have a rather strong negative correlation. This is consistent with a scenario in which these neighbour separations are associated with the presence of ‘wrong’, i.e. unusually-sized, particles, for a given region of the crystal. They may, for instance, result from an interstitial defect, with a particularly small particle squeezing itself in amongst a region of average or slightly large on-lattice particles.

The relationship between and in the crystal is therefore due to the need to minimise lattice distortions over small regions of the crystal, and the fact that unusual neighbour separations will tend to be associated with particles of unusual size for a given region distorting the crystal structure. We stress that this effect is in addition to the overall preference for a reduced crystal polydispersity. Conceivably, the crystal could have reduced its polydispersity precisely in order to avoid having to enact the local ordering observed. In fact, for the dynamically grown crystal here, a little of both effects seems to be required: the crystal reduces its polydispersity, but still ‘cares’ about the local distribution of particle size in its lattice. Since the crystals grown in our dynamical simulations are not necessarily at equilibrium, it would be interesting to compare these results to those from an equilibrated polydisperse crystal. We intend to pursue this question in future work.


Both the homogeneous (gas-free, point 2 in FIG. 1) and G-L separating (gas-mediated, point 1 in FIG. 1) fluid regions show similar behaviour with respect to the -dependence of and , so we will not distinguish between them for the purposes of this analysis. In the fluid, appears to be approximately in phase quadrature with , so that minima or maxima in appear, respectively, halfway up or halfway down the slope around a peak in .

We interpret this in terms of shells around a nominal test particle , in analogy to the explanation of similar oscillations in terms of shells near a hard wall in Refs. Pagonabarraga et al. (2000); Buzzacchi et al. (2004). Consider that a single shell of surrounding particles corresponding to a peak in may be roughly divided into those which are nearer than average and those which are further away than average. In order to most efficiently fill space, particles which are closer should be small, and those which are further away should be large. Put another way, particles which are further away should be further away because they are big; placing small particles further away than is required by their size would waste space, as it were, reducing the packing efficiency and increasing the free energy of the fluid. Given the stronger structural signal from larger particles Weysser et al. (2010), the data are dominated by the case where particle is big. Therefore, the smaller, closer, particles in the shell will on average contribute negatively to , resulting in a minimum therein at a value of slightly less than the peak in . Conversely, the larger, further away particles contribute positively, giving a peak in just outside the peak in .

In the DFT study of Ref. Pagonabarraga et al. (2000), an observed phase quadrature between the size distribution near a wall and the mean density profile (roughly corresponding to our and functions) was explained in these terms. Our findings are therefore significant in showing that similar ‘local size segregation’ appears to be present in the bulk fluid, rather than being solely the result of spatial inhomogeneity. We note for completeness that measuring in a cubic system with no crystal template (so that our system is truly homogeneous and isotropic) had no significant effect on the behaviour of ; our measurements do seem to capture the behaviour of a bulk fluid.

Thus, the observed relationship between and in the fluid is a manifestation of the need for efficient packing in the fluid, coupled with the tendency for larger particles to contribute more strongly to the structuring signal. The fact that local size segregation previously observed near a hard wall appears to be manifest also in a bulk fluid is very interesting, and would seem to merit further investigation. In particular, we aim to investigate in future work whether this same effect is present in an equilibrium polydisperse fluid.

v.4 Equilibrium insights

Although our simulations are inherently not at equilibrium, it is important to consider the effects of polydispersity on the equilibrium phase diagram in interpreting the results. The most detailed investigations into polydisperse equilibria are the theoretical modelling of Sollich and collaborators, based on the moment free energy method Fasolo and Sollich (2004); Sollich and Wilding (2011). A polydisperse equilibrium phase diagram is not presently available for the square well interaction considered here. Still, the hard sphere calculations in Ref. Fasolo and Sollich (2004) yield important qualitative insights. Within the hard sphere crystal-fluid coexistence region, even small parent polydispersity is associated with reduced polydispersity in the crystal, as observed in our dynamical simulations. For higher parent polydispersity and volume fraction, there is the possibility at equilibrium of multiple crystalline phases. However, the qualitative picture within the crystal-fluid coexistence region (so that the crystal can reject particles into the coexisting fluid), at the polydispersities we consider here, is of a fluid coexisting with a single crystal phase, of reduced polydispersity relative to the parent Fasolo and Sollich (2004). Additionally, when such multiple crystal phases do enter into the phase diagram, they may not be kinetically attainable on experimental timescales Liddle et al. (2011) – hence the question of how a single crystalline phase grows would remain important.

We note also that, within the moment free energy method, the present best assumption is that the equilibrium crystal(s) have substitutionally-disordered FCC structure, as is assumed for our crystal template. Against the backdrop of overall substitutional disorder (i.e. an FCC-like structure, as opposed to an alloy), we found (Section V.3) that some local ordering is detectable, i.e. the crystal ‘notices’ particle size in deciding the local distribution of particles. For high polydispersity, it has been noted Sollich and Wilding (2011) that substitutionally-ordered alloy structures such as may become preferable. No calculations have yet dealt with this question, however. Therefore, for the present degree of polydispersity, with the best current equilibrium information, our templating strategy seems realistic, in that it allows the expected equilibrium crystal to form.

v.5 Summary

The results above show that the crystal phase is able to reduce its polydispersity relative to the fluid even as it is growing. In addition, we have measured a local radial size correlation function to quantify the nature of local size ordering in each phase. Both the overall form of and its relationship to must change qualitatively between the fluid and crystal phases. This indicates a further, local ordering process in addition to the overall preference of the crystal to reduce its polydispersity. To our knowledge, this phenomenon has not previously been observed, and gives new insight into the detailed structural features of polydisperse phases. We have also found that local size segregation, previously found in a polydisperse fluid near a wall, seems to have a close analogue in the structuring of a bulk polydisperse fluid.

We have proposed that the necessary fractionation and rearrangement processes must take place at the crystal interface and therefore, relying on self diffusion near the interface, are enhanced by the presence of the gas layer. If true, this mechanism would explain the very slow crystallisation in the gas-free case, where no such layer exists.

The simulations provide an interesting intermediate case, in which the crystal growth in the gas-mediated case is initially faster before being overtaken by the gas-free case. In light of our hypothesis, this may be because the gas-free growth is relatively slower until the interfacial fluid is sufficiently depleted to allow adequate diffusion near the interface. In the gas-mediated growth, a gas layer is formed immediately, so that growth is initially faster. However, once the gas-free system has sufficiently depleted the interfacial fluid, the interfacial diffusion is enhanced and it now grows the fastest, the gas-mediated system lagging behind for the reasons discussed in Section IV.3. In contrast, for the gas-free growth is so slow (presumably due to the greater fractionation required) that the interfacial fluid is not depleted enough for such a crossover to happen on the simulated timescale.

We stress that the ideas outlined in this section constitute only one potential explanation of the results we observed in Section IV. Further work will be required to determine other factors (e.g. a quantitative calculation of polydispersity-induced changes to the equilibrium phase diagram of square well particles), and to more rigorously test the dynamical explanation outlined here, which we expect will remain an important part of any later explanation. Nonetheless it is clear that, whatever the mechanisms, the influence of metastability on crystal growth kinetics is far from trivial, and that its role can be qualitatively switched by the presence of a very mild degree of polydispersity.

Vi Conclusions

We have studied the kinetics of crystal growth using a model system in which a single crystal grows from an amorphous fluid. By varying the interaction strength and the polydispersity, we were able to investigate the separate and combined effects of metastable gas-liquid separation and polydispersity on the crystal growth process. Taking advantage of the pre-determined growth direction of the crystal, we have used one-dimensional concentration profiles to clearly display the phase ordering scenarios when the system is inside or outside the metastable G-L binodal. Inside the binodal, we observed the formation of a gaseous layer coating the crystal, in agreement with approximate free energy curves and with previous experimental observations of three-phase ordering. The diffusive growth of the crystal-gas-liquid ‘split interface’ was modelled theoretically, and we found approximate agreement between theoretical and simulation interface growth profiles. This agreement validated our use of two approximations in the theory, namely that the interfaces are locally equilibrated to the densities given by the phase diagram (FIG. 1) and that particle transport through the gas layer is approximately ideal gas-like due to the gas’s low density.

The dynamic influence of this gaseous layer, in the monodisperse case, was to impede growth by slowing particle transport to the crystal interface, suggesting that, while metastable G-L separation can enhance crystal nucleation Anderson and Lekkerkerker (2002); Fortini et al. (2008), it may (in the ideally monodisperse case at least) have the opposite effect on the subsequent growth of the crystal if the free energy landscape is such that the metastable liquid cannot coexist locally with the crystal, as was the case here. However, introducing a small amount of polydispersity, corresponding roughly to that found in a ‘near-monodisperse’ colloidal system, we found that the gaseous layer instead strongly enhanced growth.

The crystals grown show an overall reduction in polydispersity versus the parent, consistent with existing equilibrium theory and equilibrium simulations. This fractionation is often proposed as an explanation for slow crystallisation in polydisperse systems, but to our knowledge has not previously been measured in the dynamical crystal growth of polydisperse spherical colloids. We have postulated that this fractionation process, requiring self diffusion, is facilitated by the gas layer, perhaps explaining why this layer enhances crystal growth in the polydisperse case, rather than impeding it as in the monodisperse case.

Additionally, we postulated and observed in detail a process of local size ordering in crystal formation. We defined and measured a radial size correlation function , whose appearance resembles the radial distribution function in each phase. We have described and explained the relationship between these two functions in both the crystal and fluid phases, showing that both the form of and its relationship to must be qualitatively altered as the crystal is grown. The findings for the fluid phase seem to indicate that local size segregation previously observed near a hard wall also appears in the structuring of a bulk fluid.

These measurements indicate a local ordering process, in addition to the overall preference of the crystal to reduce its polydispersity relative to the fluid, which would also be enhanced by the presence of the gas layer. Aside from this, further study of the detailed form of , the dependence on system parameters, equilibration, etc. could provide a useful new angle on the influence of polydispersity on local crystal structure. For instance, preliminary work in this direction (not shown here) shows that such correlations persist even when the parent polydispersity is so low that a dynamically-grown crystal does not reduce its polydispersity at all relative to the parent, i.e. overall fractionation does not take place.

It must be stressed that our proposed explanation requires further investigation and probably does not constitute the full story. Other factors, discussed in Sections V.4 and V.5, must be considered – further work will be required to establish their influence, if any, and to rigorously test the explanation we outlined here. We used two state points (in and ) in order to focus on the effects of multiple polydispersities at these points – it would be beneficial in future work to explore the phase diagram further, for instance by ‘exiting’ the gas-liquid binodal in another direction (e.g. increasing/decreasing as well as varying ). We have, though, mentioned simulations run at lower , further from the gas-liquid critical point, which suggest that the effect of the gas layer here is not dependent on critical phenomena. Also, as with other comparable simulations, hydrodynamic interactions (HI) have been neglected; more complex and expensive simulations would be required to quantify possible effects of HI on the dynamics observed here.

We have presented a detailed dynamical study of crystal growth scenarios in the presence of metastability and polydispersity, and shown how our findings relate to the free energy landscape of the system. Whatever emerges from future work, the results demonstrate the importance of both metastability and polydispersity for soft matter phase transition kinetics and, moreover, that these two factors interact in a complex and previously unknown manner, the causes of which remain to be investigated further.

Thanks are due to P. Bartlett and C. P. Royall for their interest and suggestions, and in particular to N. B. Wilding for a useful discussion regarding local size ordering. J. J. W. was supported by an Engineering and Physical Sciences Research Council Doctoral Training Grants award.


  1. The Bates distribution is a sum of random numbers uniformly distributed on the unit interval – in our implementation .


  1. V. Prasad, D. Semwogerere,  and E. R. Weeks, J. Phys. Condens. Matter 19, 113102 (25pp) (2007).
  2. L. Onsager, Ann. N.Y. Acad. Sci. 51, 627 (1949).
  3. R. M. L. Evans, D. J. Fairhurst,  and W. C. K. Poon, Phys. Rev. Lett. 81, 1326 (1998).
  4. R. M. L. Evans, J. Chem. Phys. 114, 1915 (2001).
  5. P. Sollich and N. B. Wilding, Phys. Rev. Lett. 104, 118302 (2010).
  6. P. Sollich and N. B. Wilding, Soft Matter 7, 4472 (2011).
  7. M. Fasolo and P. Sollich, Phys. Rev. E 70, 041410 (2004).
  8. N. B. Wilding and P. Sollich, Europhys. Lett. 67, 219 (2004).
  9. R. Fantoni, D. Gazzillo, A. Giacometti,  and P. Sollich, J. Chem. Phys. 125, 164504 (2006).
  10. P. Bartlett and P. B. Warren, Phys. Rev. Lett. 82, 1979 (1999).
  11. E. Zaccarelli, C. Valeriani, E. Sanz, W. C. K. Poon, M. E. Cates,  and P. N. Pusey, Phys. Rev. Lett. 103, 135704 (2009).
  12. S. R. Williams, C. P. Royall,  and G. Bryant, Phys. Rev. Lett. 100, 225502 (2008).
  13. S. Martin, G. Bryant,  and W. van Megen, Phys. Rev. E 67, 061405 (2003).
  14. S. Martin, G. Bryant,  and W. van Megen, Phys. Rev. E 71, 021404 (2005).
  15. H. J. Schöpe, G. Bryant,  and W. van Megen, J. Chem. Phys. 127, 084505 (2007).
  16. S. M. Liddle, T. Narayanan,  and W. C. K. Poon, J. Phys. Condens. Matter 23, 194116 (2011).
  17. P. B. Warren, Phys. Chem. Chem. Phys. 1, 2197 (1999).
  18. R. M. L. Evans and C. B. Holmes, Phys. Rev. E 64, 011404 (2001).
  19. V. J. Anderson and H. N. W. Lekkerkerker, Nature 416, 811 (2002).
  20. N. Dixit and C. Zukoski, J. Colloid Interface Sci 228, 359 (2000).
  21. C. Haas and J. Drenth, J. Phys. Chem. B 104, 368 (2000).
  22. C. Haas and J. Drenth, J. Cryst. Growth 196, 388 (1999).
  23. J. J. Williamson and R. M. L. Evans, Phys. Rev. E 86, 011405 (2012).
  24. The Bates distribution is a sum of random numbers uniformly distributed on the unit interval – in our implementation .
  25. H. Liu, S. Garde,  and S. Kumar, J. Chem. Phys. 123, 174505 (2005).
  26. P. R. ten Wolde, M. J. Ruiz-Montero,  and D. Frenkel, Faraday Discuss. 104, 93 (1996).
  27. A. Stukowski, Modelling Simul. Mater. Sci. Eng. 18, 015012 (2010).
  28. W. C. K. Poon, F. Renth, R. M. L. Evans, D. J. Fairhurst, M. E. Cates,  and P. N. Pusey, Phys. Rev. Lett. 83, 1239 (1999).
  29. W. C. K. Poon, F. Renth,  and R. M. L. Evans, J. Phys.: Condens. Matter 12, A269 (2000).
  30. F. Renth, W. C. K. Poon,  and R. M. L. Evans, Phys. Rev. E 64, 031402 (2001).
  31. N. F. Carnahan and K. E. Starling, J. Chem. Phys. 51, 635 (1969).
  32. K. R. Hall, J. Chem. Phys. 57, 2252 (1972).
  33. G. Stell, Physica 29, 517 (1963).
  34. A. Trokhymchuk, I. Nezbeda, J. Jirsak,  and D. Henderson, J. Chem. Phys. 123, 024501 (2005).
  35. Y. Choi, T. Ree,  and F. H. Ree, J. Chem. Phys. 95, 7548 (1991).
  36. R. M. L. Evans, W. C. K. Poon,  and F. Renth, Phys. Rev. E 64, 031403 (2001).
  37. R. M. L. Evans and W. C. K. Poon, Phys. Rev. E 56, 5748 (1997).
  38. R. M. L. Evans, W. C. K. Poon,  and M. E. Cates, Europhys. Lett. 38, 595 (1997).
  39. M. Tirado-Miranda, C. Haro-PŽrez, M. Quesada-PŽrez, J. Callejas-Fern‡ndez,  and R. Hidalgo-çlvarez, J. Colloid Interface Sci. 263, 74 (2003).
  40. R. M. L. Evans and D. J. Fairhurst, Colloid Polym. Sci. 282, 766 (2004).
  41. D. V. Byelov, M. C. D. Mourad, I. Snigireva, A. Snigirev, A. V. Petukhov,  and H. N. W. Lekkerkerker, Langmuir 26, 6898 (2010).
  42. F. Weysser, A. M. Puertas, M. Fuchs,  and T. Voigtmann, Phys. Rev. E 82, 011504 (2010).
  43. I. Pagonabarraga, M. E. Cates,  and G. J. Ackland, Phys. Rev. Lett. 84, 911 (2000).
  44. M. Buzzacchi, I. Pagonabarraga,  and N. B. Wilding, J. Chem. Phys. 121, 11362 (2004).
  45. A. Fortini, E. Sanz,  and M. Dijkstra, Phys. Rev. E 78, 041402 (2008).
This is a comment super asjknd jkasnjk adsnkj
The feedback cannot be empty
Comments 0
The feedback cannot be empty
Add comment

You’re adding your first comment!
How to quickly get a good reply:
  • Offer a constructive comment on the author work.
  • Add helpful links to code implementation or project page.