# Kinetically driven ordered phase formation in binary colloidal crystals

## Abstract

The aggregation of binary colloids of same size and balanced charges is studied by Brownian dynamics simulations for dilute suspensions. It is shown that, under appropriate conditions, the formation of colloidal crystals is dominated by kinetic effects leading to the growth of well-ordered crystallites of the sodium-chloride (NaCl) bulk phase. These crystallites form with very high probability even when the cesium-chloride (CsCl) phase is more stable thermodynamically. Global optimization searches show that this result is not related to the most favorable structures of small clusters, that are either amorphous or of CsCl structure. The formation of the NaCl phase is related to the specific kinetics of the crystallization process, which takes place by a two-step mechanism. In this mechanism, dense fluid aggregates form at first and then crystallization follows. It is shown that the type of short-range order in these dense fluid aggregates determines which phase is finally formed in the crystallites. The role of hydrodynamic effects in the aggregation process is analyzed by Stochastic Rotation Dynamics - Molecular Dynamics simulations, finding that these effects do not play a major role in the formation of the crystallites.

## I Introduction

Colloidal suspensions have been widely studied for a variety of applications. From the point of view of basic science, a growing interest stems from the fact that colloidal suspensions are many-particle systems that are easier to observe than molecular systems. Aggregation and nucleation phenomena of colloids are therefore simpler to analyze, also by means of real-space techniques. Moreover, the interactions between colloids can be tuned in a controlled way so that a variety of structures with different properties can be experimentally produced Leunissen et al. (2005). Recently, much attention has been devoted to binary colloids in which two types of colloidal particles (types A and B in the following) can acquire opposite charges in suspension, interacting mainly through screened electrostatic forces. Binary colloids can thus aggregate López-López et al. (2006); Rollie and Sundmacher (2008); Leb (), in a process which is usually called heteroaggregation. The competition between attractive and repulsive forces gives rise to a rich phenomenology, which is the subject of active investigation. In fact, binary colloids have been shown to be able to form either compact colloidal crystals with a good degree of long-range order Leunissen et al. (2005); Bartlett and Campbell (2005); Hynninen et al. (2006), or colloidal gels characterized by a network of fractal-like aggregates Sanz et al. (2008); Piechowiak et al. (2012); Russell et al. (2012), the latter being more tenuous than in single-component gels Russell et al. (2012).

From the theoretical/computational point of view, notable attention has been devoted both to the determination of the equilibrium properties of binary suspensions López-López et al. (2006); Sanz et al. (2008); Bier et al. (2010); Pavaskar et al. (2012), with the calculation of phase diagrams of representative model systems, and to their aggregation kinetics López-López et al. (2006); Sanz et al. (2007, 2008); Cerbelaud et al. (2010a); Piechowiak et al. (2010). However, several aspects of heteroaggregation are still to be understood, especially for what concerns the basic mechanisms that determine the structures obtained in the actual aggregation processes.

An important issue in this colloidal aggregation is to determine to what extent the structures that are actually formed in the growth of colloidal crystals reflect their thermodynamic equilibrium structures. This has been investigated in the case of the heteroaggregation of colloids interacting through screened electrostatic forces by Sanz et al. Sanz et al. (2007), in the case of high volume fraction . In their simulations, Sanz et al. used the forward-flux sampling method Allen et al. (2005) to determine nucleation rates. They found that highest nucleation rate was associated to the formation of a substitutionally disordered fcc crystal, which is neither the lowest free-energy phase nor the phase characterized by the lowest free-energy barrier for nucleation from the melt. In fact, the lowest free-energy structure is of cesium-chloride (CsCl) type. Sanz et al. attributed this finding to the fact that sub-critical small clusters have disordered fcc structures, so that growth continues keeping this kind of structure because the transition to the CsCl structure is kinetically inhibited.

