Molecular dynamics simulation of binary hard sphere colloids near the glass transition

Abstract

This study reports on molecular dynamics simulations of concentrated binary hard sphere colloids near the glass transition. We present the three-dimensional mean square displacements for systems that increase in density, into the glass transition region. The molecular dynamics simulations probe wider time scales than the experimentally assessable range. They cover the short time scale that corresponds to ballistic motion, the intermediate scale that corresponds to caged Brownian motion, and the long time diffusion described by the Einstein–Smoluchowski limit. We report results for α₂, the non-Gaussian parameter, which captures the details of the mean square displacements around the time that, on average, rearrangements of the particle cages occur. The dynamic temporal correlations were examined in terms of the four point susceptibility. The implications for ageing are also discussed.

---

1. Introduction

Hard sphere colloids are an excellent model for studying the cooperative behaviour in systems with simple and well-defined interactions. They can help provide important insights about the structural and transport properties of suspensions including Brownian dynamics, rheology, phase behaviour and stability. Binary mixtures of hard spheres were shown to crystallize in the same lattice geometry as metal alloys if the ratio between the two types of particles is identical to that of the metal atoms.^1–4 These extremely important experiments showed that excluded volume effects (i.e., hard sphere interactions) are dominant over a wide range of length-scales, and implied that the physical conclusions and insights obtained from studying hard sphere colloids may have a much greater validity, than previously expected.

The dynamics of hard spheres at high densities is dominated by multiparticle correlations. As the overall volume fraction increases the individual particle motions become increasingly hindered until the system approaches either a phase-transition or a disordered glassy state. The latter often results when the system is not a monodisperse colloid. In a recent paper⁵ Narumi et al. presented results for the dynamics of a binary colloid system at densities that are close to the glass transition. The particle sizes and composition of the binary colloid were specifically chosen so as to prevent the system from easily crystallizing. The individual particle motions were tracked by confocal microscopy. The main observations in this paper were the spatial and temporal heterogeneities as the system approaches glass transition. This was attributed to a cooperative motion of whole domains within the suspension. Similar conclusions were also reached by other authors.^6–12

In the present paper we report on a Molecular Dynamics (MD) study of a hard sphere mixture that is similar to the dense suspensions studied experimentally by Narumi et al.⁵ The MD approach allows us an examination of a greater time interval and provides a more detailed description of the particle motions, both at very short times when the particle motion is in the ballistic regime, as well as at long times when the particle breaks out of the “cage” of the surrounding hard spheres. The MD computations do not account for the presence of the solvent, which is qualitatively correct when the dynamics is dominated by particle–particle instead of particle–solvent molecules collisions. Replacing the MD by Brownian dynamics that includes the solvent (through a friction kernel) is problematic for systems near glass transition as discussed at the end of this article.

The paper is organized as follows, the next section briefly describes the details of the MD procedure, Section 3 presents the Results and discussion and Section 4 summarizes the Conclusions.

2. Molecular dynamics procedure

The system of interest is chosen to be like that of Narumi et al.⁵ who studied a binary colloid suspension. The components making up the colloid consisted of colloidal hard spheres of diameters 3.10 and 2.36 μm, or a size ratio of small over large of 0.7613. Both the large (L) and the small (S) components displayed a polydispersity of 5%. To mimic the polydispersity we selected to subdivide the large component into three very similar sizes, and similarly the small component into three sizes, such that the polydispersity of each component was 5%. In what follows, we will usually simply refer to the large and small components, and not make any distinctions between the sub-classes. The number ratio, small to large (N_S/N_L), used is that of Narumi et al.,⁵ namely 1.56.

The MD procedure is an event-driven approach based on binary collisions between the molecules such that both the momentum and kinetic energy are conserved. When two particles i and j with initial velocities v_i and v_j collide they obtain new velocities v^′_i and v^′_j that are given by

where r_ij = r_i − r_j and v_ij = v_i − v_j.

No explicit solvent molecules are present in the computations and thus all collisions take place between the colloidal particles. This is an approximation, of course, but this choice does not lead to a qualitatively different result. The Brownian motion in our simulation is entirely due to collisions between the particles and does not account for the random collisions with the solvent molecules that are present in the experimental suspension. The absence of an explicit low molecular weight solvent overemphasizes the very short-time ballistic motion, but the long time purely diffusion behaviour should be very similar since it is driven by a fluctuation–dissipation balance, which is similar the presence bath of small molecules. An example of this behaviour can be found in the work of Scala et al.,¹³ who introduced event-driven Brownian dynamics for hard spheres. They demonstrated that the density dependence of the diffusion constant in the liquid phase is, in fact, the same for molecular dynamics and Brownian dynamics when the diffusion coefficient is scaled by its appropriate low-density value (see also ref. 14–17).

