Crystallization of anisotropic colloids with a Yukawa potential
CCrystallization of anisotropic colloids with a Yukawa potential
Fu-Jun Lin , Jing-jing Liao , and Bao-quan Ai ∗ Guangdong Provincial Key Laboratory of Quantum Engineering and Quantum Materials,Guangdong-Hong Kong Joint Laboratory of Quantum Matter,School of Physics and Telecommunication Engineering, South China Normal University,Guangzhou 510006, China. School of Science, Jiangxi University of Science and Technology, Ganzhou 341000, China.
Abstract
Crystallization in a dense suspension of anisotropic spherical colloidal particles with aYukawa potential is numerically investigated in a two-dimensional plane. It is foundthat a strong anisotropy can hinder the particles from crystallizing, while a weakanisotropy but super-strong coupling facilitates colloids to freeze into a hexagonalcrystal. Different criterions are employed to describe the phase transition, one can findthat a competition between anisotropic degree and coupling strength shall widened thetransition region in the phase diagram, where the heterogeneous structures coexist,which render as a quasi-platform stretched across the probability distribution curve ofthe local order parameter. Our study maybe helpful for the experiments relating to thecrystallizing behavior in statistical physics, materials science and biophysical systems. ∗ Email: [email protected] . INTRODUCTION
Phase transition is a universal phenomenon in nature as well as in industry, and playsan essential role in statistical physics, materials science, chemistry and biophysics.Crystallization represents the prime example of a phase transition and hence gets a lotof attention. Understanding the crystallization kinetics and microscale processes arenot only important for fundamental research, but also for the applications associatedwith crystalline materials [1] and biotechnology [2,3].Quantitative crystallization examinations are mostly based on the models ofrepulsive hard spheres (discs) [4-8], charged spheres [9-11], or attractive spheres[12,13]. As far as we know, it is challenging experimentally to investigate the particlesdynamics on an atomic/molecular scale, due to the small size. Therefore, the colloidalparticles, are usually viewed as the large atoms with tailorable size, shape andinteractions [14], have been extensively studied to explore the physics ofcrystallization [15-20], because their thermal motion can be visualized by opticalmicroscopy [21] and tracked by image processing [22]. The systematic studies cansignificantly enhance understanding and help controlling of the crystallizationprocesses. A flowing crystalline structure, named “ rheocrystal ” , spontaneously formsby removing the particles in experimental work, supported by numerical simulations[23]. By doping colloidal glasses with small amounts of self-propelled particles, thecrystallization dynamics can be speeded up [6]. Additionally, the molecular dynamicssimulation of equal-sized hard-spheres systems, showing conventional nucleation andgrowth of crystals at concentrations near melting and crossing over to a spinodal-likeregime, and a Monte Carlo ‘constrained aging’ method slows crystallization [24]. It isworth to mentioning that an excellent work on crystallization of monodisperseself-propelled colloidal particles at sufficiently high densities has been done by Bialket.al. [15], which provides us the inspiration for doing the current work. Furthermore,the crystallization in system of binary hard spheres with different sizes have beenreported in experiments [25,26] and theory [27,28], the results indicate that the smallspheres can fit between large spheres to stabilize binary crystals. Bommineni et.al. [29]firstly perform the spontaneous formation of AlB , NaZn and Laves phases withswap simulations, they point out that the binary hard sphere crystals grow robustlyand reproducibly in standard event-driven-molecular dynamics, and the particles areintegrated faster into the crystal due to the enhanced diffusion.Most of the works on particle crystallization mentioned above have focused onspherical or isotropic particles. However, few natural particles have the perfectsymmetry. Anisotropy is an intrinsic property, including the shape (e.g., ellipsoids,rods, or cubes), interaction (e.g., janus) or diffusion (e.g., polar particles) anisotropy.Thus relevant researches are necessary and valuable. The needlelike ellipsoidalparticles under the action of an external potential exhibit complex motion and destroythe directed transport [30,31]. Triblock Janus colloids, are viewed as potentialuilding blocks for functional materials, their inherent anisotropy could induceself-assembly at a certain condition, have been used to measured the crystal growthkinetics [32]. A variety of crystalline structures have been discovered by varying theshape anisotropy of spheroids [33]. Predictably, simulations [34] of a generalizedVicsek model, report a mutual exclusion of the polar and the structural order. We aremotivated by these works to understand how the ability of particles to form crystallinestructures is affected by the anisotropy. In the paper, we consider a spherical colloidalparticle with different mobilities in orthogonality direction, and the interactionbetween particles is determined by the repulsive Yukawa potential. II. MODEL AND METHODS
FIG. 1. Schematic diagram illustrating the setup forming the basis of our calculations. (a) Modelof colloids moving in a two-dimensional box. (b) Representation of an anisotropic particle. Theaxes x and y are body-fixed coordinates, and the axes x and y are the lab-frame coordinates. In this section, we mainly review the derivation of the equations for Brownianmotion of an anisotropic spherical particle in two-dimension. Due to the particleanisotropy, the rotational and translational motion are always coupling in thelab-frame coordinates, which makes the analysis of relevant issues very complicated.For simplicity, we shall initially describe the motion of particle in body-fixedcoordinates where the rotational and translational motion are decoupled. Then, theparticle’s position vector S (t) of its center of mass at a given time t can bedecomposed in the body frame as ( x , y ), corresponding to the coordinates ( x , y ) in the lab frame. θ ( t ) is the angle between the x axis of the lab frame and the x axis of the body frame. Thus, it can be easily obtained from FIG.1(b) that theLangiven equations of overdamped motion of the anisotropic particle in the bodyframe are [35,36] x y x F t F t tt , (1) y x y F t F t tt , (2) t tt . (3)Where , are the mobilities along x and y axis, respectively. isrotational mobility. The noise i in the body frame has mean zero and correlations ,
2( ) ( ) , , 1, 2, 3,
Bi j i ji k Tt t t t i j (4)with T is temperature and B k is the Boltzmann constant. x F and y F are thecomponents of the force along the x and y axis of lab frame. For convenience, theenergy, length and time are normalized by B k T , and D , respectively.Here is the number density and D is the bare diffusion coefficient. For amany-particle system, particles are considered to interact pairwise through therepulsive Yukawa potential: r eu r r , with screening length , and thedimensionless coupling strength / B V k T , with the bare potential strength V .Thus, the force acts on the i th particle is i i i j i j u F r r . We now translate theEqs.(1)-(3) into that in lab-frame coordinates by means of a straightforward rotationof coordinates, the displacement relationship between the two frames is represented asfollowing: cos sin x x y , (5) sin cos y x y . (6)By a series of manipulations, the Langiven equations in the lab frame of overdampedmotion of the anisotropic particle can be obtained as following: cos 2 sin 2 ( ), x y x t F t F t tt (7) ( ) cos 2 sin 2 , y x y t F t F t tt (8) t tt (9)Where is the average of the mobility, and determines the anisotropy of the particle, indicates that the particle isisotropic, and while , the particle is prone to diffuse along x axis. Thenoise ( ) i t in the lab frame has mean zero and satisfies [35,36]: ( ) ( ) 2 , t t D t t (10) , ( ) ( ) 2 , , 1, 2, i j B ij t t k T t t i j (11)with cos 2 sin 2sin 2 cos 2 ij ij , (12)where B D k T is the rotational diffusion coefficient, describing the angularfluctuation. The superscripts of statistical averages in Eq.(11) mean over which thenoises are averaged, and the subscript quantity is kept fixed.In this paper, we focus on the influence of the particles’ anisotropic degree on thecrystallization behavior. Defining / 0,1 to describe the anisotropy, Thediffusion coefficient of the particle are B D k T and B D k T D respectively. We model a system with N spherical particles of radius r moving in a 2D box of size x y L L with periodicboundary conditions, and numerically investigate the dynamic behavior of thecrystallization by integration of the Langevin equations (7)-(9) using the second-orderstochastic Runge-kutta algorithm. III. NUMERICAL RESULTS AND DISCUSSION
In our simulations, we choose the same parameters as in
Ref .[15]: the totalparticle number N =1936, the particle radius r the rotational diffusioncoefficient D and the translational diffusion coefficient D the inversescreening length / 2 / 3 x y L L is chosen such that the particles cancrystallize into the perfect hexagonal crystal. The initial position and orientation ofparticles are random. The integration time step dt is chosen to be smaller than and the total integration time as great than , which ensure that the system reachesa steady state and the simulation results will not depend on the integration time andtime step.The structural transitions between the different phases can be characterized bythe global bond-orientational order parameter [37] ij N ii j i q i q eN , (13)where i is the set of the six nearest neighbors of the i th particle and ij is anglebetween the vector from particle i to j and an predefined direction. For a perfectcrystal whereas in a disordered phase Following Schweigert et al.[38], the liquid-to-solid transition can be identified by a jump of the order parameterbove a value of (see also Ref. [15]). Meanwhile, a dynamical criterionfor crystalline process of the anisotropic particles is determined by the abrupt drop ofthe long-time diffusion coefficient lim it D t r (14)with ( ) ( ) (0) i i i t t r r r . The value of 0.086 for D has been proved by Löwen (1996)to be ‘universal’, which does not matter whether freezing occurs continuously via ahexagonal phase or is a conventional first-order transition[39], thus we employ thiscriterion to describe the dynamical phase transition. FIG. 2. Cooling curves for (a) the bond-orientational parameter and (b) the long-timediffusion coefficient D versus particles’ anisotropic degree for different . The crossingswith the dashed horizontal lines define the position of the structural transition S ( )and the dynamical freezing D ( D =0.086), respectively.FIG. 3. Cooling curves for (a) the bond-orientational parameter and (b) the long-timediffusion coefficient D versus coupling strength at . We firstly examine the effect of anisotropic degree of particles on the phasetransition by monitoring the global bond-orientational order parameter and thelong-time diffusion coefficient D for different in FIG.2. At the beginning of theimulation, a particle configuration with random particle positions and orientations isemployed. From FIG.2(a) we find that for a given , the particles can freeze into anordered crystalline phase in the small region, especially for isotropic particles( = 0) , has its maximum value. With increasing , the Brownian diffusionalong a specific axis of a particle becomes increasingly dominant, the random motionhinders system from assembling an ordered structure. Therefore, the order parameterdrops to below = 0.45 at S , indicating a loss of the long-rangeorientational order in the system, and the system transiting into the liquid phase[40].Figure 2(b) exhibits dynamical phase transition for different , the diffusioncoefficient increases abruptly near the phase transition point D which gives thelower bound to liquid region, and even exceeds that of a free passive Brownianparticle. Obviously, the phase transition point is significantly shifted to the large with increasing, which is determined by the competition between the particles’ FIG. 4. Snapshots of particle configurations for (left column), 0.5 (middle column) and0.9 (right column), the rows correspond to constant , from top to bottom = 200(0.2886,0.2806, 0.2014), 700(0.7436, 0.6160, 0.2636) and 1400(0.8627, 0.8128, 0.2982). Particles arecolored according to their q . nisotropic degree and the coupling strength. As shown in FIG.3, with increasing,the corresponding temperature of the system is decreased, which leads to the densecolloids gradually freeze into an ordered configuration and the diffusion coefficientdrops toward zero.To effectively distinguish liquidlike from ordered regions, we utilize an orderparameter j i q q i q j per particle to describe the local behavior ofits neighbours [41]. The snapshots of particle configurations for , 0.5 and 0.9with different are illustrated in FIG.4. Clearly, for either small coupling strength(e.g. ) or large anisotropic degree (e.g. ), the system always stayliquid with small order parameter. Whereas for the isotropic particles ( ), thesystem crystallize into hexagonal structures with under strong coupling. FIG. 5. (a) Time dependence of the Lindermann parameter for = 1000. (b) Phase diagram inthe plane. To assess in more detail the phase transition, we consider a process of meltingstarting from a perfect hexagonal crystal, introduce the Lindemann-like parameter[42] ( ) ( )( ) 2 i jL t tt r r , (15)as a criterion to decide the upper bound to solid region. Where the subscripts i and j donate two particles that are initially neighbors, and the lattice spacing of thehexagonal crystal This criterion states that the meltingcommences once the vibrational displacement of a particle exceeds a certain fractionof the lattice spacing [15]. In our case, both the small coupling strength and largeparticles anisotropic degree can induce particles to vibrate sharply, the former hasbeen explicited in
Ref .[15]. We now mainly focus on the time dependence of theLindemann-like parameter for different at a given , e.g., = 1000 as shownin FIG.5(a). It is found that, in the liquid region, the curves of Lindemann-likeparameter are abruptly divergent over time. Especially for > 0.65, due to thexistence of intensive diffusion in a special direction, the particles easily escape fromtheir lattice position, and the crystal structure is destroyed. When = 0.6, one canfind a quasi-plateau with Lindemann-like parameter L , we defined this value L as melting point. Based on these above, the Phase diagram is mapped in FIG.5(b).Clearly, the curves of S and D mostly coincide in all parameter space, whichindicates that the two criterions we employed are compatible with our model. Notethat in the region with large and small , the particles crystalized into hexagonalstructure, whereas in the small or large region, the bond-orientational order isdestroyed due to long-time diffusion. Moreover, there is a transition regime betweenliquid and crystal, which is characterized by a high structural order and low butnonvanishing diffusion [15]. It is clear that the transition regime widens in the regionof moderate parameter space, which we regard as the result of the competitionbetween the particles’ anisotropic degree and coupling strength, as what shows inFIG.2. FIG. 6. Probability distributions of q (a) for different at = 0.5 and (b) at = 0 (solidline), 0.5(dashed line) and 0.9(dot dashed line) for four different global values. We begin by briefly reviewing our results to characterize the solid, liquid andtransition regions have been shown in FIG.5(b). Figure.6(a) illustrates the probabilitydistribution for q with different at = 0.5. One can find that the system is lessstructure in liquid region( = 200), while in the solid region ( = 1400), the particlesare well-crystallized with a few “bubbles” due to the long-range bond-orientationalorder. Notablely, For the transition region ( = 700), there is a broad peak probabilitydistribution for q , corresponding to a coexistence phase as the second snapshot ofmiddle column shown in FIG.4. Figure 6(b) displays the probability distributions of q at = 0, 0.5 and 0.9 for four different global value. As what has mentionedthat when = 0.9, the system is always disordered ( < 0.3), thus there is only onecurve. One can easily find that the probability distributions of q are less dependenton . It should be emphasized, a quasi-platform stretched across the distributioncurve of = 0.45, which means that the transition regime is occupied by a variety ofeterogeneous structures. IV. CONCLUDING REMARKS
To summarize, we have numerically investigated the crystallization in a densesuspension of spherical colloids with different mobilities in an orthogonal direction.Starting form a random particle configuration, we integrate the Langivin equationsusing the secondorder stochastic Runge-Kutta algorithm in a two dimensional boxwith periodic boundary condition. By employing different criterions, the phasediagram in the plane is obtained. It is found that the anisotropy hindersparticles from frozing into a crystalline structure. Especially for a large or small regime, the system always keeps liquid phase. However, the competition between and widened the transition region in the moderate parameter space, in whichthe suspension is overall ordered but with a lot of heterogeneous structures.Furthermore, our results indicate that the anisotropic colloidal spheres in a densesuspension can be frozen into perfect hexagonal crystals in the case of weakanisotropy and strong coupling. ACKNOWLEDGMENTS
This work was supported in part by the National Natural Science Foundation ofChina (Grant No.12075090, 11905086, 11575064), the GDUPS (2016), and theNatural Science Foundation of Guangdong Province (Grant No. 2017A030313029),and the Major Basic Research Project of Guangdong Province (Grant No.2017KZDXM024), and the Natural Science Foundation of Jiangxi Province (GrantNo. 20192BAB212006 and No. GJJ191598), and High-level Scientific ResearchFoundation for the introduction of talent of Jiangxi University of Science andTechnology.[1] D. Li, H. Zhou and I. Honma, Nat. Mater. , 65 (2004).[2] H. Liu, S. K. Kumar, J. F. Douglas, Phys. Rev. Lett. , 018101 (2009).[3] H. X. Lin, S. Lee, L. Sun, M. Spellings, M. Engel, S. C. Glotzer and C. A. Mirkin,Science , 931 (2017).[4] E. Zaccarelli, C. Valerians, E. Sanz, W. C. K. Poon, M. E. Cates, P.N. Pusey,Phys. Rev. Lett. , 135704 (2009).[5] E. Sanz, C. Valeriani, E. Zaccarelli, W. C. K. Poon, P.N. Pusey, M.E. Cates, Phys.Rev. Lett. , 215701 (2011).[6] R. Ni, M. A. Cohen Stuart, M. Dijkstra and P. G. Bolhuis. Soft Matt. , 6609(2014).7] G. Briand, O. Dauchot, Phys. Rev. Lett. , 098004 (2016).[8] A. Mori, Crystals, , 102 (2017).[9] S. Auer, D. Frenkel, J. Phys-Condens. Matter , 7667, (2002).[10] D. Reinke, H. Stark, H. H. von Gr¨unberg, A. B. Schofield, G. Maret and U.Gasser, Phys. Rev. Lett. , 074505 (2011).[12] V. J. Anderson, H. N. W. Lekkerkerker, Nature, , 811 (2002).[13] W. C. K. Poon, J. Phys. Condens. Matter , R859 (2002).[14] D. Frenkel, Science, , 65 (2002).[15] Z. Wang, A. M. Alsayed, A. G. Yodh and Y. Han, J. Chem. Phys. , 154501(2010).[16] J. Bialk´e, T. Speck and H. L¨owen, Phys. Rev. Lett. , 168301 (2012).[17] T. Palberg, J. Phys. Condens. Matter, ,7110 (2015).[19] B. Li, D. Zhou, Y. Han, Nat. Rev. Mater. , 15011 (2016).[20] D. P. Singh, U. Choudhury, P. Fischer and A. G. Mark, Adv. Mater. , 1701328(2017).[21] A. Kose, M. Ozaki, K. Takano, Y. Kobayashi, S. J. Hachisu, Colloid Interface Sci. , 330 (1973).[22] J. C. Crocker, D. G. Grier, J. Colloid Interface Sci. 179, 298 (1996).[23] G. Briand, M. Schindler and O. Dauchot, Phys. Rev. Lett. , 208001 (2018).[24] C. Valeriani, E. Sanz, E. Zaccarelli, W. C. K. Poon, M. E. Cates, P. N. Pusey, J.Phys. Condens. Matter , 194117 (2011).[25] A. B. Schofield, P. N. Pusey, and P. Radcliffe, Phys. Rev. E , 031407 (2005).[26] N. Schaertl, D. Botin, T. Palberg, and E. Bartsch, Soft Matter , 5130 (2018).[27] A. P. Hynninen, L. Filion, and M. Dijkstra, J. Chem. Phys. , 064902 (2009).[28] M. Dijkstra, Adv. Chem. Phys. , 35 (2014).[29] P. K. Bommineni, M. Klement, and M. Engel, Phys. Rev. Lett. , 218003(2020).[30] B. Q. Ai and J. C. Wu, J. Chem. Phys. , 062151 (2018).[32] W. F. Reinhart, A. Z. Panagiotopoulos, J. Chem. Phys. , 124506 (2018).[33] W. Jin, H. K. Chan, Z. Zhong, Phys. Rev. Lett. , 248002 (2020).[34] C. A. Weber, C. Bock, and E. Frey, Phys. Rev. Lett. , 168301 (2014).[35] Y. Han, A. M. Alsayed, M. Nobili, J. Zhang, T. C. Lubensky, and A. G. Yodh,Science, , 626 (2006); Y. Han, A. Alsayed, M. Nobili, and A. G. Yodh, Phys. Rev.E , 011403 (2009).[36] R. Grima and S. N. Yaliraki, J. Chem. Phys. , 084511 (2007); R. Grima, S. N.aliraki, and M. Barahona, J. Phys. Chem. B , 5380 (2010).[37] B. I. Halperin and D. R. Nelson, Phys. Rev. Lett. , 121(1978).[38] I. V. Schweigert, V. A. Schweigert, and F. M. Peeters, Phys. Rev. Lett. , 5293(1999).[39] H. Löwen, Phys. Rev. E , R29 (1996).[40] P. Hartmann, G. J. Kalman, Z. Donko and K. Kutasi, Phys. Rev. E , 026409(2005).[41] W. Lechner and C. Dellago, J. Chem. Phys. , 114707 (2008).[42] K. Zahn and G. Maret, Phys. Rev. Lett.85