In this article we perform Brownian dynamics simulations of crystal nucleation and growth in the same model system as Sanz et. al. Sanz et al. (2007, 2008), consisting of binary spherical colloids of same radius and absolute value of the charge, interacting through a screened electrostatic potential of Yukawa-like form. This potential gives a valid description of the interactions between colloids in index-matched solvents Leunissen et al. (2005); Bier et al. (2010). For a given colloid volume fraction , the equilibrium behavior of colloids interacting by this potential can be described as a function of two parameters: , where is the well depth, and , where is the Debye screening length. At variance with Refs. Leunissen et al. (2005); Sanz et al. (2007, 2008), we consider a quite dilute system, with a volume fraction , which is below the percolation threshold.

Our simulations show that aggregates corresponding to metastable colloidal crystals form with very high probability under appropriate conditions. While in the case of high volume fractions treated in Ref. Sanz et al. (2007) the metastable phase was disordered, in our dilute system we find that a very well ordered metastable phase form. Specifically, we show that crystallites with the sodium-chloride (NaCl) structure grow for values of at which the cesium-chloride (CsCl) crystal structure is the most stable from the thermodynamic point of view. The formation of the NaCl-type crystallites is rationalized in the framework of the two-step mechanism of crystal nucleation Lutsko and Nicolis (2006); Vekilov (2010); Schilling et al. (2010), in which a dense colloidal fluid acts as a precursor for crystal nucleation. We show specifically that the type of short-range order in this dense fluid is crucial for determining which crystal phase is going to be formed. We show also that the NaCl phase formation is not related to energetically favorable structures of small clusters, the latter being either amorphous or with CsCl structure. Colloidal crystals with the NaCl structure were obtained with colloids of different sizes Vermolen et al. (2009). Our results demonstrate these crystals can grow in systems with balanced size and charge. We analyze also the importance of hydrodynamic effects in the nucleation of the NaCl phase by performing Stochastic Rotation Dynamics - Molecular Dynamics (SRD-MD) simulations Malevanets and Kapral (2000); Padding and Louis (2006); Tomilov et al. (2012), which show that these effects do not play a major role.

## Ii Model and Methods

In our simulations the interaction potential between the colloids is of the form

(1) |

where is a constant, which gives the value of the potential at contact, and the sign of the charge depends on the colloid and being either of type A or B. In the following we use . This expression of the potential is valid for . For , the hard-wall repulsion is approximated by a power law increasing as Sanz et al. (2008). We have checked that, in infinite crystals, the Madelung energy favors the CsCl phase over the NaCl phase for all values of , as pointed out in Ref. Leunissen et al. (2005). The Madelung energy difference between the two phases decreases with decreasing . In Ref. Leunissen et al. (2005) it was also shown that, at finite temperature, entropic contributions increase the stability of the NaCl bulk phase, which should however become more stable than the CsCl phase only below for of the magnitude considered here.

In the following, we use three different simulation techniques, Brownian Dynamics (BD) simulations, Basin Hopping (BH) global optimization searches, and SRD-MD simulations.

Our BD simulations are performed by numerically solving the Langevin equation by means of the algorithm described in Ref. Piechowiak et al. (2012). We choose the dimensionless time unit as , where is the diffusion coefficient of isolated colloids in the suspension. is the time taken by an isolated colloid to diffuse across a distance equivalent to its diameter. is given by Einstein’s relation, in which the solvent friction and the colloid mass are chosen as in Ref. Piechowiak et al. (2012). The Langevin equation is solved by using time steps in the range . The number of colloids in the simulations is 200, half being positively and half negatively charged. The colloids are initially placed in a cubic box of volume with random positions. Periodic boundary conditions are applied. Simulations are stopped after a time , which is sufficient to let all colloids form a single aggregate of crystalline structure. For comparison, in the system of Ref. Piechowiak et al. (2012), this would amount to s.

The lowest-energy structures of colloidal clusters have been searched for by a global optimization procedure employing the BH method Wales and Doye (1997), which has proven to be very efficient in the global optimization of atomic clusters. This method consists of a Metropolis Monte Carlo procedure in which the energy is locally minimized at each step. BH has been applied recently also to colloidal clusters in Ref. Cerbelaud et al. (2010b), where the computational procedure is described (see also Rossi and Ferrando (2009) for a more complete description of the method). We note here that we have used three different kinds of elementary moves to generate new configurations in the Monte Carlo procedure. These are the Brownian move, the shake move Rossi and Ferrando (2009), and a combined move, made of an uniform compression, in which the whole cluster is compressed by a factor 0.8 - 0.9, followed by a mild shake, i.e. of small random displacements of all colloids around their positions. After each move, the cluster is locally relaxed to generate the new configuration, which is then either accepted or refused according to the Metropolis rule. In all cases, 5 simulations of 10 steps each have been performed.