Two sets of computations were performed. The first one was for a total of 40 000 hard spheres, while the second was for 4000. The larger number of particles allows for obtaining better statistics with accurate averages for the mean square displacements. The lower number of particles was used to capture spatial and temporal heterogeneities in the collective dynamics of the particles.

3. Results and discussion

3.1. Equilibrium structure and radial distribution function

The equilibrium structure of liquids and complex fluids in general is characterized by the radial distribution function (RDF), which illustrates the short range liquid order induced by the repulsive interactions whose amplitude decays at longer length scales.¹⁸ Our system of interest, consisting of a binary hard sphere mixture centred on two distinct particle sizes, gives rise to three radial distribution functions according to the particle types, i.e., 11, 12, and 22. A typical result for the RDF is shown in Fig. 1. It shows the small–small, small–large and large–large distributions functions. As mentioned above, strictly speaking our simulations consist of six different size species, to account for the experimentally observed polydispersity. We present the RDFs for the two main sizes corresponding to the peaks in the distribution. The rest of the hard spheres, of less common sizes that are present in the mixture, however, still have an effect on the shape of each RDF by introducing the lower satellite peaks in the vicinity of the main ones. The main peaks positions correspond to sum of the radii of the particular interacting species. The distance on the horizontal is scaled with the diameter of the larger particles in our simulations. The RDFs shown in Fig. 1 can be directly compared to the ones reported by Narumi et al. (see Fig. 2 of ref. 5). There is excellent agreement between peak locations and peak heights.

Fig. 1 Radial distribution functions for small–small (2–2, dotted line), small–large (2–1, solid line) and large–large (1–1, dashed line) cases. Only the RDFs two main sizes corresponding to the peaks in the distribution are shown. The distance on the horizontal is scaled with the diameter of the larger particles σ₁.

3.2. Particle dynamics and mean square displacements

MD allows tracing a significantly greater range to times compared to experiments. There is no limit for the shortest times, and the high end, limited by computer resources, is well beyond the times that can be assessed by experiment.⁵Fig. 2 presents results for simulations for a total of 40 000 particles. Each run consists of about 2.5 × 10¹⁰ collisions. All six species (the main two components and the satellite particles that account for the polydispersity) are plotted. Hence, the different curves in each figure correspond to different size fraction of the mixture starting from the smallest (top curves) and ending with the largest (lowest curves). The initial linear slopes in all plots shown in the figure correspond to ballistic motion of each particle which is described by the equationwhere r_i is the displacement of particle i with mass m_i, kT is the thermal energy and t is time. In this regime most particles simply move freely along straight lines. At longer times the ballistic motion no longer characterizes the particle motion, instead multiple collisions with neighbouring particles (or molecules) lead to Brownian motion, which for moderately diluted systems is given by

〈r_i²〉 = 6Dt (3.2)

with D denoting the diffusion coefficient. A complete analysis describing both regimes as well as the transition between the two types of motion was given by Uhlenbeck and Ornstein¹⁹ (see also ref. 20–22 or 23). Hence the Uhlenbeck–Ornstein process predicts two distinct regimes of particle motion: ballistic and Brownian. Therefore, when plotted in logarithmic scale the 〈r²〉 vs. t dependence consists of two straight lines with different slopes bound by a smooth transition region without any inflection. Our results in Fig. 2 however do not fully conform to this behaviour. The data for ϕ = 0.53 (see Fig. 2a) are very close to the Uhlenbeck–Ornstein model, but a more careful inspection reveal a that there is a trend toward formation of an inflection point at intermediately long times. This is even more pronounced at higher volume fractions (see Fig. 2b–d), where one may distinguish three different regimes of particle motion. The free particle ballistic motion shows that 〈r²〉 rapidly increases with time (see Fig. 3a). At longer times, however, the particle dynamics is dominated by large number of inter-particle collisions. The particle exhibits diffusion-like dynamics, but 〈r²〉 does not increase much because of the obstruction due to a “cage” of surrounding particles (see Fig. 3b). This cage “traps” the tracer particle and the range if the motion is of the order of the cage inner dimensions. At sufficiently long times the particle would eventually break free from this particular cell and move to other locations in space that correspond to larger displacements and 〈r²〉 will further increase when that happens (see Fig. 3c). That is why at long times the mean square displacement 〈r²〉 is greater. However, the higher the volume fraction the longer is the time that a particle spends trapped in the cage from the surrounding neighbours. This is evident from Fig. 2d–f. In fact for ϕ = 0.62 (Fig. 2f) and higher, no increasing trend in 〈r²〉 was obtained for the times covered by the simulations. If the system has undergone a glass transition, each particle may be forever locked within the cage of the surrounding nearest neighbours.

