Hydrodynamic Correlations slow down Crystallization of Soft Colloids
aa r X i v : . [ c ond - m a t . s t a t - m ec h ] A p r Hydrodynamic Correlations slow down Crystallization of Soft Colloids
D. Roehm, ∗ S. Kesselheim, † and A. Arnold ‡ University of Stuttgart, Allmandring 3, 70569 Stuttgart, Germany (Dated: September 10, 2018)Crystallization is often assumed to be a quasi-static process that is unaffected by details of particletransport other than the bulk diffusion coefficient. Therefore colloidal suspensions are frequentlyargued to be an ideal toy model for experimentally more difficult systems such as metal melts. In thisletter, we want to challenge this assumption. To this aim, we have considered molecular dynamicssimulations of the crystallization in a suspension of Yukawa-type colloids. In order to investigatethe role of hydrodynamic interactions (HIs) mediated by the solvent, we modeled the solvent bothimplicitly and explicitly, using Langevin dynamics and the fluctuating Lattice Boltzmann method,respectively. Our simulations show a dramatic reduction of the crystal growth velocity due to HIseven at moderate hydrodynamic coupling. A detailed analysis shows that this slowdown is dueto the wall-like properties of the crystal surface, which reduces the colloidal diffusion towards thecrystal surface by hydrodynamic screening. Crystallization in suspensions therefore differs stronglyfrom pure melts, making them less useful as a toy model than previously thought.
PACS numbers: 64.70.D-, 47.57.E-, 64.75.Xc, 07.05.Tp
Crystallization and nucleation of undercooled meltsare often studied using model systems of charged col-loids in solution [10, 17, 18], such as polystyrene (PS)or polymethylmethacrylat (PMMA) spheres suspendedin water [7, 15, 26]. The solvent gives rise to hydrody-namic interactions (HIs) between the colloidal particlesthat are not present in, e.g., metal melts. The influenceof HIs on the dynamical properties of colloidal suspen-sions has been extensively studied in recent years [23, 24].L¨owen et al. [17, 19] showed that the ratio of the long-time to short-time self-diffusion coefficients has a univer-sal value along the fluid freezing line. Recent studies byPesche [25] and N¨agele [20] of quasi-2D dispersions showthat HIs have an impact on the self-diffusion functionin these soft-sphere suspensions. However, since crystalgrowth happens on much larger time scales, it is com-monly believed to be a quasi-static process that is unaf-fected by HIs. For example, in classical nucleation theory(CNT) [14], hydrodynamic interactions are assumed toenter only through the effective diffusion constant of theattaching colloids, which can be measured convenientlyin the bulk liquid. Also in most computer simulationstudies of nucleation, hydrodynamic interactions are ne-glected to avoid the high computational costs.In the following, we will show with the help of com-puter simulations that HIs do have a remarkable influenceon the dynamics of the crystallization in a colloidal sus-pension. Even at moderate coupling, we found the crystalgrowth velocity to be reduced by a factor of three. Sim-ilar findings have been reported by Schilling et al. usingdifferent simulations techniques [27].To investigate the influence of hydrodynamic interac-tions on crystal growth, we studied the crystallizationof Yukawa-type particles confined between two planarwalls. The influence of hydrodynamic correlations oncrystallization was evaluated by performing both sim- ulations including and excluding HIs, by employing afluctuating lattice Boltzmann method [6] and Langevindynamics [12], respectively. We prepare our systemsas an undercooled liquid and let the system crystallize.Due to the presence of the confining walls, the nucle-ation barrier is sufficiently small, so that we do notneed special rare event sampling techniques. All simu-lations were performed using the MD simulation package
ESPResSo [4, 5, 28].As interparticle potential we used a screened Coulombinteraction potential U ( r ) = l B k B T Q Q exp( − λ D r ) r , (1)where l B is the Bjerrum length, k B is the Boltzmannconstant, Q and Q give the charges on the interact-ing particles, r is distance between the particles. Therange of the potential is determined by the Debye-H¨uckelscreening length λ D . The static properties of a Yukawasystem can be characterized by two dimensionless param-eters [10] κ = wλ D and Γ = Q Q πǫ wk b T = Q Q l B w , (2)where w = (3 / (4 πρ )) / is the Wigner-Seitz radius ofthe crystal phase and ρ the particle density. Phase di-agrams of systems with Yukawa-type interactions havebeen calculated both by Monte Carlo simulations [21]and MD simulations [10, 11], which consistently foundthree regimes: a fluid phase and two different solid phaseswith BCC or FCC structure, respectively. For our simu-lations, we chose κ = 3 . γ , which is inverselyproportional to the single particle diffusion constant.In order to introduce hydrodynamic interactions with-out simulating the solvent explicitly, we used the latticeBoltzmann method [6] on a three-dimensional (3D) lat-tice with 19 velocity densities (D3Q19). In all reportedsimulations with HIs, we used a grid spacing of g = 1 . γ [2].Note that this friction constant is equivalent to the fric-tion in the Langevin dynamics, which is why we use thesame label γ . The no-slip boundaries modeling the fluid-wall interaction were realized by the link bounce backrule [33]. The strength of the hydrodynamic interactionscan be quantified by the hydrodynamic radius r H = k B T πηD , (3)where r H depends on the viscosity η of the fluid as wellas of the single particle diffusion coefficient D , whichis inversely proportional to the friction constant γ . Incontrast to other methods for including HIs, such asDPD [8, 13], LB methods allow the friction γ to be tunedindependently from the viscosity η of the fluid. The par-ticle mobility is not simply the inverse of γ , but alsodepends on the viscosity η and the lattice spacing of theLB grid due to feedback from the moving fluid [2, 31].Also, the back flow of the solvent introduces finite-sizeeffects in a system with periodic boundary conditions.The equations of motion of the Yukawa particles wereintegrated by a Velocity-Verlet integrator [12]. If nototherwise stated, the simulations were performed with16 ,
384 particles in a box of size 66 × ×
16 confinedby two planar walls located at x = 0 . x = 65 . dt = 0 .
01, the total simulation length 750 ,
000 timesteps. The same time step was also used for the LB fluidupdate, when applied. As basic length we use the meanparticle distance in the crystal phase a = 1 .
1. The timestep is dτ ∗ = 0 . τ , with τ = a p k B T /m p , where m p isthe mass of the colloids. The viscosity is η = 0 . ρ fl = 1 . ρ = 1 .
0. The friction γ varies between 0 . . ,
384 particles in a D T l D w/o HIsw HIs FIG. 1: The diffusion coefficient of the tracer particles in thebulk D Tl as a function of the single particle diffusion coeffi-cients D . (in simulation units). Red squares show results forthe system without HIs and the blue triangles for the systemswith HI. The gray dashed lines are a guide to eye. The blacklines illustrate the matching of the diffusion coefficient. Twodifferent single particle diffusion D in the Langevin dynam-ics and the LB coupling result in the same long-time diffusioncoefficient D Tl . box of size 64 × ×
16 without confining walls. Figure 1shows the measured long-time diffusion coefficient D Tl oftracer particles in the bulk calculated from the MSD, asa function of the single particle diffusion coefficient D .The simulations with and without HIs led to significantlydifferent tracer diffusion coefficients in the bulk, due tothe fact that e.g. the single particle diffusion in case ofLB is not simply D ∝ /γ [2]. Furthermore, there isno simple analytical rule to calculate D Tl in case of asoft-sphere bulk systems. In order to match the tracerparticle diffusion coefficient D Tl , we used different val-ues of γ in the Langevin dynamics and the LB coupling.The black line in Fig. 1 illustrates this matching: in or-der to obtain the same diffusion coefficient D Tl = 0 . γ = 4 . γ = 7 . γ are always the onesthat apply to the LB coupling.Using the matched tracer diffusion, we investigatedthe freezing of the undercooled fluid confined betweentwo planar walls. In order to distinguish the liquid andthe different solid phases, we used the Steinhardt or-der parameter [30]. Investigations by Moroni et al. [22]showed that especially q and q are good choices todetermine whether cubic or hexagonal structures arepresent in the system, respectively. In the following,we will focus on FCC and BCC crystal structures, forwhich the q order parameter is well suited. Due tothe strong fluctuations in our system, we applied an en-hanced averaging method for the Steinhardt order pa-rameter q , introduced by Lecher [16]. The literaturevalues are q (BCC) = 0 . q (HCP) = 0 . q (LIQ) = 0 . q of three snapshotstaken at different times during a typical simulation run.Note that the points only represent the peaks the distri-bution of q , since in the crystal, there is a strong layeringparallel to the wall. In between the peaks, the densitydrops nearly to zero in the crystal, and consequently sodoes the order parameter. As expected, the crystal startswith a HCP wall layer, followed by a BCC crystal frontthat grows with time. To evaluate the position of thecrystal front, we fitted the q peaks to a function of shape − h · arctan(( x − s ) /w ), where x is the x -position in thesimulation box, h is the height difference between q inthe liquid and in the FCC phase, w is the width of theliquid-crystal transition region and s is the position ofthe crystal front. Note that q in the liquid bulk is largerthan the literature value, since we report only the peaks,not the usual average value. q - x * q (t )q (t )q (t ) FIG. 2: The figure shows our fitting procedure for s (positionof the crystal front), from the function of q depending on x ∗ = x/a for different times t < t < t . The gray dashedlines indicate the computed front locations ( s < s < s ),while the black dashed line shows the fit for s at time t .The symbols represent the peaks of the q order parameter. After some initial time, the crystal grew very uni-formly, so that we could determine a constant growthvelocity u by a linear fitting of the front position. Fig-ure 3 shows the measured velocities u as a function of thehydrodynamic radius r H , which we varied by changingthe friction coefficient γ and applying the matching pro-cedure described above. Every measurement representsthe mean growth velocity sampled from 24 independentruns. For hydrodynamic radii r ∗ H = r H /a < . u / ( D T l a ) r * H w/o HIsw HIs FIG. 3: Growth velocity u normalized by the bulk diffusioncoefficient D Tl times the mean particle distance a as a functionof the relative hydrodynamic radius r ∗ H = r H /a . The redsquares show the results for simulations without HIs, whichare almost independent of r ∗ H within the error bars. The bluetriangles represent the results for simulations with HIs, whichshow a strong decay of the growth velocity with increasinghydrodynamic radius. The gray dashed lines are guidelinesto eye. a is the mean particle distance in the crystal phase, theinfluence of HI is almost negligible as one would expect.But already in case of moderate ratios 0 . < r ∗ H < . r ∗ H = 0 .
25. In case of noHIs the normalized growth velocity is virtually constant,with a decay for small frictions due to improper couplingto the thermostat.In order to elucidate what causes this difference, weanalyzed the particle diffusion in the system relative tothe actual position of the crystal front. To accomplishthis, we binned our system along the growth directioninto bins of width b = 2 a and determined the long-timediffusion coefficient D COMx of the center of mass in thedirection of growth, which can be seen as a measure forthe transport of particles towards the growing crystalfront. In Fig. 4 D COMx relative to the position of thecrystal front x rel , normalized by the center of mass dif-fusion in the bulk D COMlv of the Langevin simulation isshown. The front of the crystal is located at x rel = 0,while the pure bulk fluid phase is located at x rel = 7 andthe crystal phase is present for x rel <
0. As expected,the long-time diffusion coefficients for the center of massin the crystalline region are almost zero, rise in the re-gion of the crystal front, and settle off to the liquid bulkvalue far away of the crystal front. The left-hand side ofFig. 4 shows the values for low r ∗ H = r H /a = 0 .
025 ratio,which are virtually the same for both systems: with and D x C O M / D l v C O M x rel (a) w/o HIsw HIs -6 -4 -2 0 2 4 6( )(b) FIG. 4: Long-time diffusion coefficient of the center of massin the direction of growth. Every value is the long-time diffu-sion coefficient of the center of mass for particles in a bin ofsize b = 2 a normalized by the value of the Langevin simula-tion (w/o HIs). The front of the crystal is located at x rel = 0,while x rel = 7 is located in the pure bulk fluid phase. (a) Thefirst two cases (I: red circles, II: blue rhombus) have a low r ∗ H = r H /a = 0 .
025 ratio and show virtually similar behavior.(b) For moderate r ∗ H = 0 .
25, case IV (blue triangles) with HIsshows a different behavior from case III (red squares) withoutHIs, especially in the region in front of the crystal phase. without HIs. However, on the right-hand side r ∗ H = 0 . et al. [9] and von Hansen etal. [32].Our simulations show that hydrodynamic interactionshave a strong influence on crystallization, even at mod-erate hydrodynamic radii. Similar finding has been re-ported by Schilling et al. [27] as well. The effects arisemainly on the particle transport towards the crystalfront, which are in particular important for nucleationprocesses. This puts the common assumption into doubtthat hydrodynamic interactions can be ignored whenstudying crystallization or nucleation in suspensions. Atleast in Yukawa suspensions, these processes do not seemto be quasi-static, and the often drawn analogy to truemelts might not be true.We thank the staff at the Institute for ComputationalPhysics for useful discussions and acknowledge financial support from the German Science Foundation (DFG)through SFB716 and the cluster of excellence SimTech. ∗ Electronic address: [email protected] † Electronic address: [email protected] ‡ Electronic address: [email protected][1] R. Adhikari, K. Stratford, M.E. Cates, and A.J. Wag-ner. Fluctuating Lattice Boltzmann.
Euro. Lett. , 71:473,2005.[2] P. Ahlrichs and B. D¨unweg. Simulation of a single poly-mer chain in solution by combining lattice Boltzmannand molecular dynamics.
J. Chem. Phys. , 111(17):8225–8239, 1999.[3] H.C. Andersen. Molecular dynamics simulations at con-stant pressure and/or temperature.
J. Chem. Phys. ,72:2384, 1980.[4] A. Arnold, O. Lenz, S. Kesselheim, R. Weeber,F. Fahrenberger, D. Roehm, P. Koˇsovan, and C. Holm.ESPResSo 3.1 — Molecular Dynamics Software forCoarse-Grained Models. In M. Griebel and M. A.Schweitzer, editors,
Meshfree Methods for Partial Dif-ferential Equations VI , volume 89 of
Lecture Notes inComputational Science and Engineering , pages 1–23.Springer, 2013.[5] A. Arnold, B.A. Mann, H. Limbach, and C. Holm.ESPResSo – An Extensible Simulation Package for Re-search on Soft Matter Systems. In Kurt Kremer andVolker Macho, editors,
Forschung und wissenschaftlichesRechnen 2003 , volume 63 of
GWDG-Bericht , pages 43–59. Gesellschaft f¨ur wissenschaftliche Datenverarbeitungmbh, G¨ottingen, Germany, 2004.[6] B. D¨unweg and A.J.C. Ladd. Lattice boltzmann simu-lations of soft matter systems. In
Advanced ComputerSimulation Approaches for Soft Matter Sciences III , vol-ume 221 of
Advances in Polymer Science , pages 89–166.Springer, Berlin, Germany, 2009.[7] A. Engelbrecht, R. Meneses, and H.J. Sch¨ope. Heteroge-neous and homogeneous crystal nucleation in a colloidalmodel system of charged spheres at low metastabilities.
Soft Matter , 7(12):5685–5690, 2011.[8] P. Espa˜nol and P. Warren. Statistical mechanics of dissi-pative particle dynamics.
Europhys. Lett. , 30:191, 1995.[9] M.I.M. Feitosa and O.N. Mesquita. Wall-drag effect ondiffusion of colloidal particles near surfaces: a photoncorrelation study.
Phys. Rev. A , 44(10):6677, 1991.[10] S. Hamaguchi, R.T. Farouki, and D.H.E. Dubin. Phasediagram of yukawa systems near the one-component-plasma limit revisited.
J. Chem. Phys. , 105:7641, 1996.[11] S. Hamaguchi, R.T. Farouki, and D.H.E. Dubin. Triplepoint of yukawa systems.
Phys. Rev. E , 56(4):4671, 1997.[12] E.J. Hinch. Application of the langevin equation to fluidsuspensions.
J. Fluid Mech. , 72:499–511, 1975.[13] P.J. Hoogerbrugge and J.M.V.A. Koelman. Simulatingmicroscopic hydrodynamic phenomena with dissipativeparticle dynamics.
Euro. Lett. , 19(3):155–160, 1992.[14] V. Kalikmanov. Classical nucleation theory.
Nuc. The-ory , pages 17–41, 2013.[15] R. Klein and P. Meakin. Universality in colloid aggrega-tion.
Nature , 339, 1989.[16] W. Lechner and C. Dellago. Accurate determination of crystal structures based on averaged local bond order pa-rameters. arXiv preprint arXiv:0806.3345 , 2008.[17] H. L¨owen. Dynamical criterion for two-dimensional freez-ing.
Phys. Rev. E , 53(1):29–32, 1996.[18] H. L¨owen. Possibilities of phase separation in colloidalsuspensions.
Physica A , 235:129–141, 1997.[19] H. L¨owen, T. Palberg, and R. Simon. Dynamical cri-terion for freezing of colloidal liquids.
Phys. Rev. Lett. ,70(10):1557–1560, 1993.[20] M.G. McPhie and G. N¨agele. Long-time self-diffusion ofcharged colloidal particles: Electrokinetic and hydrody-namic interaction effects.
J. Chem. Phys. , 127:034906,2007.[21] E.J. Meijer and D. Frenkel. Melting line of yukawa sys-tem by computer simulation.
J. Chem. Phys. , 94:2269,1991.[22] D. Moroni, P.R. Ten Wolde, and P.G. Bolhuis. Interplaybetween structure and size in a critical crystal nucleus.
Phys. Rev. Lett. , 94(23):235703, 2005.[23] G. Naegele. On the dynamics and structure of charge-stabilized suspensions.
Phys. Rep. , 272(5–6):215–372,1996.[24] J. T. Padding and A. A. Louis. Hydrodynamic in-teractions and brownian forces in colloidal suspensions:Coarse-graining over time and length scales.
Phys. Rev.E , 74(3):031402, 2006.[25] R. Pesch´e, M. Kollmann, and G. N¨agele. Brownian dy-namics study of dynamic scaling and related freezingcriteria in quasi-two-dimensional dispersions.
J. Chem. Phys. , 114:8701, 2001.[26] D.N. Petsev and N.D. Denkov. Diffusion of charged col-loidal particles at low volume fraction: theoretical modeland light scattering experiments.
J. Coll. Interf. Sci. ,149(2):329–344, 1992.[27] M. Radu and T. Schilling. Solvent hydrodynamics af-fect crystal nucleation in suspensions of colloidal hard-spheres. arXiv preprint arXiv:1301.5592 , 2013.[28] D. Roehm and A. Arnold. Lattice Boltzmann simulationson GPUs with ESPResSo.
EPJ -ST , 210:89–100, 2012.[29] U.D. Schiller.
Thermal fluctuations and boundary condi-tions in the lattice Boltzmann method . PhD thesis, Jo-hannes Gutenberg-Universit¨at Mainz, 2008.[30] P.J. Steinhardt, D.R. Nelson, and M. Ronchetti. Bond-orientational order in liquids and glasses.
Phys. Rev. B ,28(2):784, 1983.[31] O.B. Usta, A. J.C. Ladd, and J.E. Butler. Lattice-boltzmann simulations of the dynamics of polymer so-lutions in periodic and confined geometries.
J. Chem.Phys. , 122:094902, 2005.[32] Y. von Hansen, M. Hinczewski, and R.R. Netz. Hy-drodynamic screening near planar boundaries: Effectson semiflexible polymer dynamics.
J. Chem. Phys. ,134(23):235102, 2011.[33] D.P. Ziegler. Boundary conditions for lattice boltz-mann simulations.