Hydrodynamic effects have been taken into account by means of SRD-MD simulations. SRD-MD is a combination of Stochastic Rotation Dynamics and standard Molecular Dynamics that allows taking into account long range hydrodynamic effects in the study of colloidal suspensions Malevanets and Kapral (2000); Padding and Louis (2006). This technique has been chosen for its relative efficiency and low computational cost, to check whether hydrodynamic effects could affect the nucleation process in our system. SRD-MD is based on the use of fake fluid particles that mimic the fluid properties (essentially its momentum). The simulation box is divided into smaller cubic cells. In each cell the fluid particles exchange momentum thanks to a rotation of their velocities. For an accurate description of the method we refer to the literature Malevanets and Kapral (2000); Padding and Louis (2006); Tomilov et al. (2012). Here we report only the choices of the most important parameters for our system. The linear size of the cells is taken as half of the radius of the colloids. For the key parameters of the SRD-MD simulations we made the same choice proposed by Padding and Louis Padding and Louis (2006): the number of fluid particles per cell equal to 5, the dimensionless mean free path equal to 0.1 and the rotation angle equal to 90. The parameters have been chosen to reproduce the same diffusion coefficient of isolated colloids in both BD and SRD-MD Tomilov et al. (2012). Recently it has been proved that these choices are able to correctly reproduce the essential features of the hydrodynamic interactions in a colloidal suspension Tomilov et al. (2012). For what concerns the energetic model, we have chosen the same Yukawa potential used in BD simulations for colloid-colloid interactions, while for colloid-fluid interactions we have chosen a soft-sphere potential

(2) |

for and 0 otherwise. We have chosen =12 and . The value of has been taken smaller than the colloidal radius () to avoid spurious depletion attractions between colloids Padding and Louis (2006).

## Iii Results and Discussion

Let us consider the aggregation of the colloids for , at which the NaCl crystal phase is never found to be the most stable Leunissen et al. (2005); Bier et al. (2010) from the thermodynamic point of view. The results concerning final aggregates in BD simulations (see Table 1) show that in a large majority of cases (15 in 20 simulations) structures of NaCl type are formed. The remaining simulations give either mixed or CsCl structures, the latter being found only in two cases. Representative snapshots of NaCl, mixed and CsCl structures are shown in Fig. 1. Also the simulations for give a large majority of NaCl-type aggregates. As increases the CsCl phase becomes more and more energetically favorable. Therefore, we expect that there should be a critical value at which the probability of forming both types of aggregates is the same. From the results in Table 1, it turns out that is close to 3.3, a value which is significantly larger than the at which the bulk NaCl and CsCl phases coexist in thermodynamic equilibrium Leunissen et al. (2005). Increasing further leads to the formation of CsCl-type aggregates with an overwhelming probability. On the other hand, if is decreased to 2.55, NaCl-type crystallites are always formed.

NaCl | mixed | CsCl | ||
---|---|---|---|---|

5% | 2.55 | 10 | 0 | 0 |

5% | 3.00 | 15 | 3 | 2 |

8% | 3.00 | 8 | 1 | 1 |

5% | 3.30 | 10 | 4 | 6 |

5% | 3.75 | 1 | 0 | 9 |

5% | 5.00 | 0 | 0 | 10 |

5% | 9.00 | 0 | 0 | 10 |

What are the driving forces leading to the formation of these aggregates of the metastable NaCl bulk phase?

A possible explanation of this result may originate from the fact that the aggregation process starts with the formation of very small colloidal clusters that subsequently grow and coalesce, as shown in the simulation of the aggregation in a similar system Piechowiak et al. (2012). Therefore, the structure of the final aggregate may preserve memory of the equilibrium structures of small aggregates. In order to investigate this point we have searched for the lowest-energy structures of colloidal clusters in the size range from AB to AB. The main structural motifs singled out for the values reported in Table 1 are shown in Fig. 2, where a few selected sizes are considered (AB, AB, AB, AB, AB, AB, and AB). Depending on , the dominant motif may change as follows.

