PHYSICAL REVIEW E 89, 032412 (2014)
Wall-induced phase transition controlled by layering freezing Huijun Zhang,1,2 Shuming Peng,2 Xinggui Long,2 Xiaosong Zhou,2 Jianhua Liang,2 Chubin Wan,1 Jian Zheng,2,3 and Xin Ju1,* 1 School of Mathematics and Physics, University of Science and Technology Beijing, Beijing 100083, China Institute of Nuclear Physics and Chemistry, China Academy of Engineering Physics, Mianyang 621900, China 3 School of Nuclear Science and Technology, University of Science and Technology of China, Hefei 230026, China (Received 30 October 2013; revised manuscript received 27 December 2013; published 31 March 2014) 2
Molecular dynamics simulations of the Lennard-Jones model are used to study phase transitions at a smooth surface. Our motivation is the observation that the existence of an attractive wall facilitates crystallization. To investigate how this wall influences phase transitions, the strength of wall-particle interaction is varied in our studies. We find that the phase behavior depends on the strength parameter α, i.e., the ratio between wall-particle and the particle-particle attraction strength. Three critical values of the ratio, namely, αp , αw , and αc , are used to define the qualitative nature of the phase behaviors at a smooth surface. Some interesting phenomena due to the increase of α are observed. First, a set of close-packed planes, i.e., {111} planes in fcc structures or {0001} planes in hcp structures, are “rotated” from intersecting to parallel to the wall when α = αp ; second, the layering phase transition close to the wall antecedes that of the bulk when α = αw . Finally, the first-order phase transition in the first two layers is supplanted by a continuous phase transition when α = αc , which to some extent can be treated as a quasi-two-dimensional process. We find that bulk freezing always discontinuously occurs through a first-order phase transition, and seems to be isolated from the freezing process occurring close to the attractive surfaces. Moreover, during the heating process, we observe minimal dependence at a strongly attractive surface. DOI: 10.1103/PhysRevE.89.032412
PACS number(s): 81.10.−h, 64.70.dg, 68.35.Rh, 64.70.dj
I. INTRODUCTION
Crystallization has a pivotal function in many natural and industrial processes, such as the production of pharmaceuticals and functional materials, formation of ice and frost, and formation of minerals and bone tissues. These highordered and three-dimensional crystals form spontaneously from the disordered matrix phase by the earliest rare event, i.e., nucleation [1–3]. This phase transition almost always begins with the formation of the new phase at a surface [4–6]. Understanding the effects of the smooth surfaces on the phase transition is still limited because of experimental confinements, such as harsh conditions [7] and accurate observation techniques [8]. However, these limitations do not bound computer simulations. In this article, we report on a study of the Lennard-Jones (LJ) systems that are adequate to reproduce the thermodynamic and kinetic behaviors of classical fluids. Phase transition at a surface presents various scenarios depending on the strength of wall-particle interactions [9–12]. For strongly attractive surfaces, wettability induces highly ordered layers near the wall [13,14]. Such highly quasilong-range translational order and long-range orientational order of layers can, to some extent, be treated as a twodimensional (2D) system. The generally accepted theory about 2D solid-liquid phase transition is the Kosterlitz-ThoulessHalperin-Nelson-Young (KTHNY) theory [15], which states that melting in 2D systems may involve two continuous transitions. The first is from solid to a hexatic phase [16] through the unbinding of dislocation pairs, and the second is from the hexatic phase to liquid through the unbinding of disclination pairs [17]. Although the mechanism described in
*
Corresponding author: [email protected]
1539-3755/2014/89(3)/032412(7)
the KTHNY theory seems rather general and appealing, it is still controversial whether this is the only process occurring in the 2D melting [15,18]. Page and Sear [11] recently used Monte Carlo simulations of LJ models to research phase transition at a surface. They found that freezing in bulk is controlled by prefreezing (freezing that occurs above the melting point). Moreover, prefreezing could abolish the nucleation barrier to bulk freezing at a strongly attractive surface. Nevertheless, they have not determined whether prefreezing proceeds through the KTHNY mechanism. Radhakrishnan et al. [9,16] reported that a continuous phase transition occurs through a hexatic phase at strongly attractive surfaces of pores, and this phase transition could involve several layers. Grivova et al. [12] studied slit pore systems that accommodate single or double layers. They found that a continuous phase transition exists only in a monolayer of confined particles. They also stated that the controversy concerning the observation of a hexatic phase derives from the different order parameters chosen and from the observation techniques employed. Nevertheless, the slit pore systems they used are affected by confinement effects [19]. These studies cannot clearly explain phase behaviors at a surface without confinement. In the present work we observe, via numerical simulations, that several distinct scenarios regarding the phase transition of the layers near the wall depend on the strength of the wall-particle attraction. These scenarios are defined as “parallel to the wall,” “complete wetting,” and “continuous transition.” We also find that continuous transition can emerge only in the first two layers. All these findings may better elucidate phase transitions at a smooth surface. II. SIMULATION DETAILS
We investigate the crystallization process of a fluid in the presence of a smooth surface by means of molecular dynamics (MD) simulations. The pair interaction between fluid particles 032412-1
©2014 American Physical Society
ZHANG, PENG, LONG, ZHOU, LIANG, WAN, ZHENG, AND JU
PHYSICAL REVIEW E 89, 032412 (2014)
is modeled with the standard Lennard-Jones (LJ) potential, σ 12 σ 6 , (1) − U (r) = 4ε r r where r is the distance between particles, σ is the diameter of a particle, and ε is the strength of interaction between particles. The wall-particle interaction is obtained from the integral of the classical 12-6 LJ interaction over an infinite surface [20,21], 3 2 σ 9 σ UW (z) = εW − , (2) 15 z z where εW is the strength of the interaction, and z is the distance between a particle and the wall. Well depth for wall-particle interaction is approximated as εw . Both LJ potentials are truncated at 4.7σ . Herein, we use ε as the unit of energy, σ as the unit of length, and ε/kB as the unit of temperature (kB is the Boltzmann constant). We refer to a homogeneous system [22] wherein 3400 particles in a cubic box (Lx = Ly = Lz = 17.62σ ) provide a reduced density of 0.62. The attractive wall is set at z = 0 while a reflective wall is placed at z = Lz . Periodic boundary conditions are applied along the x and y axis. In the simulations, the Nose-Hoover thermostat is used to control temperature in the canonical ensemble (NVT). A strength parameter (SP) α (α = εW /ε) that is the ratio of the wall-particle to the particle-particle well depth of interaction potential, is defined to describe the walls. We vary α (keeping ε fixed) from 1 to 12 to study the wall-induced phase transition. For each α, the reduced temperature T is decreased from 0.83 to 0.43 by an interval of 0.04, and 2.0 × 105 and 1.0 × 105 time steps are, respectively, used to reach equilibrium and sample in each interval. To analyze the structures of crystalline phases in simulations, we employ the common neighbor analysis (CNA) technique [23]. This method identifies particles by considering the number and connectivity of the neighbors shared by two neighboring particles. Finally, we implement our simulations using the LAMMPS software package [24]. III. RESULTS AND DISCUSSION
We first plot the curves of reduced potential energy U versus reduced temperature T , as shown in Fig. 1. For 1< = α< = 4, we observe sharp drops U in the potential energy (U 2800) as T decreases, which manifests a first-order phase transition [19,21]. For strongly attractive surfaces, the onset T at which the energy drop U takes place shifts towards higher values, implying the facilitation of freezing. In addition, the transition range (i.e., the range of reduced T in which U is observed) extends from one to several reduced temperature intervals when α exceeds 5. This behavior indicates that the phase transition mechanism at a strongly attractive wall differs from that at a weakly attractive wall. Four possible hypotheses can explain the widening of the transition range: (i) During the freezing process particles are trapped into a kinetically metastable phase without crystallization; hence, no liquid-solid phase transition occurs, and the sudden drop in U disappears. (ii) Particles freeze through a continuous phase transition, similar to the prefreezing phenomenon at the strongly attractive wall [11]
FIG. 1. (Color online) Reduced potential energy U as a function of reduced temperature T for increased SP from α = 1 to 12. The black dashed line is the result of a reference system without the foreign wall. With cooling from T = 0.83 to 0.43, a sudden drop (U 2800) can be observed in U for 1< = α< = 4, manifesting a first-order phase transition. These drops shrink as α continuously increases, which is confirmed by the extension of the temperature range.
and the wetting case predicted by classical nucleation theory [25]. (iii) A cascade of layering first-order phase transitions completes the liquid-solid transition; consequently, the sudden drops are flattened by several drop parts. (iv) The liquid-solid transition involves both continuous and discontinuous phase transition steps that cause the gradual drop. To identify the correct scenario among our hypotheses, we have analyzed the structural properties of the LJ system after the temperature quench, for 1< = α< = 12. By tracking the trajectories from T = 0.83 to 0.43, we observe that the system is always converted into a solid composed of close-packed planes (CPPs). This observation of CPPs is consistent with the simulations of Malley et al. [26]. Moreover, the wall prefers parallel CPPs when SP changes from α = 1 to 2 (Fig. 2); this scenario is specified as “parallel to the wall.” The LJ walls can create a potential well that could trap particles. Thus, a deeper well depth corresponds to more trapped particles. The parallel CPPs [Fig. 2(b)] have more trapped particles than the intersecting state [Fig. 2(a)], which may explain why the
FIG. 2. (Color online) Simulation snapshots at reduced temperature T = 0.43. Yellow (light gray) and green (gray) balls represent fcc and hcp particles, respectively, while liquidlike particles are not shown. The black walls represent the LJ surfaces. The final crystals consist of fcc and hcp form layers, and the layers cross the wall for α = 1 (a) and parallel to the wall for α = 2 (b). The individual images are created using Atomeye [27].
032412-2
WALL-INDUCED PHASE TRANSITION CONTROLLED BY . . .
PHYSICAL REVIEW E 89, 032412 (2014)
stronger attractive wall prefers parallel CPPs. Consequently, a critical value αp between 1 and 2 exists, above which the CPPs are always parallel to the wall [Fig. 2(b)]. Moreover, for α = 1, CPPs with dominant hcp arrays cross the wall at an approximately 62° angle [Fig. 2(a)], and this fixed degree serves as the optimum choice in crystal nucleation and growth of pure hcp structure [28]. We infer that this fixed angle may preferably change to 70.5° when the dominance reverses. This inference is due to the fact that the 70.5° angle is the intrinsic angle consisting of two sets of CPPs in perfect fcc structures [28,29]. Considering that the wall-particle attraction quickly decays to zero as the distance from the wall increases, the phase behaviors of the layers near the wall may essentially differ from those in the bulk. Thus, we mainly examine phase behaviors in
the first three layers. The layering radial distribution function (RDF) [17] gxy (r ) and the bulk RDF g(r) are then investigated. The bulk RDF g(r) is defined as g(r) =
dn(r) , ρ
(3)
where dn(r) presents the average local number at distance r from a fixed particle and ρ is the overall number density. The layering RDF gxy (rxy ) is defined as
i j
i,i=j δ rxy δ rxy − rxy gxy (rxy ) = , (4) ρxy where ρxy is the layering number density and rxy is the distance between the particles in the xy plane. Comparison of the
FIG. 3. (Color online) Layering RDFs. The sixth layer RDFs for α = 2, 4, 5, and 6 are shown in (a), (b), (d), and (f) to represent the bulk situation. (c) and (e) are layering RDFs in the second layer for α = 5 and 6, respectively, and the insets are layering RDFs in the first layer correspondingly. These RDFs all exhibit characteristic features of liquids or solids depending on temperature. In each of them, an abrupt jump occurs from liquid to solid as the temperature reaches the critical value marked in red (gray), indicating the occurrences of a first-order phase transition. When α increases to 6, the transition range of bulk and that of the first two layers are not synchronized. 032412-3
ZHANG, PENG, LONG, ZHOU, LIANG, WAN, ZHENG, AND JU
PHYSICAL REVIEW E 89, 032412 (2014)
FIG. 4. (Color online) RDFs for α = 7 and 8. The insets are parts of the first and third layers. For α = 7, the first two layers are discontinuously converted into solid through a first-order transition (a). For α = 8, a continuous phase transition exclusively occurs in the first two layers (c). The freezing of the sixth and third layers always proceeds through the first-order phase transition (b) and (d). The reduced transition temperatures are marked in red (gray).
RDFs of the sixth layer and those of the bulk reveals that they exhibit the same features (see Supplemental Material [30]). Then, the RDFs of the sixth layer are used to substitute the bulk. In Fig. 3, we observe abrupt jumps in the peaks as the reduced temperature decreases. Before the jump, RDF peaks consist of a prominent peak and a rapidly decaying peak before decaying to unit 1. This behavior exhibits the isotropic nature of liquids. Conversely, when the temperature reaches the critical value (red marks in Fig. 3), the RDFs show a long-range order that is the feature of the crystal. Thus, these jumps are of a first-order phase transition representing liquid-solid conversion. Moreover, the order of the system that has just finished the jump weakens as α increases. Such an order-limited system needs further evolution to reach the final stable state because of the sluggishness [31,32]. When α reaches 6, the freezing of the first two layers and that of the bulk do not synchronize with temperature [Figs. 3(e) and 3(f)], and the so-called “complete wetting” occurs [13,14]. Hence, a critical value αw must exist between 5 and 6. This value distinguishes the phase transition of the first two layers from that of the bulk. Furthermore, the reduced temperatures for the bulk freezing steadily increase and then stabilize at T = 0.67. Such rises in temperature correspond to the shift of transition ranges in the energy profiles (Fig. 1). Treatment of the first two layers as a quasi-2D system, with their freezing anteceding the bulk, is yet to be verified.
Since the particular phase behaviors in the first two layers depend on the SP α, the layering phase behaviors at the more attractive walls are then investigated. Figure 4 shows still remarkable jumps in RDF peaks at α = 7 [Figs. 4(a) and 4(b)]. Nevertheless, these jumps occur within different temperature ranges. In the first two layers, peaks jump from T = 0.79 to 0.75, corresponding to a small drop (U 995) in U as shown in Fig. 1. In the third and in the bulk layers, peaks jump from T = 0.75 to 0.71 and from T = 0.71 to 0.67, respectively. Notably when α reaches 8, these jumps disappear in the first and second layers [Fig. 4(c)], indicating that freezing in the first two layers occurs through a continuous phase transition. We specifically name this scenario as “continuous transition.” However, the third and the bulk RDFs still exhibit the features of the first-order phase transition [Fig. 4(d)]. Therefore, a critical value αc must exist between 7 and 8, beyond which the first two layers continuously freeze. This continuous phase transition differs from the mechanism of phase behavior in three dimensions and can be treated as a quasi-2D phase transition [15]. Page and Sear [11] stated that the energy barrier for prefreezing disappears at the strongly attractive surface, consistent with our observation (continuous transition). The difference is that in our study the continuous transition is specifically limited in the first two layers. Although the wall strength they used is absolutely different from ours, the SP α is really close to ours. This finding
032412-4
WALL-INDUCED PHASE TRANSITION CONTROLLED BY . . .
PHYSICAL REVIEW E 89, 032412 (2014)
FIG. 5. (Color online) NDPs of the first three layers for α = 7 (a) and α = 8 (b). They all exhibit excellent layering at the lower temperature. At each reduced temperature, decreased distance from the wall corresponds to sharper layers. For α = 7, the first three layers all exhibit an abrupt change, and their critical temperatures are T = 0.75, 0.75, and 0.71, respectively. By contrast, for α = 8, continuous changes replace the abrupt states in the first two layers, whereas the third layer retains a discontinuous transition at 0.71.
may imply that the SP α we introduced can be used as a general parameter to describe phase transition at a smooth surface in LJ systems. To further confirm the existence of continuous transition, we examine the number density profile (NDP). In Fig. 5, we can see that peaks become sharper as the reduced temperature decreases. The sharper peaks represent the higher-ordered layers parallel to the wall, and the particles in such layers may approach the 2D situations. Generally, the extending of peaks can be ascribed to thermal fluctuations. Therefore, the wide peaks in the third layer imply that it is relatively far away from the 2D situation. For α = 7 [Fig. 5(a)], the peaks narrow and rise abruptly in the first two layers as the reduced temperature decreases from 0.79 to 0.75. Moreover, these sudden changes are delayed to the next reduced temperature interval (0.75–0.71) in the third layer. By contrast, when α = 8, gradual changes of peaks exclusively occur in the first two layers [Fig. 5(b)], confirming the continuous transition. These observations are consistent with Fig. 4. The mobility of the particles in the first two layers is then considered. The most straightforward method is to estimate the number of particles in each layer. These results are plotted in Fig. 6. The curves undergo a continuous increase at the early stage and then abruptly jump to a constant value for α = 2 and 7 [Figs. 6(b) and 6(c)]. These abrupt increases of the particles in the layers correspond to the occurrences of the first-order phase transitions for which replenishments of particles from the bulk are necessary. Hypothetically, when the attractive force from the wall is sufficiently strong (limiting condition), more particles are confined within a narrow region [inset of Fig. 6(d)] and ought to form a relatively high-ordered array. In this case, as the temperature decreases, the highly ordered layers may solidify without changing their relatively ordered configuration to some extent. In addition, this transition may not need more replenishment from neighboring particles, and a continuous phase transition thus occurs. Actually, curves in α = 8 show a continuously small rise, and are consistent with
this limiting assumption [Fig. 6(d)]. However, the relatively huge and continuous increases for α = 1 [Fig. 6(a)] do not represent a continuous phase transition. Evidently, such increases should be attributed to the freezing with nonparallel layers near the wall [Fig. 2(a) and inset of Fig. 6(a)]. Finally, reverse procedures (melting) are performed using the final freezing results. We select the midpoint (Ttm ) of transition ranges for both freezing and melting processes because of the phase transitions across a relatively wide temperature range. As shown in Fig. 7, three regions (each later region is included in the former) marked in green, gray, and yellow correspond to the scenarios when SP (α) reaches the critical value αp , αw , and αc , respectively. We observe the considerable hysteresis between freezing and melting procedures. Such hysteresis diminishes as α increases, which represents the reduction of the energy barrier to the bulk freezing and the weakening of the first-order transition. Nevertheless, this diminution terminates when α reaches αc , which signifies that the bulk nucleation barrier cannot be thoroughly abolished [13]. In addition, phase behaviors above αc (yellow region) appear to be similar because Ttm hardly varies with α. In this region, α is sufficiently large to make the first two layers behave as a quasi-2D system. Such quasi-2D system possesses a denser and more ordered array of particles [Fig. 5(d) and inset of Fig. 6(d)], which ensures their lower mobility in the out-of-plane direction. In addition, the wall-particle interaction decays to a minute value in the bulk. Therefore, bulk freezing is reasonably expected to proceed in an isolated manner even through it starts from the wall. Several runs are performed for both the freezing and the melting processes. Data for freezing are coincident, whereas those for melting exhibit large fluctuations. These phenomena can be due to the fact that freezing always starts from the surface [insets of Figs. 6(a) and 6(d)] but melting begins from the exterior where fcc and hcp CPPs array with distinct sequence in each run. We therefore conclude that freezing depends on the walls, whereas melting depends on the
032412-5
ZHANG, PENG, LONG, ZHOU, LIANG, WAN, ZHENG, AND JU
PHYSICAL REVIEW E 89, 032412 (2014)
FIG. 6. (Color online) Number of particles in the first two layers. Arrows represent the transition tendency (vertical for abrupt increase and oblique for gradual increase). For abrupt changes, we mark the critical reduced temperature in red (gray). The insets are the snapshots of early stages for α = 1 and 8, and the details are the same as those in Fig. 2.
structures. This conclusion can be supported by Fig. 7 in which freezing is sensitive but the melting is insensitive to α. To test the finite size effects, a larger simulation box (Lx = Ly = Lz = 29.37σ ) is used. Results show few finite size effects (Fig. 7). Recalling the previous hypotheses, the first and the second are excluded, whereas the third and the fourth characterize the true
FIG. 7. (Color online) Midpoint of transition range as a function of SP α. Freezing and melting procedures are shown. The black dashed line denotes the triple point. Three regions are shown in green (light gray), gray (gray), and yellow (white) for the three scenarios described as parallel to the wall, complete wetting, and continuous transition, respectively. The open marks represent the results of the larger systems used to investigate the finite size effects.
scenarios of phase transitions at strongly attractive surfaces in our systems. IV. CONCLUSIONS
We have performed computer simulations to study the liquid-solid transition of a LJ system in the presence of an attractive surface. We identify different scenarios for phase transitions depending on the value of the parameter α, i.e., the ratio between the wall-particle and the particle-particle attraction strength. For weakly attractive walls, crystallization occurs at the walls with nonparallel CPPs. As the wall strength increases to the first critical value αp , “parallel to the wall” emerges. Throughout the simulation box crystallization occurs at the same temperature via a first-order transition when α < αw . When α = αw , the first two layers near the wall exhibit advanced freezing, and complete wetting occurs. Furthermore, for α>αc , the first two layers undergo freezing not through first-order phase transition, but through a continuous phase transition, which is also known as continuous transition This phenomenon has been observed to take place always in the first two layers, while in the bulk freezing still occurs via a first-order transition. Hence, the two layers can be treated as a quasi-2D system. Their liquid-solid transition may continuously occur through a hexatic phase, because Wierschem and Manousakis [33] stated that phase transition in slit pore with a soft core potential is expected to proceed through the KTHNY mechanism. However, we cannot determine whether the KTHNY mechanism can explain our results because distinguishing a hexatic phase from others requires a very large
032412-6
WALL-INDUCED PHASE TRANSITION CONTROLLED BY . . .
system size and a very long time, let alone the imperfection of the current techniques [33]. In the third layer, liquid-solid transition always proceeds through first-order phase transition and can occur earlier than the bulk freezing. This finding implies that the third layer may serve as the transient regime between the quasi-2D layers and the bulk phase. We also find that a considerable energy barrier to the bulk freezing always exists because the hysteresis cannot be completely abolished.
[1] P. R. ten Wolde and D. Frenkel, Science 277, 1975 (1997). [2] R. P. Sear, J. Chem. Phys. 129, 164510 (2008). [3] T. H. Zhang and X. Y. Liu, Angew. Chem., Int. Ed. 48, 1308 (2009). [4] A. Cacciuto, S. Auer, and D. Frenkel, Nature 428, 404 (2004). [5] M. Heni and H. L¨owen, Phys. Rev. Lett. 85, 3668 (2000). [6] W. S. Xu, Z. Y. Sun, and L. J. An, J. Chem. Phys. 132, 144506 (2010). [7] L. M. Ghiringhelli, C. Valeriani, E. J. Meijer, and D. Frenkel, Phys. Rev. Lett. 99, 055702 (2007). [8] P. R. ten Wolde, M. J. Ruiz-Montero, and D. Frenkel, J. Chem. Phys. 104, 9932 (1996). [9] R. Radhakrishnan, K. E. Gubbins, and M. SliwinskaBartkowiak, J. Chem. Phys. 116, 1147 (2002). [10] C. Alba-Simionesco, B. Coasne, G. Dosseh, G. Dudziak, K. E. Gubbins, R. Radhakrishnan, and M. Sliwinska-Bartkowiak, J. Phys.: Condens. Matter 18, R15 (2006). [11] A. J. Page and R. P. Sear, Phys. Rev. E 80, 031605 (2009). [12] N. Gribova, A. Arnold, T. Schilling, and C. Holm, J. Chem. Phys. 135, 054514 (2011). [13] S. Auer and D. Frenkel, Phys. Rev. Lett. 91, 015703 (2003). [14] M. Dijkstra, Phys. Rev. Lett. 93, 108303 (2004). [15] K. J. Strandburg, Rev. Mod. Phys. 60, 161 (1988). [16] R. Radhakrishnan, K. E. Gubbins, and M. SliwinskaBartkowiak, Phys. Rev. Lett. 89, 076101 (2002). [17] A. H. Marcus and S. A. Rice, Phys. Rev. E 55, 637 (1997). [18] J. G. Dash, Rev. Mod. Phys. 71, 1737 (1999). [19] D. Koga, H. Tanaka, and X. C. Zeng, Nature 408, 564 (2000). [20] C. Y. Lee, J. A. McCammon, and P. J. Rossky, J. Chem. Phys. 80, 4448 (1984). [21] P. Kumar, S. V. Buldyrev, F. W. Starr, N. Giovambattista, and H. E. Stanley, Phys. Rev. E 72, 051503 (2005). [22] In the homogeneous system, some characteristic nuclei could form from a homogeneous matrix. To ensure that a stand-alone
PHYSICAL REVIEW E 89, 032412 (2014) ACKNOWLEDGMENTS
This work was supported by Institute of Nuclear Physics and Chemistry (INPC), China Academy of Engineering Physics. We thank Nadezhda Gribova and Li Mao for the fruitful discussions about the layering radial distribution function. We also thank Jing Xu from Peking University for assistance with the manuscript.
[23] [24] [25] [26] [27] [28]
[29] [30]
[31] [32]
[33]
032412-7
nucleus matures without spatial restriction in most directions, we choose 3400 particles in a cubic box (with length L = 17.62). Such unimpeded nucleation and growth can provide us its formation mechanism at the atomic level. J. D. Honeycutt and H. C. Andersen, J. Phys. Chem. 91, 4950 (1987). S. Plimpton, J. Comput. Phys. 117, 1 (1995). R. P. Sear, J. Phys.: Condens. Matter 19, 033101 (2007). B. O’Malley and I. Snook, Phys. Rev. Lett. 90, 085702 (2003). Atomeye: Atomistic Configuration Viewer, http://li.mit.edu/ Archive/Graphics/A/. A. J. Page and R. P. Sear, J. Am. Chem. Soc. 131, 17550 (2009). Four different sets of close-packed planes (CPPs) exist in the fcc structure, but only one exists in the hcp state. We find that two sets of such CPPs in perfect fcc structures mutually intersect with an approximate angle of 70.53°, whereas the CPPs and the {1 − 10 − 1} planes in the hcp structures provide an angle of 62.06°. Such two angles are consistent with the work of Page and Sear [11]. Hence, we infer that these two intrinsic angles have a significant function in crystal nucleation and growth, which can help fcc (with 70.53° ridge) and hcp (with 62.06° ridge) blocks to form and grow without restriction. B. W. van de Waal, Phys. Rev. Lett. 67, 3263 (1991). See Supplemental Material at http://link.aps.org/supplemental/ 10.1103/PhysRevE.89.032412 for comparison of RDFs of the sixth layer and those of the bulk. W. A. Caldwell, J. H. Nguyen, B. G. Pfrommer, F. Mauri, S. G. Louie, and R. Jeanloz, Science 277, 930 (1997). H. Cynn, C. S. Yoo, B. Baer, V. Iota-Herbei, A. K. McMahan, M. Nicol, and S. Carlson, Phys. Rev. Lett. 86, 4552 (2001). K. Wierschem and E. Manousakis, Phys. Rev. B 83, 214108 (2011).
Wall-induced phase transition controlled by layering freezing Huijun Zhang,1,2 Shuming Peng,2 Xinggui Long,2 Xiaosong Zhou,2 Jianhua Liang,2 Chubin Wan,1 Jian Zheng,2,3 and Xin Ju1,* 1 School of Mathematics and Physics, University of Science and Technology Beijing, Beijing 100083, China Institute of Nuclear Physics and Chemistry, China Academy of Engineering Physics, Mianyang 621900, China 3 School of Nuclear Science and Technology, University of Science and Technology of China, Hefei 230026, China (Received 30 October 2013; revised manuscript received 27 December 2013; published 31 March 2014) 2
Molecular dynamics simulations of the Lennard-Jones model are used to study phase transitions at a smooth surface. Our motivation is the observation that the existence of an attractive wall facilitates crystallization. To investigate how this wall influences phase transitions, the strength of wall-particle interaction is varied in our studies. We find that the phase behavior depends on the strength parameter α, i.e., the ratio between wall-particle and the particle-particle attraction strength. Three critical values of the ratio, namely, αp , αw , and αc , are used to define the qualitative nature of the phase behaviors at a smooth surface. Some interesting phenomena due to the increase of α are observed. First, a set of close-packed planes, i.e., {111} planes in fcc structures or {0001} planes in hcp structures, are “rotated” from intersecting to parallel to the wall when α = αp ; second, the layering phase transition close to the wall antecedes that of the bulk when α = αw . Finally, the first-order phase transition in the first two layers is supplanted by a continuous phase transition when α = αc , which to some extent can be treated as a quasi-two-dimensional process. We find that bulk freezing always discontinuously occurs through a first-order phase transition, and seems to be isolated from the freezing process occurring close to the attractive surfaces. Moreover, during the heating process, we observe minimal dependence at a strongly attractive surface. DOI: 10.1103/PhysRevE.89.032412
PACS number(s): 81.10.−h, 64.70.dg, 68.35.Rh, 64.70.dj
I. INTRODUCTION
Crystallization has a pivotal function in many natural and industrial processes, such as the production of pharmaceuticals and functional materials, formation of ice and frost, and formation of minerals and bone tissues. These highordered and three-dimensional crystals form spontaneously from the disordered matrix phase by the earliest rare event, i.e., nucleation [1–3]. This phase transition almost always begins with the formation of the new phase at a surface [4–6]. Understanding the effects of the smooth surfaces on the phase transition is still limited because of experimental confinements, such as harsh conditions [7] and accurate observation techniques [8]. However, these limitations do not bound computer simulations. In this article, we report on a study of the Lennard-Jones (LJ) systems that are adequate to reproduce the thermodynamic and kinetic behaviors of classical fluids. Phase transition at a surface presents various scenarios depending on the strength of wall-particle interactions [9–12]. For strongly attractive surfaces, wettability induces highly ordered layers near the wall [13,14]. Such highly quasilong-range translational order and long-range orientational order of layers can, to some extent, be treated as a twodimensional (2D) system. The generally accepted theory about 2D solid-liquid phase transition is the Kosterlitz-ThoulessHalperin-Nelson-Young (KTHNY) theory [15], which states that melting in 2D systems may involve two continuous transitions. The first is from solid to a hexatic phase [16] through the unbinding of dislocation pairs, and the second is from the hexatic phase to liquid through the unbinding of disclination pairs [17]. Although the mechanism described in
*
Corresponding author: [email protected]
1539-3755/2014/89(3)/032412(7)
the KTHNY theory seems rather general and appealing, it is still controversial whether this is the only process occurring in the 2D melting [15,18]. Page and Sear [11] recently used Monte Carlo simulations of LJ models to research phase transition at a surface. They found that freezing in bulk is controlled by prefreezing (freezing that occurs above the melting point). Moreover, prefreezing could abolish the nucleation barrier to bulk freezing at a strongly attractive surface. Nevertheless, they have not determined whether prefreezing proceeds through the KTHNY mechanism. Radhakrishnan et al. [9,16] reported that a continuous phase transition occurs through a hexatic phase at strongly attractive surfaces of pores, and this phase transition could involve several layers. Grivova et al. [12] studied slit pore systems that accommodate single or double layers. They found that a continuous phase transition exists only in a monolayer of confined particles. They also stated that the controversy concerning the observation of a hexatic phase derives from the different order parameters chosen and from the observation techniques employed. Nevertheless, the slit pore systems they used are affected by confinement effects [19]. These studies cannot clearly explain phase behaviors at a surface without confinement. In the present work we observe, via numerical simulations, that several distinct scenarios regarding the phase transition of the layers near the wall depend on the strength of the wall-particle attraction. These scenarios are defined as “parallel to the wall,” “complete wetting,” and “continuous transition.” We also find that continuous transition can emerge only in the first two layers. All these findings may better elucidate phase transitions at a smooth surface. II. SIMULATION DETAILS
We investigate the crystallization process of a fluid in the presence of a smooth surface by means of molecular dynamics (MD) simulations. The pair interaction between fluid particles 032412-1
©2014 American Physical Society
ZHANG, PENG, LONG, ZHOU, LIANG, WAN, ZHENG, AND JU
PHYSICAL REVIEW E 89, 032412 (2014)
is modeled with the standard Lennard-Jones (LJ) potential, σ 12 σ 6 , (1) − U (r) = 4ε r r where r is the distance between particles, σ is the diameter of a particle, and ε is the strength of interaction between particles. The wall-particle interaction is obtained from the integral of the classical 12-6 LJ interaction over an infinite surface [20,21], 3 2 σ 9 σ UW (z) = εW − , (2) 15 z z where εW is the strength of the interaction, and z is the distance between a particle and the wall. Well depth for wall-particle interaction is approximated as εw . Both LJ potentials are truncated at 4.7σ . Herein, we use ε as the unit of energy, σ as the unit of length, and ε/kB as the unit of temperature (kB is the Boltzmann constant). We refer to a homogeneous system [22] wherein 3400 particles in a cubic box (Lx = Ly = Lz = 17.62σ ) provide a reduced density of 0.62. The attractive wall is set at z = 0 while a reflective wall is placed at z = Lz . Periodic boundary conditions are applied along the x and y axis. In the simulations, the Nose-Hoover thermostat is used to control temperature in the canonical ensemble (NVT). A strength parameter (SP) α (α = εW /ε) that is the ratio of the wall-particle to the particle-particle well depth of interaction potential, is defined to describe the walls. We vary α (keeping ε fixed) from 1 to 12 to study the wall-induced phase transition. For each α, the reduced temperature T is decreased from 0.83 to 0.43 by an interval of 0.04, and 2.0 × 105 and 1.0 × 105 time steps are, respectively, used to reach equilibrium and sample in each interval. To analyze the structures of crystalline phases in simulations, we employ the common neighbor analysis (CNA) technique [23]. This method identifies particles by considering the number and connectivity of the neighbors shared by two neighboring particles. Finally, we implement our simulations using the LAMMPS software package [24]. III. RESULTS AND DISCUSSION
We first plot the curves of reduced potential energy U versus reduced temperature T , as shown in Fig. 1. For 1< = α< = 4, we observe sharp drops U in the potential energy (U 2800) as T decreases, which manifests a first-order phase transition [19,21]. For strongly attractive surfaces, the onset T at which the energy drop U takes place shifts towards higher values, implying the facilitation of freezing. In addition, the transition range (i.e., the range of reduced T in which U is observed) extends from one to several reduced temperature intervals when α exceeds 5. This behavior indicates that the phase transition mechanism at a strongly attractive wall differs from that at a weakly attractive wall. Four possible hypotheses can explain the widening of the transition range: (i) During the freezing process particles are trapped into a kinetically metastable phase without crystallization; hence, no liquid-solid phase transition occurs, and the sudden drop in U disappears. (ii) Particles freeze through a continuous phase transition, similar to the prefreezing phenomenon at the strongly attractive wall [11]
FIG. 1. (Color online) Reduced potential energy U as a function of reduced temperature T for increased SP from α = 1 to 12. The black dashed line is the result of a reference system without the foreign wall. With cooling from T = 0.83 to 0.43, a sudden drop (U 2800) can be observed in U for 1< = α< = 4, manifesting a first-order phase transition. These drops shrink as α continuously increases, which is confirmed by the extension of the temperature range.
and the wetting case predicted by classical nucleation theory [25]. (iii) A cascade of layering first-order phase transitions completes the liquid-solid transition; consequently, the sudden drops are flattened by several drop parts. (iv) The liquid-solid transition involves both continuous and discontinuous phase transition steps that cause the gradual drop. To identify the correct scenario among our hypotheses, we have analyzed the structural properties of the LJ system after the temperature quench, for 1< = α< = 12. By tracking the trajectories from T = 0.83 to 0.43, we observe that the system is always converted into a solid composed of close-packed planes (CPPs). This observation of CPPs is consistent with the simulations of Malley et al. [26]. Moreover, the wall prefers parallel CPPs when SP changes from α = 1 to 2 (Fig. 2); this scenario is specified as “parallel to the wall.” The LJ walls can create a potential well that could trap particles. Thus, a deeper well depth corresponds to more trapped particles. The parallel CPPs [Fig. 2(b)] have more trapped particles than the intersecting state [Fig. 2(a)], which may explain why the
FIG. 2. (Color online) Simulation snapshots at reduced temperature T = 0.43. Yellow (light gray) and green (gray) balls represent fcc and hcp particles, respectively, while liquidlike particles are not shown. The black walls represent the LJ surfaces. The final crystals consist of fcc and hcp form layers, and the layers cross the wall for α = 1 (a) and parallel to the wall for α = 2 (b). The individual images are created using Atomeye [27].
032412-2
WALL-INDUCED PHASE TRANSITION CONTROLLED BY . . .
PHYSICAL REVIEW E 89, 032412 (2014)
stronger attractive wall prefers parallel CPPs. Consequently, a critical value αp between 1 and 2 exists, above which the CPPs are always parallel to the wall [Fig. 2(b)]. Moreover, for α = 1, CPPs with dominant hcp arrays cross the wall at an approximately 62° angle [Fig. 2(a)], and this fixed degree serves as the optimum choice in crystal nucleation and growth of pure hcp structure [28]. We infer that this fixed angle may preferably change to 70.5° when the dominance reverses. This inference is due to the fact that the 70.5° angle is the intrinsic angle consisting of two sets of CPPs in perfect fcc structures [28,29]. Considering that the wall-particle attraction quickly decays to zero as the distance from the wall increases, the phase behaviors of the layers near the wall may essentially differ from those in the bulk. Thus, we mainly examine phase behaviors in
the first three layers. The layering radial distribution function (RDF) [17] gxy (r ) and the bulk RDF g(r) are then investigated. The bulk RDF g(r) is defined as g(r) =
dn(r) , ρ
(3)
where dn(r) presents the average local number at distance r from a fixed particle and ρ is the overall number density. The layering RDF gxy (rxy ) is defined as
i j
i,i=j δ rxy δ rxy − rxy gxy (rxy ) = , (4) ρxy where ρxy is the layering number density and rxy is the distance between the particles in the xy plane. Comparison of the
FIG. 3. (Color online) Layering RDFs. The sixth layer RDFs for α = 2, 4, 5, and 6 are shown in (a), (b), (d), and (f) to represent the bulk situation. (c) and (e) are layering RDFs in the second layer for α = 5 and 6, respectively, and the insets are layering RDFs in the first layer correspondingly. These RDFs all exhibit characteristic features of liquids or solids depending on temperature. In each of them, an abrupt jump occurs from liquid to solid as the temperature reaches the critical value marked in red (gray), indicating the occurrences of a first-order phase transition. When α increases to 6, the transition range of bulk and that of the first two layers are not synchronized. 032412-3
ZHANG, PENG, LONG, ZHOU, LIANG, WAN, ZHENG, AND JU
PHYSICAL REVIEW E 89, 032412 (2014)
FIG. 4. (Color online) RDFs for α = 7 and 8. The insets are parts of the first and third layers. For α = 7, the first two layers are discontinuously converted into solid through a first-order transition (a). For α = 8, a continuous phase transition exclusively occurs in the first two layers (c). The freezing of the sixth and third layers always proceeds through the first-order phase transition (b) and (d). The reduced transition temperatures are marked in red (gray).
RDFs of the sixth layer and those of the bulk reveals that they exhibit the same features (see Supplemental Material [30]). Then, the RDFs of the sixth layer are used to substitute the bulk. In Fig. 3, we observe abrupt jumps in the peaks as the reduced temperature decreases. Before the jump, RDF peaks consist of a prominent peak and a rapidly decaying peak before decaying to unit 1. This behavior exhibits the isotropic nature of liquids. Conversely, when the temperature reaches the critical value (red marks in Fig. 3), the RDFs show a long-range order that is the feature of the crystal. Thus, these jumps are of a first-order phase transition representing liquid-solid conversion. Moreover, the order of the system that has just finished the jump weakens as α increases. Such an order-limited system needs further evolution to reach the final stable state because of the sluggishness [31,32]. When α reaches 6, the freezing of the first two layers and that of the bulk do not synchronize with temperature [Figs. 3(e) and 3(f)], and the so-called “complete wetting” occurs [13,14]. Hence, a critical value αw must exist between 5 and 6. This value distinguishes the phase transition of the first two layers from that of the bulk. Furthermore, the reduced temperatures for the bulk freezing steadily increase and then stabilize at T = 0.67. Such rises in temperature correspond to the shift of transition ranges in the energy profiles (Fig. 1). Treatment of the first two layers as a quasi-2D system, with their freezing anteceding the bulk, is yet to be verified.
Since the particular phase behaviors in the first two layers depend on the SP α, the layering phase behaviors at the more attractive walls are then investigated. Figure 4 shows still remarkable jumps in RDF peaks at α = 7 [Figs. 4(a) and 4(b)]. Nevertheless, these jumps occur within different temperature ranges. In the first two layers, peaks jump from T = 0.79 to 0.75, corresponding to a small drop (U 995) in U as shown in Fig. 1. In the third and in the bulk layers, peaks jump from T = 0.75 to 0.71 and from T = 0.71 to 0.67, respectively. Notably when α reaches 8, these jumps disappear in the first and second layers [Fig. 4(c)], indicating that freezing in the first two layers occurs through a continuous phase transition. We specifically name this scenario as “continuous transition.” However, the third and the bulk RDFs still exhibit the features of the first-order phase transition [Fig. 4(d)]. Therefore, a critical value αc must exist between 7 and 8, beyond which the first two layers continuously freeze. This continuous phase transition differs from the mechanism of phase behavior in three dimensions and can be treated as a quasi-2D phase transition [15]. Page and Sear [11] stated that the energy barrier for prefreezing disappears at the strongly attractive surface, consistent with our observation (continuous transition). The difference is that in our study the continuous transition is specifically limited in the first two layers. Although the wall strength they used is absolutely different from ours, the SP α is really close to ours. This finding
032412-4
WALL-INDUCED PHASE TRANSITION CONTROLLED BY . . .
PHYSICAL REVIEW E 89, 032412 (2014)
FIG. 5. (Color online) NDPs of the first three layers for α = 7 (a) and α = 8 (b). They all exhibit excellent layering at the lower temperature. At each reduced temperature, decreased distance from the wall corresponds to sharper layers. For α = 7, the first three layers all exhibit an abrupt change, and their critical temperatures are T = 0.75, 0.75, and 0.71, respectively. By contrast, for α = 8, continuous changes replace the abrupt states in the first two layers, whereas the third layer retains a discontinuous transition at 0.71.
may imply that the SP α we introduced can be used as a general parameter to describe phase transition at a smooth surface in LJ systems. To further confirm the existence of continuous transition, we examine the number density profile (NDP). In Fig. 5, we can see that peaks become sharper as the reduced temperature decreases. The sharper peaks represent the higher-ordered layers parallel to the wall, and the particles in such layers may approach the 2D situations. Generally, the extending of peaks can be ascribed to thermal fluctuations. Therefore, the wide peaks in the third layer imply that it is relatively far away from the 2D situation. For α = 7 [Fig. 5(a)], the peaks narrow and rise abruptly in the first two layers as the reduced temperature decreases from 0.79 to 0.75. Moreover, these sudden changes are delayed to the next reduced temperature interval (0.75–0.71) in the third layer. By contrast, when α = 8, gradual changes of peaks exclusively occur in the first two layers [Fig. 5(b)], confirming the continuous transition. These observations are consistent with Fig. 4. The mobility of the particles in the first two layers is then considered. The most straightforward method is to estimate the number of particles in each layer. These results are plotted in Fig. 6. The curves undergo a continuous increase at the early stage and then abruptly jump to a constant value for α = 2 and 7 [Figs. 6(b) and 6(c)]. These abrupt increases of the particles in the layers correspond to the occurrences of the first-order phase transitions for which replenishments of particles from the bulk are necessary. Hypothetically, when the attractive force from the wall is sufficiently strong (limiting condition), more particles are confined within a narrow region [inset of Fig. 6(d)] and ought to form a relatively high-ordered array. In this case, as the temperature decreases, the highly ordered layers may solidify without changing their relatively ordered configuration to some extent. In addition, this transition may not need more replenishment from neighboring particles, and a continuous phase transition thus occurs. Actually, curves in α = 8 show a continuously small rise, and are consistent with
this limiting assumption [Fig. 6(d)]. However, the relatively huge and continuous increases for α = 1 [Fig. 6(a)] do not represent a continuous phase transition. Evidently, such increases should be attributed to the freezing with nonparallel layers near the wall [Fig. 2(a) and inset of Fig. 6(a)]. Finally, reverse procedures (melting) are performed using the final freezing results. We select the midpoint (Ttm ) of transition ranges for both freezing and melting processes because of the phase transitions across a relatively wide temperature range. As shown in Fig. 7, three regions (each later region is included in the former) marked in green, gray, and yellow correspond to the scenarios when SP (α) reaches the critical value αp , αw , and αc , respectively. We observe the considerable hysteresis between freezing and melting procedures. Such hysteresis diminishes as α increases, which represents the reduction of the energy barrier to the bulk freezing and the weakening of the first-order transition. Nevertheless, this diminution terminates when α reaches αc , which signifies that the bulk nucleation barrier cannot be thoroughly abolished [13]. In addition, phase behaviors above αc (yellow region) appear to be similar because Ttm hardly varies with α. In this region, α is sufficiently large to make the first two layers behave as a quasi-2D system. Such quasi-2D system possesses a denser and more ordered array of particles [Fig. 5(d) and inset of Fig. 6(d)], which ensures their lower mobility in the out-of-plane direction. In addition, the wall-particle interaction decays to a minute value in the bulk. Therefore, bulk freezing is reasonably expected to proceed in an isolated manner even through it starts from the wall. Several runs are performed for both the freezing and the melting processes. Data for freezing are coincident, whereas those for melting exhibit large fluctuations. These phenomena can be due to the fact that freezing always starts from the surface [insets of Figs. 6(a) and 6(d)] but melting begins from the exterior where fcc and hcp CPPs array with distinct sequence in each run. We therefore conclude that freezing depends on the walls, whereas melting depends on the
032412-5
ZHANG, PENG, LONG, ZHOU, LIANG, WAN, ZHENG, AND JU
PHYSICAL REVIEW E 89, 032412 (2014)
FIG. 6. (Color online) Number of particles in the first two layers. Arrows represent the transition tendency (vertical for abrupt increase and oblique for gradual increase). For abrupt changes, we mark the critical reduced temperature in red (gray). The insets are the snapshots of early stages for α = 1 and 8, and the details are the same as those in Fig. 2.
structures. This conclusion can be supported by Fig. 7 in which freezing is sensitive but the melting is insensitive to α. To test the finite size effects, a larger simulation box (Lx = Ly = Lz = 29.37σ ) is used. Results show few finite size effects (Fig. 7). Recalling the previous hypotheses, the first and the second are excluded, whereas the third and the fourth characterize the true
FIG. 7. (Color online) Midpoint of transition range as a function of SP α. Freezing and melting procedures are shown. The black dashed line denotes the triple point. Three regions are shown in green (light gray), gray (gray), and yellow (white) for the three scenarios described as parallel to the wall, complete wetting, and continuous transition, respectively. The open marks represent the results of the larger systems used to investigate the finite size effects.
scenarios of phase transitions at strongly attractive surfaces in our systems. IV. CONCLUSIONS
We have performed computer simulations to study the liquid-solid transition of a LJ system in the presence of an attractive surface. We identify different scenarios for phase transitions depending on the value of the parameter α, i.e., the ratio between the wall-particle and the particle-particle attraction strength. For weakly attractive walls, crystallization occurs at the walls with nonparallel CPPs. As the wall strength increases to the first critical value αp , “parallel to the wall” emerges. Throughout the simulation box crystallization occurs at the same temperature via a first-order transition when α < αw . When α = αw , the first two layers near the wall exhibit advanced freezing, and complete wetting occurs. Furthermore, for α>αc , the first two layers undergo freezing not through first-order phase transition, but through a continuous phase transition, which is also known as continuous transition This phenomenon has been observed to take place always in the first two layers, while in the bulk freezing still occurs via a first-order transition. Hence, the two layers can be treated as a quasi-2D system. Their liquid-solid transition may continuously occur through a hexatic phase, because Wierschem and Manousakis [33] stated that phase transition in slit pore with a soft core potential is expected to proceed through the KTHNY mechanism. However, we cannot determine whether the KTHNY mechanism can explain our results because distinguishing a hexatic phase from others requires a very large
032412-6
WALL-INDUCED PHASE TRANSITION CONTROLLED BY . . .
system size and a very long time, let alone the imperfection of the current techniques [33]. In the third layer, liquid-solid transition always proceeds through first-order phase transition and can occur earlier than the bulk freezing. This finding implies that the third layer may serve as the transient regime between the quasi-2D layers and the bulk phase. We also find that a considerable energy barrier to the bulk freezing always exists because the hysteresis cannot be completely abolished.
[1] P. R. ten Wolde and D. Frenkel, Science 277, 1975 (1997). [2] R. P. Sear, J. Chem. Phys. 129, 164510 (2008). [3] T. H. Zhang and X. Y. Liu, Angew. Chem., Int. Ed. 48, 1308 (2009). [4] A. Cacciuto, S. Auer, and D. Frenkel, Nature 428, 404 (2004). [5] M. Heni and H. L¨owen, Phys. Rev. Lett. 85, 3668 (2000). [6] W. S. Xu, Z. Y. Sun, and L. J. An, J. Chem. Phys. 132, 144506 (2010). [7] L. M. Ghiringhelli, C. Valeriani, E. J. Meijer, and D. Frenkel, Phys. Rev. Lett. 99, 055702 (2007). [8] P. R. ten Wolde, M. J. Ruiz-Montero, and D. Frenkel, J. Chem. Phys. 104, 9932 (1996). [9] R. Radhakrishnan, K. E. Gubbins, and M. SliwinskaBartkowiak, J. Chem. Phys. 116, 1147 (2002). [10] C. Alba-Simionesco, B. Coasne, G. Dosseh, G. Dudziak, K. E. Gubbins, R. Radhakrishnan, and M. Sliwinska-Bartkowiak, J. Phys.: Condens. Matter 18, R15 (2006). [11] A. J. Page and R. P. Sear, Phys. Rev. E 80, 031605 (2009). [12] N. Gribova, A. Arnold, T. Schilling, and C. Holm, J. Chem. Phys. 135, 054514 (2011). [13] S. Auer and D. Frenkel, Phys. Rev. Lett. 91, 015703 (2003). [14] M. Dijkstra, Phys. Rev. Lett. 93, 108303 (2004). [15] K. J. Strandburg, Rev. Mod. Phys. 60, 161 (1988). [16] R. Radhakrishnan, K. E. Gubbins, and M. SliwinskaBartkowiak, Phys. Rev. Lett. 89, 076101 (2002). [17] A. H. Marcus and S. A. Rice, Phys. Rev. E 55, 637 (1997). [18] J. G. Dash, Rev. Mod. Phys. 71, 1737 (1999). [19] D. Koga, H. Tanaka, and X. C. Zeng, Nature 408, 564 (2000). [20] C. Y. Lee, J. A. McCammon, and P. J. Rossky, J. Chem. Phys. 80, 4448 (1984). [21] P. Kumar, S. V. Buldyrev, F. W. Starr, N. Giovambattista, and H. E. Stanley, Phys. Rev. E 72, 051503 (2005). [22] In the homogeneous system, some characteristic nuclei could form from a homogeneous matrix. To ensure that a stand-alone
PHYSICAL REVIEW E 89, 032412 (2014) ACKNOWLEDGMENTS
This work was supported by Institute of Nuclear Physics and Chemistry (INPC), China Academy of Engineering Physics. We thank Nadezhda Gribova and Li Mao for the fruitful discussions about the layering radial distribution function. We also thank Jing Xu from Peking University for assistance with the manuscript.
[23] [24] [25] [26] [27] [28]
[29] [30]
[31] [32]
[33]
032412-7
nucleus matures without spatial restriction in most directions, we choose 3400 particles in a cubic box (with length L = 17.62). Such unimpeded nucleation and growth can provide us its formation mechanism at the atomic level. J. D. Honeycutt and H. C. Andersen, J. Phys. Chem. 91, 4950 (1987). S. Plimpton, J. Comput. Phys. 117, 1 (1995). R. P. Sear, J. Phys.: Condens. Matter 19, 033101 (2007). B. O’Malley and I. Snook, Phys. Rev. Lett. 90, 085702 (2003). Atomeye: Atomistic Configuration Viewer, http://li.mit.edu/ Archive/Graphics/A/. A. J. Page and R. P. Sear, J. Am. Chem. Soc. 131, 17550 (2009). Four different sets of close-packed planes (CPPs) exist in the fcc structure, but only one exists in the hcp state. We find that two sets of such CPPs in perfect fcc structures mutually intersect with an approximate angle of 70.53°, whereas the CPPs and the {1 − 10 − 1} planes in the hcp structures provide an angle of 62.06°. Such two angles are consistent with the work of Page and Sear [11]. Hence, we infer that these two intrinsic angles have a significant function in crystal nucleation and growth, which can help fcc (with 70.53° ridge) and hcp (with 62.06° ridge) blocks to form and grow without restriction. B. W. van de Waal, Phys. Rev. Lett. 67, 3263 (1991). See Supplemental Material at http://link.aps.org/supplemental/ 10.1103/PhysRevE.89.032412 for comparison of RDFs of the sixth layer and those of the bulk. W. A. Caldwell, J. H. Nguyen, B. G. Pfrommer, F. Mauri, S. G. Louie, and R. Jeanloz, Science 277, 930 (1997). H. Cynn, C. S. Yoo, B. Baer, V. Iota-Herbei, A. K. McMahan, M. Nicol, and S. Carlson, Phys. Rev. Lett. 86, 4552 (2001). K. Wierschem and E. Manousakis, Phys. Rev. B 83, 214108 (2011).
