The role of the attractive forces in a supercooled liquid
aa r X i v : . [ c ond - m a t . s o f t ] F e b The role of the attractive forces in a supercooled liquid
Søren Toxvaerd
DNRF centre “Glass and Time,” IMFUFA,Department of Sciences, Roskilde University,Postbox 260, DK-4000 Roskilde, Denmark (Dated: February 8, 2021)
Abstract
Molecular Dynamics simulations of crystallization in a supercooled liquid of Lennard-Jones par-ticles with different range of attractions shows that the inclusion of the attractive forces fromthe first, second and third coordination shell increases the trend to crystallize systematic. Thebond-order Q in the supercooled liquid is heterogeneously distributed with clusters of particleswith relative high bond-order for a supercooled liquid, and a systematic increase of the extent ofheterogeneity with increasing range of attractions. The onset of crystallization appears in such acluster, which together explains the attractive forces influence on crystallization. The mean squaredisplacement and self-diffusion constant exhibit the same dependence on the range of attractions inthe dynamics and shows, that the attractive forces and the range of the forces plays an importantrole for bond-ordering, diffusion and for crystallization. . INTRODUCTION Ever since Berni Alder [1] in 1957 performed the first Molecular Dynamics (MD) simu-lation of a hard sphere system with crystallization there has been a general understandingof, that crystallization of a simple liquid is given by the harsh repulsive short range forces,and the MD simulation with particles with the more realistic Lennard-Jones (LJ) potentialsupported the assumption [2, 3]. Not only did the thermodynamics of a LJ system agreedwith the corresponding behaviour of a system of noble gas atoms, but the tendency to crys-tallize for LJ systems with and without attractive forces is also very similar (Figure 1). Thissimilarity is explained with, that there is an overall agreement between the radial distribu-tion function g ( r ) for LJ systems with- and without the attractive forces, as well as with theradial distribution function for a hard-sphere system. These rather closely similarities in theradial distributions of the particles have given reason to the well-established “perturbationtheory” [4, 5], where the thermodynamic and dynamic behavior of a system is obtainedfrom systems of purely repulsive particles by mean field corrections for contributions fromthe attractive forces.Here we analyse the role of the attractive forces on the supercooled state and the crys-tallization by Molecular Dynamics (MD) simulations of LJ systems with different range ofattractions. The simulations show, that the attractive forces play an important role in asupercooled liquid. They increase the bond-order in the supercooled liquid, given by Q [6, 7], and the tendency to crystallize. II. CRYSTALLIZATION OF THE SUPERCOOLED LIQUIDS
The systems with different range of attractions are cooled down from the liquids at thestate point ( ρ l , T c ) =(1.095, 2.25) for a liquid in equilibrium with fcc solid. The ( ρ, T ) phasediagrams for LJ systems with- and without attractive forces are shown in Figure 1. The T ( ρ )curves for coexisting liquids and fcc solids are obtained thermodynamically [8, 9]. At ( ρ l , T c )=(1.095, 2.25) a liquid with only repulsive forces and the LJ liquid with (full) attractionscrystallize at the same state point.The MD systems consist of N = 80000 LJ particles, where the forces are “cutted and2 fluid fcc crystal T e m pe r a t u r e Density R ad i a l d i s t r i bu t i on g (r) Radial distance r
FIG. 1: The liquid-solid phase diagram for a LJ and a WCA system. LJ: green is liquid inequilibrium with gas; red: liquid in equilibrium with fcc solid: blue. WCA: black is fluid inequilibrium with fcc solid: light blue. The systems are cooled down from liquids at ( ρ, T ) =(1.095,2.25): black point to ( ρ, T ) =(1.095, 1.25): red point. Inset: the radial distribution functions g(r)at (1.095, 1.25). Red: supercooled liquid, and with black points for WCA. Blue: fcc crystal andwith blue points WCA. shifted” to zero at different particle distances greater than r c [10]. The simulations (unitlength l ∗ : σ ; unit time t ∗ : σ p m/ǫ , for computational details see [10]) are performed forfour different values of r c : 3.5, 2.3, 1.41 and 2 / , respectively. Only the strong repulsive LJforces are included in the dynamics for the short cut at r c = 2 / , and this system appearsin the literature with the name WCA [5]. For r c = 1 .
41 also the attractive forces fromparticles in the first coordination shell are included, for r c = 2 . r c = 3 . g ( r ) for the two systemsin the supercooled state (red point in the Figure). With red is g ( r ) for the supercooledliquid and the black points are g ( r ) for the WCA system. The blue curve in the inset is g ( r ) for a LJ fcc crystal and with blue points for fcc WCA. The g ( r ) for the WCA systemsare shown with points for illustrative reasons because the differences between g ( r ) for LJand WCA are small. The overall similarity between the two systems g ( r ) implies, that themean pair-distributions in the supercooled state with coordination shells around a particles3 P o t en t i a l ene r g y pe r pa r t i c l e Time
FIG. 2: Energy per particle as a function of log − time after the quench to the supercooled state.With black lines are for the systems with only repulsive forces (WCA). The six blue curves are for r c =1.41 (with the attractive forces in the first coordination shell included); green six curves are for r c =2.3 (with attractive forces also from second coordination shell); red six curves r c =3.50 (withattractive forces also from third coordination shell). are almost identical for the two systems, and that the mean effects of the attractive forceson pressure, energy and free energy can be obtained as mean field contributions.The latent heat is released when a supercooled liquid crystallizes spontaneously, and theenergy decreases and the temperature increases without a thermostat. In [11] the effect ofa thermostat on the spontaneous crystallizations in the big MD supercooled systems wasinvestigated by performing ensemble simulations with- and without a thermostat and withthe conclusion, that the intensive MD thermostat, as expected, had no effect on the onset ofcrystallization. The present ( N V T ) simulations are with a thermostat by which the latentheat is removed smoothly, and the energy per particle decreases during the spontaneouscrystallization at the constant supercooled temperature (Figure 2).The systems are cooled down from liquids at the state point where the liquids are inthermodynamic equilibrium with fcc solid, and at the point on the coexisting phase lineswhere the lines for the two systems crosses each other, by which the degree of supercooling
T /T c = 1 . / .
25 = 0 .
556 at the constant density ( ρ l = 1 . log − time evolution of the energies per particle for systems in thesupercooled state and for different values of r c . The mean field energies from particles in4 P r obab ili t y d i s t r i bu t i on o f Q Q FIG. 3: Probability distribution P( Q ) of bond-order Q in the supercooled state. Red: LJsupercooled liquid with r c = 3 .
5; black: supercooled WCA fluid; blue: fcc LJ with r c = 3 .
5; lightblue: WCA fcc. Inset: log (P) in the interval Q ∈ [0 . , . the interval [ r c , .
5] are added to the functions u ( t ) = u ( t, r c ) + 2 πρ Z . r c g ( r ) u LJ ( r ) r dr, (1)and the energies u ( t ) for different cuts of the forces are almost equal before the onset ofcrystallization in accordance with the perturbation theory. The time evolution are shownwith a logarithmic time scale. With black lines are the WCA systems with only repulsiveforces, and they were simulated ∆ t = 22000 (2.2 × time steps). However, they remainedin the supercooled state without crystallization. (The WCA systems were crystallized ata lower temperature T =1.15.) The six blue curves are for r c =1.41 with the attractiveforces in the first coordination shell included in the force, and they crystallized within thetime interval [1200 , r c =2.3 withattractive forces also from the second coordination shell, and they crystallized within thetime interval [150 , r c =3.50 with attractive forces also fromthe third coordination shell, and they crystallized within the time interval [80 , . < r c on the stability of the supercooled statewas, however, not investigated due to a lack of computer facilities for these very demandingsimulations.(The energies after the crystallization are rather different. In general a hard sphere systemas well as a LJ system crystallizes with polymorphism [15–17] to different polycrystallinefcc states with different mean energy per particles [11], as also seen in Figure 2.)5 M ean nu m be r o f c l u s t e r s w i t h N pa r t i c l e s N particles C l u s t e r s w i t h N N particles
FIG. 4: Distributions (logarithmic) of clusters of N particles with Q > .
25 in the supercooledLJ liquids. Red: for the LJ system with r c =3.5; green: r c =2.3; blue: r c =1.41; black: r c = 2 / (WCA). Inset: The discrete distributions of big clusters (with points). The sensitivity of the range of attractions to the ability to crystallize is surprising giventhat the pair distributions of the different systems are very similar in the supercooled stateas well as in the crystalline state. In the classical nucleation theory the size of the criticalnucleus is the size, where the gain in free energy by an increase of particles in the crystalnucleus equals the cost of the increasing surface free energy, and these excess free energiesshould not be sensitive to the range of attractions due to the similarities in g ( r ). But thedistribution of bond-order Q for the particles is sensitive to the range of the attractiveforces.A supercooled LJ liquid is characterized by a heterogeneous distribution of bond-order,given by Q [11]. Here we will argue, that it is the attractive forces impact on the extent of theheterogeneity of the bond-order, which causes the difference in the tendency to crystallize.The distribution P ( Q ) of bond order Q ( i ) for the particles i in the LJ supercooled state isshown in Figure 3 and with an inset, which shows the log − distributions in the bond-orderinterval for which the particles with these bond-order are heterogeneous distributed. As in[11] there is an overlap in the distributions in supercooled liquid and in fcc crystal in this P ( Q ) interval Q ∈ [0 . , . Q > .
25 were heterogeneous distributed and with some particles withbond-order Q > .
35. Furthermore the critical crystal nucleus appeared in such a domainand with mean bond-order < Q > ≈ .
38, which is significantly less than the bond-order inthe fcc crystal at the same state point. 6 N pa r t i c l e s Time Q Time
FIG. 5: An example of the number N ( t ) of particles in the biggest cluster of particles i with Q ( i ) > .
25. The figure shows N ( t ) at the onset of crystallization for a system with r c = 2 . N c =0.67 and < Q > =0.3938. The number of clusters with N particles with a bond-order 0 . < Q are shown in thenext figure. (A particle i in a cluster with Q ( i ) > .
25 is close to ( r ij < .
41) at least oneother particle j in the cluster.) The distributions are obtained for the supercooled state asthe mean of 200 independent determinations, and the inset shows the number of discrete andrare events of bigger clusters. The figure shows two things. For the first there is a crucialdifference between the purely repulsive force system (WCA) and the systems with attractiveforces which all exhibit big clusters with high liquid bond-order. And secondly, the insetshows, that although the three distributions with different range of attractions looks prettysimilar, there appears occasionally a much bigger clusters for the systems with long rangeattractions.The critical nucleus were determined as described in [11]. Figure 5 shows a representativeexample of the time evolution of the number of particles N ( t ) in the biggest cluster, andwith the mean bond-order in the inset of the figure. The estimated critical sizes < N c > with the bond-order < Q > for the simulations are: r c = 1 .
41 : < N c > = 73 ± < Q > = 0 . ± . r c = 2 .
30 : < N c > = 73 ± < Q > = 0 . ± . r c = 3 .
50 : < N c > = 73 ± < Q > = 0 . ± . < Q > = 0 .
390 in the critical nucleus is much less than in an7rdered fcc crystal and the size of critical nuclei is the same for the three ranges of attractions.The extension of domains with relative high bond-order varies, however, with the range ofthe cut (Figure 4), and the critical nuclei appear in a domain with high bond-order. So theextend of the heterogeneous distribution of high bond-order in the supercooled liquid andthereby the probability to obtain a critical nucleus can explain the observed differences inthe tendency to crystallize.The existence of “dynamic heterogeneity” in supercooled liquids have been known fora long time [12–14], and in [11] is was linked to the bond-ordering in subdomains. If sothe viscosity and particle diffusion should be different for domains with relative low bond-order compared with subdomains with relative high bond-order. This is, however, difficultto determine directly because the domains are not permanent and particles change bond-order with time. But the overall particle diffusion reveal the differences and the main effectof the attractive forces on diffusion and viscosity. The next figure gives the mean squaredisplacements of a particle in the supercooled state and for different range of attractions. Thefigure shows that there is a difference in the slopes end thereby the self-diffusion constants D for the different range of attractions. The self-diffusion constants are WCA: D =0.01315; r c = 1 . D =0.01033; r c = 2 . D =0.00991; r c = 3 . D =0.01004. The inset is themean square displacements in logarithmic scales, and it shows that the short time “ballisticregime” is similar for all four systems and given by the strong repulsive forces. The behaviourof the particle diffusion with respect to the range of the attractions can be explained by,that the domains with relative high bond-order slow down the mobility of a particle in thesedomains, i.e. that the mobility is heigh in domains where the bond-order is small and smallin domains with relative high bond-order, where the particles are tied together weakly.The attractive forces dynamic effect in supercooled states has also been obtained formixtures [18, 19] III. CONCLUSION
In conclusion
The radial distribution functions g ( r ) for the supercooled liquid as well asfor the fcc solid are insensitive to the range of the attractions, and hence the free energiesper particle (chemical potential) are also insensitive. Consistent with this observation, so isthe size of the critical nuclei; but nonetheless, the tendency to crystallize depends on the8 M ean s qua r e d i s p l a c e m en t time FIG. 6: Mean square displacements at (
T, ρ ) = (1.25, 1.095). Black: for r c = 2 / (WCA); blue: r c =1.41; green: r c =2.3; and red: r c =3.5, respectively. Inset: In logarithmic scales. range of attraction. The systems in the supercooled state exhibit, however, heterogeneousdistributed particles with a relative high bond-order for the supercooled state, and the extentof the heterogeneity is enhanced mainly from the attractions from the particles within thefirst coordination shell, but also the particles from the second- and third coordination shellsincrease the number of domains with relative high bond-order. In accordance with thisobservation, the systems crystallize much more easily for the systems with attractions andin a systematic way so that six crystallizations of particles with attractions from particleswithin three coordination shells crystallized ≈ eight times faster than six systems withattractions only from the first coordination shell, whereas the systems with only repulsiveforces did not crystallized at the supercooled state point.This behavior of the heterogeneous bond-order distribution is consistent with the wellknown “dynamic heterogeneity” in supercooled liquids, and the self-diffusion for the particleswith different range of attractions supports the hypothesis. The attractions slow down theself-diffusion, the main effect comes from the attractions within the first coordination shell,but also the longer-range attractions affect the diffusion. So in summary the attractiveforces enhance the extent of the domains with high bond-order, slow down the particlediffusion and catalyze the crystallization. The sensitivity of the crystallization to the rangeof attractions makes it difficult to compare nucleation rates obtained by simulations withexperimentally determined nucleation rates.9his work was supported by the VILLUM Foundation’s Matter project grant No. 16515. [1] B. J. Alder and T. E. Wainwright, J. Chem. Phys. , 1208 (1957).[2] L. Verlet Phys. Rev. , 98 (1967).[3] J. P Hansen and L. Verlet, Phys. Rev. , 151 (1969).[4] J. A. Barker and D. Henderson, J. Chem. Phys. , 4714 (1967).[5] J. Weeks, D. Chandler and H. C. Andersen, J. Chem. Phys. , 5237 (1971).[6] P. J. Steinhardt, D. R. Nelson and M. Ronchetti, Phys. Rev. B, , 784 (1983).[7] W. Lechner and C. Dellago, J. Chem. Phys. , 114707 (2008).[8] A. Ahmed and R. J. Sadus, Phys. Rev E , 061101 (2009).[9] M. A. Barroso and A. L. Ferreira, J. Chem. Phys. , 7145 (2002).[10] S. Toxvaerd and J. C. Dyre, J. Chem. Phys. , 081102 (2011).[11] S. Toxvaerd, Eur. Phys. J. B , 202 (2020).[12] M. D. Ediger, Annu. Rev. Phys. Chem. , 99 (2000).[13] A. Widmer-Cooper, P. Harrowell and H. Fynewever, Phys. Rev. Lett. , 135701 (2004).[14] A. Widmer-Cooper and P. Harrowell, Phys. Rev. Lett. , 185701 (2006).[15] B. O’Malley and I. Snook, Phys. Rev. Lett. , 085702 (2003).[16] C. Desgranges and J. Delhommelle, Phys. Rev. Lett. , 235502 (2007).[17] C. Desgranges and J. Delhommelle, J. Phys. Chem. B , 1465 (2007).[18] L. Berthier and G. Tarjus, J. Chem. Phys. , 214503 (2011).[19] S. Toxvaerd and J. C. Dyre, J. Chem. Phys. , 134501 (2011)., 134501 (2011).