For small aggregates, from AB to AB, the lowest-energy structures either belong to the CsCl motif, as structures a1, b1, and c1 in Fig. 2, or are disordered (i.e. they are neither of CsCl- nor of NaCl-type, see a2, b2, c2 in Fig. 2). CsCl-type structures prevail at , while at smaller disordered structures are more favourable. From AB on, CsCl-type structures (d1, e1, f, g1 in Fig. 2) always prevail with a single exception which is found at only. For this , the best structure of AB is of NaCl type, being a 33 cube with one vertex missing (structure e2). However, with increasing cluster size NaCl structures become more and more unfavorable. In fact, at the next magic size for NaCl structures, i.e. at size 64, where a perfect NaCl cube AB can be formed (structure g1), the CsCl structure is lower in energy for all down to 2.55.

These results show that structures of CsCl type prevail already for small aggregates, indicating that there are no evident effects that favor the NaCl structures at small sizes with respect to what is found in bulk crystals. Therefore the explanation of the growth of NaCl-type aggregates originates from a different effect, which is better understood by analyzing in details the aggregation process in the simulations.

To this end, we define two order parameters that are able to discriminate between the different crystalline phases. These parameters (referred to as and in the following) are defined as follows. is the fraction of colloids having more than 6 first neighbors. This indicates the possible presence of CsCl-type short range order, because in NaCl-type aggregates the number of first neighbors of a given colloid cannot be larger than 6. On the other hand, is related to the configuration of second neighbors, being defined as

(3) |

where is the distance between colloids and ( and being of the same type), and are the average equilibrium distance of second neighbors in the NaCl and CsCl structures, and takes into account the effect of vibrations (in our simulations we choose , a value that is more than three times smaller than ). Positive and negative values of respectively indicate the prevalence of either NaCl- or CsCl-type short-range order in the aggregate. In perfect NaCl and CsCl infinite crystals is close to 1 and -1, respectively.

Let us consider the typical behavior of a simulation for (see Fig. 3, crosses). After a rather short time, essentially all colloids attach to some aggregate and, as in most cases, a single aggregate is quickly formed (typically for ). stays very close to 0, showing that the probability of configurations in which a colloid has more than six first neighbours is very low. On the other hand, increases from the initial value of 0 to a positive value , indicating that some NaCl-type short-range order is present. However, if we look at the aggregate itself, we find configurations as those reported in the snapshot a1 of Fig. 4. The aggregate does not show any long-range order, but a liquid-like structure, in which some colloids surrounded by a NaCl-type short-range order can be singled out. After a lengthy evolution in this disordered phase, there is a rather sudden increase of , which jumps to a value of 0.4, corresponding to the transformation of the disordered aggregate into an ordered NaCl-type crystallite, as the one shown in snapshot a2 of Fig. 4.

This behavior can be compared to the case of . Again, a single aggregate is rather quickly formed. The aggregate does not show long-range order, but a somewhat prevailing CsCl-type short-range order (see b1 in Fig. 4) corresponding to a negative , and to a small but significantly positive . After some time, the sudden transition to the ordered crystallite takes place, but this time the structure is of CsCl type (b2 in Fig. 4).

This analysis shows that the final structure of the crystallite is therefore well described by the two-step mechanism of crystal nucleation Lutsko and Nicolis (2006); Vekilov (2010). Starting from diluted configurations () dense aggregates form with a metastable liquid-like structure which survives for some time. However, a transition to an ordered structure finally takes place, and the type of ordered structure is correlated, with very high probability, with the features of short-range order in the metastable liquid, as described for example by the parameter . Specifically we find that the final structure is almost always of NaCl type for disordered aggregates in which , while for the CsCl structure is found in most cases. This means that the final ordered structure can be predicted with high probability by inspecting the relative positions of the colloids in the disordered state. We note that this is not always the case in colloid crystallization, as demonstrated for hard-sphere glasses in Ref. Sanz et al. (2011).