Fig. 2 Mean square displacement vs. time for a hard sphere mixture. (a) ϕ = 0.53; (b) ϕ = 0.57; (c) ϕ = 0.58; (d) ϕ = 0.59; (e) ϕ = 0.60; (f) ϕ = 0.62. The insets in (e) and (f) show the behavior at very long times.

Fig. 3 A sketch illustrating the motion of a particle in a dense suspension. (a) Ballistic motion (no collisions); (b) short scale diffusion within a “cage”; (c) hopping between “cages”.

As the dispersion volume fractions increases the individual particle diffusivity becomes hindered, however sporadic collective motions may occur. They are related to simultaneous coordinated motions of domains containing large number of particles. Such fluctuation driven occurrences are not visible in the dynamics 40 000 particles because of the averaging over such a large system. However smaller number of particles show occasional jumps in the 〈r²〉 vs. time plots. Fig. 4 presents results obtained from a MD computation of the mean square displacement against time for binary hard sphere mixtures containing 4000 particles. The total number of collisions was again equal to 2.5 × 10¹⁰. Because of the smaller number of particles fluctuation-driven collective motions are not smoothed out by averaging and can be observed. Fig. 4a corresponds to ϕ = 0.62. Four distinct jumps in the mean square displacement 〈r²〉 are observed. The jumps propagate across all particles irrespective of their size. This effect is due to a cooperative motion of a whole domain in the process of rearranging the suspension structure.²⁴ Similar motions are detected at higher volume fractions. Fig. 4b presents data for ϕ = 0.64. This time three sudden increments in the 〈r²〉 vs. time plot are clearly visible. The relative magnitude of the jumps is smaller compared to the case shown in Fig. 4a. The reason is in the higher volume fraction, which results in less available space for particle motions and collective domain rearrangements. Still the sudden jumps in the mean square displacement exceed the noise level of the curves.

Fig. 4 Mean square displacements of each hard sphere component at long times and high densities. (a) ϕ = 0.62; (b) ϕ = 0.64.

3.3. Non Gaussian effects on the suspension dynamics

At high densities the particle movement is characterized by two broad types of movement: limited movement within a cage of neighbors, and more significant excursions that accompany the re-arrangement of those cages. If one analyzes the distribution of displacements on a timescale that corresponds to these re-arrangements, one tends to find a rather broad distribution. It has been suggested^25,26 that this broadening with respect to a Gaussian distribution can be made quantitative by the use of a so-called non-Gaussian parameter,where Δx = x(t) is the displacement in one dimension after some time t. The non-Gaussian parameter is a measure of the excess kurtosis of a distribution. For a perfect Gaussian distribution α₂ = 0. A value larger than zero is an indication that the tails of the distribution are more prominent, i.e., that a large displacement is more likely than a Gaussian distribution would predict.

In Fig. 5 we show the results for α₂ as a function of time, averaging over all particles and over multiple time origins. We display the results for both the small and the large particles. It is clear that at higher densities there is a development of a distinct maximum in α₂. This maximum shifts to larger time t as he volume fraction increases. This trend is the same for small and large particles. The maximum in α₂ as a function of time occurs earlier for the small particles (at any given volume fraction).

Fig. 5 Kurtosis parameter for the two main size species (corresponding to the peaks of the size distribution – see the text).

It is interesting to point that similar dynamic in homogeneities were also observed in binary mixtures of soft spheres.²⁴

3.4. Aging effects

