15952
J. Phys. Chem. C 2007, 111, 15952-15956
Wall-Induced Prefreezing in Hard Spheres: A Thermodynamic Perspective† Brian B. Laird* Department of Chemistry, UniVersity of Kansas, Lawrence, Kansas 66045
Ruslan L. Davidchack Department of Mathematics, UniVersity of Leicester, Leicester, LE1 7RH, United Kingdom ReceiVed: May 15, 2007; In Final Form: August 21, 2007
Using our cleaving-wall method developed previously for the study of crystal-melt interfaces (J. Phys. Chem. B 2005, 109, 17802), we compute with high precision the interfacial free energies of the hard-sphere fluid (γwf) and crystal (γwc) in contact with a structureless hard wall. The interfacial free energies are determined as functions of density and, for γwc, as functions of crystal/wall orientation [(111), (100), and (110)]. Using the current calculations, along with our previously calculated value for the hard-sphere (111) crystal-fluid 111 111 interfacial free energy, γ111 cf , we show that, at the crystal-fluid coexistence conditions, γwf - [γwc + γcf ] > 0. This inequality demonstrates that a (111) hard-sphere crystal will completely wet a hard-sphere fluid/ structureless hard wall. This conclusion from thermodynamics is consistent with recent structural evidence by Dijkstra (Phys. ReV. Lett. 2004, 93, 108303) of wall-induced prefreezing of a hard-sphere fluid in contact with a structureless hard wall. We also find evidence that a (111) crystal will wet a (110) crystal/wall interface.
1. Introduction
cos θ )
Beginning with the computer simulations of Alder and Wainwright in 1957,1 the hard-sphere reference system has played an important role in the development of our current understanding of structure and phase equilibria in simple liquids and colloids.2 While the structure and bulk phase behavior of hard-spheres and hard-sphere mixtures are relatively wellunderstood, our understanding of the hard-sphere system in confined geometries is less complete. In 1992, Courtemanche and van Swol3,4 reported evidence from simulation that a hard-sphere fluid in contact with a planar, structureless hard wall develops an equilibrium “wetting” layer of fcc crystal with a (111) orientation, even at densities below the bulk freezing density. It was later shown by Kegel,5 however, that the simulations of Courtemanche and van Swol were insufficient, because of the small system size, to distinguish between true “prefreezing”, which is a wetting of a surface by a metastable crystal phase, and “capillary freezing”, that is, lowering of the freezing density because of confinement. This uncertainty was resolved recently by Dijkstra,6 who, in a series of simulations on hard-sphere fluids in a slit-pore geometry, presented compelling structural evidence for the existence of prefreezing of a hard-sphere fluid at a hard, structureless wall. Although prefreezing in the hard-sphere/hard-wall system has been demonstrated phenomenologically in simulation, it has not been possible to verify its existence thermodynamically because of a lack of precision in the calculation of wall-fluid and wallcrystal interfacial free energies.6,7 The degree of wetting of a surface by a material is described by the wetting contact angle, θ, defined through Young’s equation,8 which for the wetting of a solid at a wall-fluid boundary is †
Part of the “Keith E. Gubbins Festschrift”. * Corresponding author.
γwf - γwc γcf
(1)
where γwf, γwc and γcf are the wall-fluid, wall-crystal, and crystal-fluid interfacial free energies, respectively. Complete wetting corresponds to a vanishing contact angle, that is, when the right-hand side of eq 1 equals or exceeds unity. (In this case, γwf gγcf + γwc and a single wall-fluid interface is thermodynamically unstable relative to growth of a layer of crystal at the wall, generating two interfaces: wall-crystal and crystal-fluid.) Because the value of cos θ is proportional to the difference γwf - γwc, the wetting contact angle is extremely sensitive to errors in the wall-fluid and wall-crystal values. In this work, we calculate γwf and γwc for the hard-sphere/ smooth hard-wall system using a modified version of our cleaving-wall method, which was developed for the calculation of crystal-melt interfacial free energies.9,10 This method generates values of γwf and γwc at crystal-fluid coexistence densities that are sufficiently precise to allow a thermodynamic prediction of the complete wetting of a (111) crystal at the hard-sphere fluid-wall interface using eq 1. In recent years, two methods have been developed to accurately calculate γcf from simulation for a diverse array of systems from hard spheres9-11 to molecules, such as succinonitrile.12 The first method, developed by the authors,9 uses specially designed cleaving walls to reversibly separate bulk solid and fluid samples and reassemble the segments to form an equilibrium system containing a solid-liquid interface. The interfacial free energy is obtained from thermodynamic integration of the separation/reassembly process. In the second method,13 the interfacial stiffness is obtained from the power spectrum of the fluctuations in interfacial position. The interfacial free energy is obtained from an analysis of the orientation dependence of the interfacial stiffness. These two methods are complementary; the magnitude of γcf is determined more precisely using the “cleaving walls” approach, but the anisotropy
10.1021/jp073756u CCC: $37.00 © 2007 American Chemical Society Published on Web 10/11/2007
Wall-Induced Prefreezing in Hard Spheres
J. Phys. Chem. C, Vol. 111, No. 43, 2007 15953
(orientation dependence) of γcf is more precisely determined from capillary fluctuations. The two methods have been applied to the hard-sphere10,11 and Lennard-Jones14,15 systems, giving values of γcf that are in agreement within the error bars. For the hard-sphere system, the most precise current values for γcf were obtained using the cleaving-wall method yielding 0.557(7), 0.571(6), and 0.592(7) kT/σ2 for the (111), (110), and (100) crystal orientations, respectively, where σ is the hard-sphere diameter, k is Boltzmann’s constant, T is the absolute temperature, and the number in parenthesis is the error (2σ) in the last digit(s) shown. For γwf and γwc, the static nature of the wall precludes the use of the capillary fluctuation method to calculate these quantities, so either thermodynamic or mechanical integration techniques must be employed. In 1984, Henderson and van Swol used a molecular-dynamics simulation to calculate γwf for the hard-sphere fluid at a flat, structureless, hard wall from the mechanical definition of surface tension16
γ)
∫-∞∞ [pn(z) - pt(z)] dz
(2)
where pn and pt are the normal and transverse components of the pressure tensor and z is the direction normal to the wall. At a bulk fluid density of F* ) Fσ3 ) 0.901, the highest density they reported, they obtain γwf ) 1.8(6)kT/σ2, where the number in parenthesis is the error in the last digit shown. (Their reported value for γwf is actually negative because they did not include the contribution of the external field to the interfacial free energy,17 which is given by pσ/2, where p is the bulk pressure.) It should be noted that eq 2 cannot be used to determine γwc or γcf because the relaxation of stress in a crystal takes place on a much longer time scale than is accessible in simulation, and the interfacial excess stress measured with eq 2 will not equal the interfacial free energy.18 In a later work, Heni and Lo¨wen19 used thermodynamic integration with square-barrier and triangular cleaving potentials to calculate γwf and γwc. At coexistence, they obtain γwf ) 1.99(18)kT/σ2 and γwc ) 1.42(10), 2.0(3), and 3.1(3) for the (111), (100), and (110) crystal orientations, respectively. (Their value for γwf was obtained by extrapolating data at lower densities, as they were unable to simulate the fluid at the coexistence density without precrystallization.) Using these values of γwf and γwc, together with the γcf values from ref 10, we find Young’s equation (eq 1) yields values for cos θ of 1.02(36), -0.03(59), and -1.95(59) for (111), (100), and (110) crystal orientations, respectively. Even with the large error bars, these values clearly indicate that hard-sphere crystals with (100) orientations will only partially wet a hard wall at coexistence, and those with (110) orientations are nonwetting. However, the precision of these calculations is not sufficient to determine if the (111) crystal will exhibit complete wetting (cos θ > 1) or partial wetting (cos θ < 1). Recently, Fortini et al.20,21 employed an exponential barrier function for the thermodynamic integration to give γwf ) 1.990(14)kT/σ2 and γwc ) 1.457(18)kT/σ2 and 2.106(21)kT/σ2, for (111) and (100) orientations, under bulk coexistence conditions. Using Young’s equation, we find that these values give cos θ ) 0.96(4) and -0.21(4), for (111) and (100), respectively. The precision of these values represent an order of magnitude improvement over the earlier calculations of Heni and Lo¨wen, but give a value of cos θ111 that, within the error bars, cannot conclusively distinguish between partial and complete wetting.
Figure 1. Illustration of cleaving-wall method with a dividing surface at z ) 0. Top: Initial insertion of walls at (zw. Because the wall at -zw only interacts with particles with z > 0 and that at +zw only interacts with particles with z < 0, the insertion of the walls at zw g σ/2 requires zero energy. Bottom: Configuration of system after the walls have been moved to z ) 0.
2. Method In order to compute the excess free energy of a hard-sphere system at a planar hard wall, we adopt the cleaving-wall approach introduced in ref 9 and reviewed in ref 22. First, separate equilibrium bulk crystal and fluid systems are constructed at coexistence conditions10,20 (Fc ) 1.0377σ-3, Ff ) 0.9393σ-3, and p ) 11.57kBTσ-3). The equilibration and production simulations are performed using the hard-sphere molecular-dynamics algorithm of Rapaport.23 Periodic boundary conditions are employed in all directions. To minimize finite size effects, we have used a relatively large system size: Lx ≈ Ly ≈ 20σ, and Lz ≈ 40σ, with about 16 000 spheres. To insert a hard wall into the bulk crystal or fluid system at a specified location, say z ) 0, two planar hard walls are placed at ( zw. The interaction of the walls with the system is defined as follows: a system sphere with coordinates (xi,yi,zi) collides with the wall at zw (-zw) when zw - zi ) σ/2 and zi e 0 (-zw - zi ) -σ/2 and zi g 0), where σ is the sphere diameter. Initially, zw > σ/2, so the walls do not interact with the system, and no work is required for their insertion. The walls are then gradually moved toward each other until they reach the position zw ) 0. At this point, the hard spheres with zi < 0 can no longer collide with the spheres at zi > 0, and the bulk fluid or crystal is transformed into a system containing an impenetrable planar wall. Figure 1 shows the initial and final system configurations. Because of the periodic boundary conditions, the final state is equivalent to a slit-pore geometry with wall separation Lz. The average pressure on the walls as a function of wall position, pw(zw), is measured during the cleaving process. The work per unit area required to insert the wall, that is, the excess free energy, is obtained by numerically evaluating the integral
15954 J. Phys. Chem. C, Vol. 111, No. 43, 2007
Laird and Davidchack TABLE 1: Results for the Hard-Sphere Fluid/Hard-Wall Interfacial Free Energy, γwf as a Function of Pressure and Density up to the Coexistencea F* ) Fσ3
p* ) p σ3/kT
γwfσ2/kT
0.9393 (c) 0.935 0.930 0.925 0.920 0.9015 0.8327 0.7639 0.573 0.382 0.191
11.571(3) (c) 11.347(4) 11.092(3) 10.843(4) 10.600(5) 9.751(6) 7.181(6) 5.307(4) 2.283(2) 0.9207(6) 0.2908(3)
1.975(2) 1.951(3) 1.925(2) 1.898(3) 1.870(3) 1.773(5) 1.446(4) 1.174(2) 0.6437(8) 0.3225(4) 0.1233(2)
a The numbers in parenthesis are the errors (2σ) in the last digit shown. The label “(c)” indicates the coexistence value.
Figure 2. Determining the interfacial free energy for the wall-fluid, γwf, and wall-(111) crystal, γwc, systems at solid-fluid coexistence densities. The top plots show the average pressure on the walls as a function of wall position. The bottom plots show excess free energy computed according to eq 3 (solid dots) and eq 4 (open circles) for different system sizes. Dashed lines represent a straight line least-squares fit. In all plots, the estimated 95% confidence intervals are smaller than the size of the symbols.
γ)
∫0σ/2 pw(zw)dzw
(3)
Because the simulation is carried out at constant volume, the insertion of the walls will lead to increasing pressure in the bulk system and thus to an overestimation of the measured excess free energy. In order to take this into account, we perform a simulation of three systems with varying sizes in the z direction: Lz, (3/2)Lz, and 3Lz, and then determine the excess free energy by extrapolating the obtained results to Lz-1 ) 0. Another way to compensate for the increased pressure in the bulk is to measure the z component of the pressure in the bulk, pzz(zw), and then re-scale the pressure on the walls as follows
γ)
∫0σ/2 pw(zw)[pzz(σ/2)/pzz(zw)] dzw
(4)
The pressure pzz(zw) is measured in the middle half of the system away from the cleaving walls and increases by 1-2% from its value pzz(σ/2) for the unperturbed bulk system as the walls are moved to position zw As illustrated in Figure 2, the results obtained with either method fall on a straight line and thus can be fitted using least-squares linear regression and extrapolated to Lz-1 ) 0. The calculated intercepts for the two methods are consistent within the error bars and thus provide precise estimates of the wall-fluid and wall-(111) crystal excess free energies. The application of the cleaving-wall method to determine γwf and γwc for the hard-sphere/hard-wall system is far simpler and more precise than the corresponding calculation of the crystalfluid interfacial free energy, γcf. Whereas the calculation of γcf requires four steps (cleaving the fluid, cleaving the liquid, rearranging the boundary conditions, and removing the cleaving walls from the crystal-fluid interface), the calculation of either γwc or γwf requires only a single step to insert the walls. 3. Results and Discussion Using the method outlined above, we have determined the wall-fluid interfacial free energy for a number of densities up
Figure 3. Wall-fluid interfacial free energy, γwf, for the hard-sphere system as a function of density. The black circles are the results of the current study, and the open diamonds, open triangles and stars are the results from refs 19, 20 and 21, and 17, respectively. The solid line shows the prediction of scaled-particle theory (SPT).
to the coexistence density, Ffc ) 0.9393σ-3. The data are summarized in Table 1 and plotted in Figure 3 together with previous estimates by Heni and Lo¨wen,19 Fortini and Dijkstra20,21 and Miguel and Jackson.17 Also shown in Figure 3 is the prediction from the scaled particle theory (SPT):24,25,19 / γwf(SPT) ) γSPT + p/CS/2
(5)
where / γSPT )-
9η2(1 + η) 2π(1 - η)3
(6)
and P/CS is the Carnahan-Starling (CS) equation of state for the hard-sphere fluid, given by
p/CS )
6η(1 + η + η2 - η3) π(1 - η)3
(7)
and η ) πFσ3/6 is the packing fraction. For better comparison with the earlier data, Figure 4 shows the data for γwf in the high-density range 0.8 < Fσ3 < 1.0. Our current values for γwf are considerably more precise than earlier estimates for this quantity. The error bars in γwf are also smaller than those obtained in our calculations of γcf, which were also calculated using the cleaving-wall approach. In the
Wall-Induced Prefreezing in Hard Spheres
J. Phys. Chem. C, Vol. 111, No. 43, 2007 15955
Figure 4. Same as Figure 3, except restricted over the reduced density range (0.8 < Fσ3 < 0.95).
TABLE 2: Results for the Hard-Sphere Crystal/Hard-Wall Interfacial Free Energy, γwc, as a Function of Pressure and Density in the Vicinity of Coexistencea F* )
Fσ3
1.0377 (c) 1.0321 1.0249 1.0178 1.0111
p* ) p
σ3/kT
11.573(3) (c) 11.355(3) 11.087(3) 10.834(3) 10.606(3)
2 γ111 wc σ /kT
2 γ100 wc σ /kT
2 γ110 wc σ /kT
1.412(2) 1.394(2) 1.370(2) 1.349(2) 1.330(2)
1.915(2)
2.86(2)
a The numbers in parenthesis are the errors (2σ) in the last digit shown. The label “(c)” indicates the coexistence value.
latter case, the hysteresis is driven by the formation of a wetting layer of crystal at the wall-fluid interface induced by the structured cleaving walls necessary for that system. In our current wall-fluid calculations for structureless walls, nucleation of such a wetting layer is not observed over the duration of the simulations; this is consistent with the observed existence of a nucleation barrier to wall-induced crystallization of hardspheres.26 The absence of wall-induced crystal nucleation allows us to calculate γwf directly at coexistence, in contrast to earlier thermodynamic integration methods19,20 that required extrapolation from sub-coexistence density values to avoid crystallization. At low densities (Fσ3 < 0.8), all methods, including SPT, give good agreement for γwf within the error bars. At high density, our results are in very good agreement with the other high-precision results of Fortini and Dijkstra,27,20,21 as well as with those of Heni and Lo¨wen, which have rather large error bars. Both our results and those of Fortini and Dijkstra are significantly larger than the SPT prediction and the lowerprecision calculations of Miguel and Jackson.17 In their paper, Miguel and Jackson17 speculate that the positive deviation of the Fortini and Dijkstra estimates from the SPT estimate (and their estimates) are due to differences in the definition of system volume. However, such differences should disappear in the limit of large system sizes and are overcome by the use of finitesize scaling to extrapolate to infinite wall separation, as we have demonstrated here. Our current values indicate that the SPT underestimates the value of the fluid-wall interfacial free energy for densities near coexistence. Table 2 shows the values of the wall-crystal interfacial free energy, γwc, at the coexistence point for the (111), (100), and (110) orientations of the fcc hard-sphere crystal at the wall. For the (111) orientation, we have also determined values at subcoexistence pressures, which will be useful in understanding
Figure 5. Estimate of the wetting layer thickness h from eq 9. The circles use the value for γcf from ref 10, and the diamonds use the value for this quantity from ref 11.
the thermodynamics of precrystallization of the hard-sphere fluid just below coexistence. The low precision of the value of γ110 wc is due to the appearance of (111) layers at the wall during the insertion of walls into the (110) system. Because the value of γ110 wc is so large, it is likely that the formation of intermediate (111) layers between the wall and the (110) system is the result of the lower free energy of such a system compared with the wall-(110) crystal system. That is, the sum of the (111)/(110) 110 grain boundary free energy and γ111 wc is smaller than γwc . Using Young’s equation eq 1, together with the earlier estimate of γcf from ref 10, we can determine value of the wetting parameter at coexistence, cos θ, to be 1.013(14), 0.111(5), and -1.57(3). These results are consistent with earlier simulations in showing that partial wetting is predicted for the (100) orientation, whereas the (110) interface is nonwetting. However, the current results for the (111) interface are of sufficient precision to show that complete wetting is predicted for this interface, although only marginally so. We can also calculate the wetting parameter using another recent estimate of the (111) crystal-fluid interfacial free energy [γcf ) 0.546(16)] obtained using the capillary fluctuation method.11 This value is slightly lower than that calculated using the cleavingwall method; however, the precision is lower, and the two estimates agree within the error bars. Using this value for γcf gives a wetting parameter for the (111) crystal orientation of cos θ ) 1.03 ( 0.03, a value also consistent with marginal complete wetting. In principle, the thickness of the prefreezing layer can be determined from knowledge of the bulk and interfacial thermodynamic properties. The difference in Gibbs free energy per unit area, ∆g, between a hard-sphere-fluid/hard-wall system with a prewetting layer of crystal of thickness h and the same system with no wetting layer can be approximated as
∆g ) Fs∆µh - ∆γ + γ0 e-h/h0
(8)
where ∆µ ) µs - µl is the chemical potential difference between the bulk solid and fluid phases, ∆γ ) γwf - [γwc + γcf], and the final term represents the asymptotic form of the interaction between the wall-crystal and crystal-fluid interfaces.28 The interface-interface interaction term, which is responsible for the logarithmic dependence of the wetting layer thickness, is, unfortunately, unknown. If this term is ignored, we can get an upper-bound estimate of the wetting layer thickness by setting
15956 J. Phys. Chem. C, Vol. 111, No. 43, 2007
Laird and Davidchack
∆g to zero, which gives
h)
∆γ Fc∆µ
(9)
If the pressure is close to the coexistence pressure, pc, we can approximate ∆µ by a first-order expansion in (pc - p):
( )
Fc - F f (p - p) ∆µ ) FfFc c c
(10)
Use of eqs 9 and 10 require the pressure dependence of the interfacial free energies. We have calculated these for the wallfluid and wall-crystal systems, but the pressure dependence of γcf is not known. However, we can estimate γcf(p) by assuming that Turnbull’s scaling relation29 holds, giving
γcf ∝ ∆Hfus F-2/3 s
dynamics, a (111) crystal layer will completely wet the wallfluid interface near coexistence. The current work also demonstrates that the cleaving-wall approach is a highly efficient and precise method for the calculation of the interfacial free energy for a static surface in contact with a fluid, crystal, or glass. Acknowledgment. One of the authors (B.B.L.) gratefully acknowledges support from the National Science Foundation under Grant CHE0316127. The computations were carried out on the University of Leicester Mathematical Modeling Centre cluster (purchased through the EPSRC strategic equipment initiative). We also wish to thank Andrea Fortini for graciously providing us with data from his Ph.D. thesis. Finally, this work has benefited greatly from the participation of B.B.L. in meetings sponsored by the Computational Materials Science Network, a U.S. Department of Energy collaboration.
(11) References and Notes
This approximation has recently been used with much success in estimating nucleation rates for undercooled liquids.30 For the hard-sphere system, the latent heat of fusion is equal to ∆(pV) (that is, ∆E ) 0) and is given by
( )
∆Hfus Ff - F c ) p N FfFc
(12)
Using the data in Tables 1 and 2, the previously calculated values of γcf, and the approximations in eqs 10-12, it is possible to estimate the wetting layer thickness using eq 9. The resulting values are plotted in Figure 5 using both the cleaving-wall method10 and fluctuation method11 estimates for γcf. The error bars in Figure 5 are derived only from the estimated errors in the input quantities and do not include the likely systematic errors from the approximations in eqs 10-12; therefore, it is likely that the actual uncertainties in h are even larger. While we can make confident statements as to the existence of complete crystal wetting in this system, the precision of the interfacial free energies, although quite high, combined with our lack of information of the interface-interface interaction and the density dependence of the crystal-melt interfacial free energy, do not allow for an accurate estimation of the thickness of the wetting layer as a function of density. In summary, we have shown that our thermodynamic integration method based on cleaving walls can yield highly precise values for the wall-fluid and wall-crystal interfacial free energies for a system of hard spheres in contact with a hard wall. The calculated free energies are precise enough that, when coupled with our earlier estimates of the crystal-melt interfacial free energy,10,11 we can conclude that, on the basis of thermo-
(1) Alder, B.; Wainwright, T. J. Chem. Phys. 1957, 27, 1208. (2) Hansen, J.; McDonald, I. Theory of Simple Liquids, 2nd ed.; Academic Press: New York, 1986. (3) Courtemanche, D.; van Swol, F. Phys. ReV. Lett. 1992, 69, 2078. (4) Courtemanche, D.; van Swol, F. Mol. Phys. 1993, 80, 861-875. (5) Kegel, W. J. Chem. Phys. 2001, 115, 6538. (6) Dijkstra, M. Phys. ReV. Lett. 2004, 93, 108303. (7) Dijkstra, M.; van Roij, R. J. Phys.: Condens. Matter 2005, 17, 83507. (8) Dietrich, S. In Phase Transitions and Critical Phenomena; Domb, C., Liebowitz, J., Eds.; Academic Press: New York, 1988; pp 1-218. (9) Davidchack, R.; Laird, B. Phys. ReV. Lett. 2000, 85, 4751. (10) Davidchack, R.; Laird, B. Phys. ReV. Lett. 2005, 94, 086102. (11) Davidchack, R.; Morris, J.; Laird, B. J. Chem. Phys. 2006, 125, 094710. (12) Feng, X.; Laird, B. J. Chem. Phys. 2006, 124, 044707. (13) Hoyt, J.; Asta, M.; Karma, A. Phys. ReV. Lett. 2001, 86, 5530. (14) Davidchack, R.; Laird, B. J. Phys. Chem. 2003, 118, 7657. (15) Morris, J.; Song, X. J. Chem. Phys. 2003, 119, 3920. (16) Kirkwood, J.; Buff, F. J. Chem. Phys. 1949, 17, 338. (17) Miguel, E. D.; Jackson, G. Mol. Phys. 2006, 104, 3717-3734. (18) Tiller, W. The Science of Crystallization: Microscopic Interfacial Phenomena; Cambridge University Press: New York, 1991. (19) Heni, M.; Lo¨wen, H. Phys. ReV. E 1999, 60, 7057. (20) Fortini, A.; Dijkstra, M. J. Phys.: Condens. Matter 2006, 18, L371. (21) Fortini, A. Ph.D. Thesis, Utrecht University, The Netherlands, 2007. (22) Laird, B.; Davidchack, R. J. Phys. Chem. B 2005, 109, 17802. (23) Rapaport, D. C. The Art of Molecular Dynamics Simulation, 2nd ed.; Cambridge University Press: New York, 2004. (24) Barker, J.; Henderson, D. ReV. Mod. Phys. 1976, 48, 587. (25) Reiss, H.; Frisch, H.; Helfand, E.; Lebowitz, J. J. Chem. Phys. 1960, 32, 119. (26) Auer, S.; Frenkel, D. J. Chem. Phys. 2003, 119, 7467. (27) Fortini, A.; Dijkstra, M.; Schmidt, M.; Wessels, P. Phys. ReV. E 2005, 71, 051403, 13 pages. (28) Heni, M.; Lo¨wen, H. J. Phys.: Condens. Matter 1002, 13, 4675. (29) Turnbull, D. J. Appl. Phys. 1950, 21, 1022. (30) Aga, R.; Morris, J.; Hoyt, J.; Mendelev, M. Phys. ReV. Lett. 2006, 96, 245701.