The two-step crystallization process in hard-sphere glasses has been studied by Monte Carlo simulations by Schilling et al. Schilling et al. (2010), at high volume fraction . They found a rather gradual development of crystallinity from the disordered state as the simulation goes on, while in our case crystallization is rather sudden.

In the disordered phase, for , the range of the repulsion between colloids with the same charge is sufficiently long to prevent the occurrence of configurations in which more than six colloids of the same sign can surround one colloid of the opposite sign. These configurations would be indeed energetically favorable, because the CsCl phase is still the most stable from this point of view, but they are very difficult to reach in an aggregation process starting from a dilute volume fraction. This is the cause of the formation of crystallites of the metastable NaCl phase. Once formed, these crystallites can survive for long lifetimes.

Rather surprisingly, we note that the formation of NaCl-type crystallites is indeed easier when this bulk phase is metastable. We have in fact tried aggregation simulations for and , a situation in which the NaCl bulk phase is more stable than the CsCl phase due to entropic effects Leunissen et al. (2005). However, the long range of the repulsion between particles with the same charge renders the aggregation itself quite difficult, so that we were not able to observe any crystallization within .

Finally we discuss the importance of hydrodynamic effects by comparing the results of BD and SRD-MD.

The SRD-MD is much heavier computationally than BD so that it is very cumbersome to simulate the entire process of two-step aggregation and reorganization, because the nucleation of the ordered phase from the liquid can occur after very long times. For this reason we have decided to restrict the comparison to the analysis of the crucial parts of the process. First of all, we have started both BD and SRD-MD from initial randomly dispersed non-overlapping configurations and analyzed the initial aggregation stage, leading to the formation of the disordered liquid-like phase, to check whether this phase forms on the same time scale and presents the same properties in both BD and SRD-MD. Then we have analyzed the formation of the ordered crystal, starting from configurations in which a NaCl seed is already well developed. These configurations have been taken from previous BD simulations. The configurations with NaCl seed have been let evolve by BD and by SRD-MD. The results for , averaged over 4 simulations in each case, are reported in Fig. 5, where the behavior of versus time is shown. The figure shows that both BD and SRD-MD produce the same kind of behavior. The liquid-like disordered phase is formed on comparable time scales in both cases, and it presents the same structure, with oscillating around 0.035-0.04. Also the growth of the NaCl seeds in the liquid occurs on the same time scale (in SRD-MD it seems only somewhat slower than in BD), producing the same kind of final structures. These results indicate that hydrodynamic effects do not play a major role in the two-step formation process of the ordered NaCl crystallites.

## Iv Conclusions

In this work we have analyzed the heteroaggregation of oppositely charged colloids by means of a combination of different simulation techniques, including Brownian Dynamics, global optimization searches and Stochastic Rotation Dynamics - Molecular Dynamics. Aggregation in dilute systems has been considered. Our results show that colloidal aggregation can be dominated by kinetic effects which can lead to the formation of very well ordered metastable crystals. These crystals significantly differ from the disordered metastable fcc phases that were previously found in Ref. Sanz et al. (2007) at higher volume fractions.

Our simulations show that, if the range of the potential is sufficiently long (), kinetic effects cause the formation of crystallites of the metastable NaCl-type bulk phase, which thus can be obtained also for binary colloids of same size and balanced charge. These crystallites form within a two-step mechanism, in which the first step is the aggregation of a dense metastable liquid whose short-range order has prevailing NaCl features. This triggers the transition to ordered NaCl-type crystallites. Therefore the most likely outcome of the aggregation process is not the crystal structure of greater thermodynamic stability, but the structure whose short-range order is closer to the short-range order in the metastable liquid. This indicates a mechanism of general character which is likely to be relevant for a variety of nucleation processes. The analysis of the most stable structures of small clusters has also shown that the formation of the NaCl phase is not related to the structures of very small aggregates, that are either amorphous or of CsCl structure. Our results show also that NaCl-type crystallites more easily form when the NaCl bulk phase is metastable or only marginally stable. Finally, the inclusion of hydrodynamic effects in the simulations has not produced major effects on the aggregation process. All these findings are important for designing strategies to obtain the desired structures of colloidal crystals in practical cases.