As the density approaches the glass transition the dynamics of the system slows down considerably, as is apparent from the mean square displacements. At (and beyond) the glass transition density the properties of the system become increasingly dependent on the history of the sample. That is, a measurement of the mean square displacements can be expected to vary with time.²⁷ In Fig. 6 we plot the mean square displacements at early and at late times, at volume fractions ϕ = 0.59 and 0.60. Both these samples were generated by compressing a system equilibrated at ϕ = 0.54. The displacements were then collected for a period of about 100 time units (mσ₂/kT)^1/2, which is long enough to see the onset of the diffusive regime. And then, similarly, for the intervals [100, 200] and at the “late” time [200, 300]. At both volume fractions the early times are markedly different from the subsequent time intervals. Comparing the results for ϕ = 0.59 and 0.60 the aging is somewhat longer for the higher density, as expected.

Fig. 6 Mean square displacements at early and at late times for (a) ϕ = 0.59 and (b) 0.60 showing the effects of aging on the overall particle dynamics.

Following Pusey et al.²⁸ we fitted out data to the expression

ln(〈r²〉/6t) ∝ γ ln(|ϕ − ϕ_g|), (3.4)

where ϕ_g is the glass transition volume fraction. The theoretical value for γ for hard spheres is 2.15.²⁹ We found that our data can be fitted with γ = 2.15 if ϕ_g = 0.592. This value is close to the one reported by Pusey et al. (ϕ_g = 0.585)²⁸ for monodisperse hard spheres. This may indicate that the location of the glass transition is only slightly sensitive to the bidispersity of the hard sphere system.

3.5. Four point susceptibility

We have determined the four point correlation (susceptibility) function χ₄, which represents the fraction of particles with correlated dynamics at given moment in time. It is defined by^30,31where

w(|r_ij − μ_j|) is an overlap function (see ref. 31 for more details).

The results for χ₄ are shown in Fig. 7. The total density is ϕ = 0.58. The curves represent the results for the larger of the two main particle species and different values for the characteristic length-scale ΔL describing the spatial inhomogeneity.^30,31 All curves show that there is no correlation for times that are either equal close to zero or very long. Hence that the dynamic temporal correlations pass through a maximum. There is a certain time that corresponds to strongest dynamic correlation, which is defined by the maximum of the curves. After that the motions becomes less correlated until χ₄ goes down to zero for long times.

Fig. 7 Four point susceptibility function for ϕ = 0.58. The different curves correspond to different displacement thresholds ΔL.

4. Summary and conclusions

In summary we performed MD analysis of the behaviour of binary hard sphere dispersions at high volume fractions. The analysis showed that as the systems approaches the glass transition three distinct time ranges can be identified: (i) short times where the particles exhibit ballistic (pre-collision) motion,²³ (ii) intermediate times corresponding to the small length-scale movements of a particle in a “cage” formed by the surrounding nearest neighbours, and (iii) long time diffusion corresponding to the Einstein–Smoluchowski limit.^32–34 The mean square displacement in case (i) is proportional to t², does not change with time for case (ii) and is proportional to t in the long-time diffusion limit, case (iii). These are shown in Fig. 2.

Our computations showed that at high densities and long times a cooperative motion of domains of particles can occur in an attempt to reach a more optimal arrangement. These appear as jump-wise increases in the mean square displacements when plotted against time – see Fig. 4.

The MD data confirmed the experimental observation⁵ that the random diffusion process has a significant non-Gaussian component (Fig. 5). This presents difficulties to in describing the complete dynamics in the framework the theory of anomalous diffusion.^23,35 Hence, other approaches are often used.^36,37 Another difficulty follows from the observation that the system ages with time (see Fig. 6 and ref. 5) and therefore is non-ergodic and not in true thermodynamic equilibrium. This implies that fluctuation–dissipation relationships relating the friction to moments of the random force distribution can only be used as quasi-static approximations or not at all. We argue that MD simulations offer an excellent tool for analysis of colloids near the glass transition.

Acknowledgements

The authors are grateful to Dr Christine Roberts for bringing ref. 5 to our attention. This work was funded by NSF (CBET 0844645) and by Sandia National Laboratories' LDRD program. Sandia National Laboratories is a multi-program laboratory managed and operated by Sandia Corporation, a wholly owned subsidiary of Lockheed Martin Corporation, for the U.S. Department of Energy's National Nuclear Security Administration under contract DEAC04-94AL85000. Finally the authors wish to thank the reviewer for the helpful comments that lead to a much improved final version of this article.

Notes and references