D. B. acknowledges support from the Programme VINCI 2011 of the French-Italian University on the subject “Studio dell’eteroaggregazione di particelle colloidali ceramiche tramite tecniche di ottimizzazione globale”.

### Footnotes

- Corresponding author, e-mail ferrando@fisica.unige.it

### References

- M. E. Leunissen, C. G. Christova, A. P. Hynninen, C. P. Royall, A. I. Campbell, A. Imhof, M. Dijkstra, R. van Roij, and A. van Blaaderen, Nature 437, 235 (2005).
- J. M. López-López, A. Schmitt, A. Moncho-Jordá, and R. Hidalgo-Álvarez, Soft Matter 2, 1025 (2006).
- S. Rollie and K. Sundmacher, Langmuir 24, 13348 (2008).
- N. I. Lebovka, Adv. Polym. Sci. (2012), DOI: 10.1007/12_2012_171.
- P. Bartlett and A. I. Campbell, Phys. Rev. Lett. 95, 128302 (2005).
- A.-P. Hynninen, M. E. Leunissen, A. van Blaaderen, and M. Dijkstra, Phys. Rev. Lett. 96, 018303 (2006).
- E. Sanz, C. Valeriani, T. Vissers, A. Fortini, M. E. Leunissen, A. van Blaaderen, D. Frenkel, and M. Dijkstra, J. Phys.: Condens. Matter 20, 494247 (2008).
- M. A. Piechowiak, A. Videcoq, R. Ferrando, D. Bochicchio, C. Pagnoux, and F. Rossignol, Phys. Chem. Chem. Phys. 14, 1431 (2012).
- E. R. Russell, J. Sprakel, T. E. Kodger, and D. A. Weitz, Soft Matter 8, 8697 (2012).
- M. Bier, R. van Roij, and M. Dijkstra, J. Chem. Phys. 133, 124501 (2010).
- G. Pavaskar, S. Sharma, and S. N. Punnathanam, J. Chem. Phys. 136, 134506 (2012).
- E. Sanz, C. Valeriani, D. Frenkel, , and M. Dijkstra, Phys. Rev. Lett. 99, 055501 (2007).
- M. Cerbelaud, A. Videcoq, P. Abélard, C. Pagnoux, F. Rossignol, and R. Ferrando, Soft Matter 6, 370 (2010a).
- M. A. Piechowiak, A. Videcoq, F. Rossignol, C. Pagnoux, , C. Carrion, M. Cerbelaud, and R. Ferrando, Langmuir 26, 12540 (2010).
- R. J. Allen, P. B. Warren, and P. R. ten Wolde, Phys. Rev. Lett. 94, 018104 (2005).
- J. F. Lutsko and G. Nicolis, Phys. Rev. Lett. 96, 046102 (2006).
- P. G. Vekilov, Nanoscale 2, 2346 (2010).
- T. Schilling, H. J. Schöpe, M. Oettel, G. Opletal, and I. Snook, Phys. Rev. Lett. 105, 025701 (2010).
- E. C. M. Vermolen, A. Kuijk, L. C. Filion, M. Hermes, J. H. J. Thijssen, M. Dijkstra, and A. van Blaaderen, Proc. Natl. Acad. Sci. 106, 16063 (2009).
- A. Malevanets and J. Kapral, J. Chem. Phys. 112, 7260 (2000).
- J. T. Padding and A. A. Louis, Phys. Rev. E 74, 031402 (2006).
- A. Tomilov, A. Videcoq, T. Chartier, T. Ala-Nissila, and I. Vattulainen, J Chem. Phys. 137, 014503 (2012).
- D. J. Wales and J. P. K. Doye, J.Phys. Chem. A 101, 5111 (1997).
- M. Cerbelaud, R. Ferrando, and A. Videcoq, J. Chem. Phys. 132, 084701 (2010b).
- G. Rossi and R. Ferrando, J. Phys. Cond. Mat. 21, 084208 (2009).
- E. Sanz, C. Valeriani, E. Zaccarelli, W. C. K. Poon, P. N. Pusey, and M. E. Cates, Phys. Rev. Lett. 106, 215701 (2011).