1. P. Bartlett, J. Phys.: Condens. Matter, 1990, 2, 4979–4989.
2. P. Bartlett and R. H. Ottewill, Langmuir, 1992, 8, 1919–1925.
3. P. Bartlett, R. H. Ottewill and P. N. Pusey, Phys. Rev. Lett., 1992, 68, 3801–3804.
4. P. Bartlett and R. H. Ottewill, J. Chem. Phys., 1990, 93, 1299–1312.
5. T. Narumi, S. V. Franklin, K. W. Desmond, M. Tokuyama and E. R. Weeks, Soft Matter, 2011, 7, 1472–1482.
6. G. Adam and J. H. Gibbs, J. Chem. Phys., 1965, 43, 139–146.
7. M. D. Ediger, Annu. Rev. Phys. Chem., 2000, 51, 99–128.
8. E. Hempel, G. Hempel, A. Hensel, C. Schick and E. Donth, J. Phys. Chem. B, 2000, 104, 2460–2466.
9. R. Richert, J. Phys.: Condens. Matter, 2002, 14, R703–R738.
10. H. Sillescu, J. Non-Cryst. Solids, 1999, 243, 81–108.
11. E. Donth, J. Polym. Sci., Part B: Polym. Phys., 1996, 34, 2881–2892.
12. C. Donati, J. F. Douglas, W. Kob, S. J. Plimpton, P. H. Poole and S. C. Glotzer, Phys. Rev. Lett., 1998, 80, 2338–2341.
13. A. Scala, T. Voigtmann and C. D. Michele, J. Chem. Phys., 2007, 126, 134109.
14. M. Tokuyama, H. Yamazaki and Y. Terada, Phys. Rev. E: Stat., Nonlinear, Soft Matter Phys., 2003, 67, 062403.
15. T. Gleim, W. Kob and K. Binder, Phys. Rev. Lett., 1998, 81, 4404–4407.
16. G. Szamel and E. Flenner, Europhys. Lett., 2004, 67, 779–785.
17. L. Berthier, G. Biroli, J. P. Bouchaud, W. Kob, K. Miyazaki and D. R. Reichman, J. Chem. Phys., 2007, 126, 184503.
18. D. A. McQuarrie, Statistical Mechanics, University Science Books, Sausalito, 2000.
19. G. E. Uhlenbeck and L. S. Ornstein, Phys. Rev., 1930, 36, 823–841.
20. N. G. van Kampen, Stochastic Methods in Physics and Chemistry, Elsevier, New York, 2 edn, 1992.
21. C. W. Gardiner, Handbook of Stochastic Methods, Springer, New York, 1985.
22. R. M. Mazo, Brownian Motion: Fluctuations, Dynamics and Applications, Clarendon Press, Oxford, 2002.
23. W. T. Coffey, Y. P. Kalmykov and J. T. Waldron, The Langevin Equation, World Scientific, London, 2nd edn, 2004.
24. H. Miyagawa, Y. Hiwatari, B. Bernu and J. P. Hansen, J. Chem. Phys., 1988, 88, 3879–3886.
25. A. Rahman, Phys. Rev., 1964, 136, A405–A411.
26. B. Bernu, J. P. Hansen, Y. Hiwatari and G. Pastore, Phys. Rev. A, 1987, 36, 4891–4903.
27. G. Brambilla, D. El Masri, M. Pierno, L. Berthier and L. Cipelletti, Phys. Rev. Lett., 2009, 102, 085703.
28. P. N. Pusey, E. Zaccarelli, C. Valeriani, E. Sanz, W. C. K. Poon and M. Cates, Philos. Trans. R. Soc., A, 2009, 367, 4993–5011.
29. W. Gotze and L. Sjogren, Rep. Prog. Phys., 1992, 55, 241–376.
30. C. Donati, S. C. Glotzer and P. H. Poole, Phys. Rev. Lett., 1999, 82, 5064–5067.
31. S. C. Glotzer, V. N. Novikov and T. B. Schrøder, J. Chem. Phys., 2000, 112, 509–512.
32. A. Einstein, Ann. Phys., 1905, 17, 549–560.
33. M. v. Smoluchowski, Ann. Phys., 1906, 21, 756–780.
34. M. v. Smoluchowski, Phys. Z., 1916, 17, 557–585.
35. E. Lutz, Phys. Rev. E: Stat., Nonlinear, Soft Matter Phys., 2001, 64, 051106.
36. F. Sciortino and P. Tartaglia, Adv. Phys., 2005, 54, 471–524.
37. A. Cavagna, Phys. Rep., 2009, 476, 51–124.

---