Size Effect on Nucleation Rate for Homogeneous Crystallization of

Article pubs.acs.org/JPCB

Size Effect on Nucleation Rate for Homogeneous Crystallization of Nanoscale Water Film Yongjun Lü,*,† Xiangxiong Zhang,‡ and Min Chen*,‡ †

School of Physics, Beijing Institute of Technology, Beijing 100081, P. R. China Department of Engineering Mechanics, Tsinghua University, Beijing 100084, P. R. China

‡

S Supporting Information *

ABSTRACT: The nucleation rate from classical nucleation theory is independent of sample size. In the past decades, several experimental and theoretical studies argued that the homogeneous nucleation rate of ice in supercooled droplets increases when the drop size is decreased. In this paper, we investigate the nucleation of ice in nanoscale water films using molecular dynamics simulations. We found that the nucleation rate of ice actually decreases when the film thickness decreases in the nanoscale regime. A theoretical model is presented to interpret the mechanism of nucleation rate decrease, which agrees well with the simulation results. The model divides films into the near-surface and the middle regions that are characterized by relatively low and high nucleation rates, respectively. The middle region dominates the nucleation process of films, whereas its effect is continuously weakened when increasing volume fraction of the near-surface region by decreasing the film size, leading to a decrease of the total nucleation rate. The structural and thermodynamic analyses indicate that the high stress induced by the surface layering slows down the diffusion and increases the nucleation barrier in the near-surface region, which is responsible for the low nucleation rate and eventually the decrease of the total nucleation rate. defined as surface-dominated processes.12 The surface of water drops is believed to facilitate nucleation when the drop size is smaller or the surface-to-volume ratio is larger than a certain threshold value. Tabazadeh et al. proposed the thermodynamic criterion that nucleation is favored by the surface when σvs − σvl < σls, where σvs, σvl, and σls denote the surface tensions of the vapor/solid, vapor/liquid, and liquid/solid interfaces, respectively. This indicates that the decrease of nucleation barrier in surface-dominated nucleation can be attributed to the partial wetting of crystalline nuclei with liquid. From Tabazadeh’s work, surface-dominated nucleation has attracted considerable interest and has much controversy. Some experimental measurements supported their assumption,13,14 while other experiments did not observe a difference in nucleation rate for water drops with different sizes.15,16 However, the sizes of water drops in experiments are mostly in micrometer scale. When the sample size is reduced to the nanoscale, it is expected that the surface or interface effect becomes more prominent, and its role in homogeneous nucleation may be important. Experimental methods for the measurement of nucleation rates mostly rely on measuring the fraction of unfrozen droplets in massive water drops as a function of observation time. The advantage of this method is the direct output of nucleation rate, whereas less information on the nucleation mechanism is obtained. In contrast, computer simulations are useful to

1. INTRODUCTION Crystallization of water is an important phenomenon involved in many fields from chemistry to biology and atmospheric science. When cooling from 273 K, water maintains a metastable liquid state named supercooled water. In the highly supercooled regime, it is hypothesized that there exists a liquid−liquid phase transition (LLPT) at 228 K.1,2 Approaching the transition from above, the thermodynamic properties of liquid water show anomalous fluctuations. Such behavior readily triggers crystallization prior to the potential LLPT, leading to the absence of supercooled water below 228 K, known as “no man’s land”.3−5 It is generally recognized that homogeneous nucleation of ice occurs at about 235 K, which is close to the liquid−liquid transition point. Homogeneous nucleation of ice can be achieved in the laboratory by decreasing drop size and increasing cooling rate and is widely observed in cirrus clouds.6−9 Nucleation rate is defined as the probability of forming a nucleus in per unit time and unit volume and is related to the free energy barrier of forming the critical nucleus. It does not depend on the sample size in classical nucleation theory (CNT).10,11 This conclusion, however, does not well account for the experimental observations of the homogeneous nucleation of ice in water drops. Comparison of the available experimental nucleation rates shows that they vary with the experimental conditions and differ by up to several orders of magnitude, which exceeds the range of measurement errors.12 Through examining the experimental data, Tabazadeh et al. found that some homogeneous nucleation events can be © 2013 American Chemical Society

Received: May 3, 2013 Revised: August 7, 2013 Published: August 12, 2013 10241

dx.doi.org/10.1021/jp404403k | J. Phys. Chem. B 2013, 117, 10241−10249

The Journal of Physical Chemistry B

Article

Figure 1. Schematic illustration of simulation cells.

maximum at about 202 K.25 It is indicated that water crystallization transforms from nucleation-dominated to growth-dominated at about 202 K from above. The temperature of 205 K was regarded as the lowest value for nucleationdominated crystallization. Second, the crystallization time is expected to greatly decrease close to 202 K (in tens of nanoseconds), which is beneficial for more effective sampling using the MFPT method. In the simulations, the x- and y-axis dimensions were kept constant at Lx = Ly = 5.42 nm, and the zaxis dimension (Lz) was varied to produce water films with different surface-to-volume ratios. The simulations were performed for nine groups with Lz ranging from 4.72 to 20.72 nm corresponding to particle numbers from 4598 to 20 086. All simulations were carried out in the NVT ensemble. The temperature was controlled with the Nosé and Hoover thermostat.26,27 Periodic boundary conditions were applied in the three directions, and a time step of 10 fs was used. Configurations are sampled every 1000 time steps for each simulation. The size of the largest cluster in each configuration (n) and the time of its first appearance (τ(n)) were recorded. More than 150 independent simulations were performed with different initial configurations for averaging τ(n). The nucleation rate was directly calculated by fitting the time of a growing nucleus versus the nucleus size from the MFPT method. As well as nucleation rate, critical nuclei size, nucleation barrier, and attachment rates of the largest clusters were also calculated (see Supporting Information for details). The locally ordered clusters were identified using the CHILL algorithm, which is special for the arrangement of the orientation of local structure and has been shown to perform well in distinguishing between hexagonal ice, cubic ice, and liquid water.28 In addition, this algorithm distinguishes an ordered intermediate phase that bridges liquid and perfect ice (see Supporting Information). The intermediate ice phase is usually observed in the region near the surface or interface and, thus, is helpful for detecting potential nucleation initiated from the surface.

investigate the thermodynamic characteristics of nucleation, for example, umbrella sampling and transition path sampling.17−20 These methods focus on the evaluation of the nucleation barrier and free-energy landscape by intensive biased sampling and then predict the nucleation rate from CNT. Alternatively, a more direct method called the mean first-passage time (MFPT) was developed by Wedekind et al.21,22 This method is straightforward to give nucleation rate and allows the prediction of the nucleation barrier and attachment rate for forming a critical nucleus. However, MFPT requires nucleation events to be produced in a single simulation, which is actually rather difficult for popular water models such as SPC/E and TIP4P.23 Recently, Molinero et al. proposed a coarse-grained water model (mW water). Although this model underestimates the melting enthalpy of ice and overestimates the ice density at the melting temperature, it reproduces the liquid density, surface tension, and melting point as well as the structure of liquid water with better accuracy than other popular water models.24,25 Importantly, it can simulate the crystallization of ice in a relatively short computation time, which is valuable for investigating the nucleation rate of ice using MFPT. In this paper, we have calculated the nucleation rate for homogeneous crystallization of water films with thicknesses ranging from 4 to 20 nm using molecular dynamics (MD) simulations. We find that the nucleation rate decreases as the film thickness is reduced. We propose a theoretical model that divides the film into two different regions in nucleation, the near-surface and the middle regions. Due to a large nucleation barrier induced by the density oscillation near the surface, the near-surface region bears the lower probability of nucleation in contrast to the middle part. With decreasing film size, the volume fraction of the near-surface region exponentially enlarges, leading to the decrease of the total nucleation rate. The theoretical prediction is in good agreement with the simulation results.

2. COMPUTATIONAL METHODS For the MD simulations, we used the mW model to investigate nucleation of ice in water films. The water films were prepared by expanding the simulation box of bulk systems that were equilibrated for 1.2 ns at 300 K by three times along the z-axis. Figure 1 shows the schematic diagram of the simulation cell. The vapor/liquid/vapor systems were then relaxed at 205 K until the crystallization was finished. There are two reasons for choosing this crystallization temperature. First, MD simulations have shown that the nucleation rate of mW water reaches a

3. RESULTS AND DISCUSSION 3.1. Nucleation Rate of Ice in Nanosized Water Film. By fitting the size of the largest clusters as a function of time of their first appearance (see Figure S1, Supporting Information), we directly obtained the probability of nuclei in a sample and in unit time R. Figure 2a shows that the probability of nuclei per unit time in the sample decreases with decreasing film thickness 10242



Article

R = f (n*)Z

exp[−β ΔG(n*)] n = n*

∫n = 0 dn exp[−β ΔG(n)]

(1)

where n* is the size of the critical nucleus and f(n*) and ΔG(n*) are the attachment rate and nucleation barrier of the critical nucleus, respectively. Z is the Zeldovich factor and β = 1/kBT. From eq 1, a decrease of nucleation rate implies a slower attachment rate or higher nucleation barrier. To clarify this, the free energy required to form an embryo of size n (ΔG(n)) was calculated using R and τ(n) (Figure S2, Supporting Information), and then, the nucleation barrier ΔG(n*) was determined from the maximum of ΔG(n), as shown in Figure 3. The nucleation barrier shows a considerable increase from

Figure 2. Nucleation rate of water films as a function of film thickness. (a) Probability of nuclei in a sample and in unit time (R). Assuming homogeneous nucleation occurs in the film, the volume nucleation rates (R/V) are plotted in (b).

Figure 3. Calculated nucleation barrier as a function of the film thickness from the MFPT method. The solid line is a polynomial fit.

2.19 kBT to 5.19 kBT accompanied by a decrease of film thickness. This is the main reason for the decrease of the nucleation rate. Compared with the values from experiments, the present nucleation barriers are an order of magnitude smaller. It is reasonable considering that the simulation temperature was controlled just above the threshold from the nucleation- to the growth-dominated crystallization. Below the threshold, the nucleation barrier disappears. The low nucleation barrier explains why the calculated nucleation rates are larger than the experimental measurements by orders of magnitude, so that nucleation events are observed within the time interval of several nanoseconds in the simulations. Another possible factor responsible for the decrease of the nucleation rate is the attachment rate of critical nuclei. This quantity can be obtained from an attachment rate versus embryo size plot using the MFPT method and can also be directly calculated using the model proposed by Auer and Frenkel as19,29 1 f (n*) = ⟨[n(t ) − n(0)]2 ⟩ (2) 2t

(Lz). In our work, all of the simulations were performed with a constant surface area. Assuming that the nucleation process is surface-dominated, the probability of nuclei R should be independent of the value of Lz. If it is volume dominated, R should vary linearly with volume. The simulations do not agree with the former but look like volume-dominated nucleation following the linear relation, R = −0.0186 + 0.0076LZ. If volume-dominated nucleation dominates the nucleation process, the corresponding nucleation rate, JV = R/V, should be independent of film thickness. In fact, this is not the case (Figure 2b). When Lz > 12.7 nm, the nucleation is dominated by volume, as is indicated by the approximate thickness independence of JV. When Lz < 12.7 nm, nucleation deviates from volume-dominated nucleation and the nucleation rate drops by more than twice for the thinnest film. This clearly demonstrates the size dependence of the nucleation rate at the nanoscale. Tabazadeh’s model based on the experimental results of micrometer water droplets predicts that the total nucleation rate is jointly determined by volume- and surfacedominated nucleation rates, with surface-dominated nucleation having a higher probability than volume-dominated nucleation due to the low nucleation barrier.12 Therefore, high nucleation rates are expected in small droplets with large surface-tovolume ratios. However, this conclusion fails to account for the present work when the sample size is decreased to the nanoscale. In terms of dynamics, nucleation is an activated process, depending on the activation barrier and diffusion. The MFPT method derives the probability of nuclei in a sample and in unit time R as the following expression22

The results produced by the two methods are in good agreement (Figures S3 and S4, Supporting Information). Here, we employed the latter because it allows us to distinguish the contributions from the near-surface and the middle regions. The calculated attachment rates of critical nuclei do not show a significant dependence on the film thickness. 3.2. Role of the Surface in Nucleation. Because of constant surface area in the simulations, decreasing film thickness means that the surface-to-volume ratio increases, and the influence of the surface may become important. In this section, we investigate the role of the surface in nucleation for films of different thicknesses. 10243



Article

Lennard-Jones liquids, and alloys.33−35 Recent MD simulations confirmed that crystalline silicon preferentially forms within low-density liquid droplets.36 This behavior results from the low barrier of the intermediate phase in the highly supercooled regime, which can be explained by the well-known Ostwald’s step rule.37 For crystallization of water, I molecules play a crucial role in stimulating nucleation, and their amount and distribution would affect the nucleation rate. Therefore, we can scale the probability of forming critical nuclei by analyzing the characteristic features of clusters containing I molecules. We averaged the number of neighboring molecules belonging to I and crystal phases for a molecule I before nucleation, which is helpful to determine the size of I molecule-centered clusters. The calculated average value was 1.62, which means that on average 2−3 I molecules aggregate together. This value does not significantly vary with decreasing film thickness, suggesting that the thickness of the film does not influence the energy landscape for forming the I phase. The density profile of I molecules along the z direction, however, shows a strong size dependence. Figure 5a shows an example of the water film with Lz = 4.72 nm compared to that of Lz = 20.72 nm. We also plotted their total density profiles along the z direction and then defined a density surface by fitting the density profiles with a hyperbolic tangent function. We define an I surface in a similar way to the density surface. The two surfaces are found

From Tabazadeh’s model, when surface-dominated nucleation is favored thermodynamically over volume dominated, the nuclei partially wet its melt. We examined all of the nucleation events obtained in our simulations and did not observe this phenomenon, which may explain why the nucleation rate decreases rather than increases with increasing surface-tovolume ratio. We averaged the location of the critical nuclei along the z direction, and the fraction profile is shown in Figure 4. The critical nuclei that are formed first are less frequent in

Figure 4. Distributions of nucleation events along the z direction. The film was divided into 5−6 slices, and the nucleation events were counted in each slice over all the trial simulations. z* is z normalized by the z-axis dimension (Lz).

the two regions close to surfaces, and most of the nuclei appear inside the water films. The profiles do not show an obvious dependence on film thickness. Figure 4 indicates that nucleation of ice in the nanoscale water films is not homogeneous: the regions near the two surfaces have a low probability of nucleation, and the middle region has a high probability. We used the CHILL algorithm to trace the structural evolution during nucleation.28 The size of the critical nuclei that is obtained by fitting τ(n) using the MFPT method is independent of the film size with a value of 48 ± 3 (Figure S5, Supporting Information). The decrease of nucleation rate in the thinner films does not result in larger critical nuclei. We have examined the structural composition of the critical nuclei and found that hexagonal and cubic ice coexist in the nuclei. The fraction of the two structures varies with the nucleus location. For critical nuclei formed in the middle region, cubic ice is preferred, whereas the nuclei are dominated by hexagonal ice near the surfaces. A hybrid of the two structures in crystallites has been reported in crystallization of bulk mW water and those confined in nanopore by MD simulations.28,30 X-ray and neutron diffraction patterns also verified that the crystallized ice from supercooled water was composed of stacked cubic and hexagonal structures.31,32 Here, we jointly identify them as crystal phases instead of distinguishing hexagonal and cubic ice. A common feature of the critical nuclei is the abundance of molecules belonging to the intermediate phase (hereafter referred to as I molecules). The I molecules are distributed on the interface between crystalline nuclei and liquid. In fact, I molecules are abundant prior to the formation of critical nuclei. In most cases, they aggregate into small clusters and are dispersed in the liquid. Once large clusters containing I molecules form, some molecules inside the cluster are identified as crystalline, which are likely to develop into critical nuclei. Similar precrystallization has been observed for colloids,

Figure 5. (a) Profiles of intermediate phase (I) atom fraction for Lz = 4.72 and 20.72 nm. The density profiles are also shown in (a). The profiles of I atoms show the characteristic features of low nucleation rate in the near-surface region. (b) I surface plotted against film thickness, which shows an inverse exponential relationship with film thickness. The dashed lines in (a) indicate the near-surface regions. 10244



Article

requires the near-surface region to be defined in advance. The definition of the near-surface region should correspond to the low nucleation rate region. As discussed above, the distribution of I molecules indicates the nucleation ability; therefore, we use the I surface to approximate the near-surface region. Correspondingly, the independent variable x in eq 5 is equivalent to the I surface width, η, and other quantities associated with the near-surface region are referred to I surface region. From CNT, κ = ρlZf, where f is the attachment rate as mentioned above and ρl is the number density of the liquid. The variable, R/V, in eq 5 defines an average nucleation rate. The Zeldovich factor (Z) and the attachment rate (f) can be calculated using the MFPT method. However, eq 7 does not distinguish the contributions from the middle and the nearsurface regions. The Zeldovich factor is defined as

to be different. Relative to the density surfaces, the positions of the I surfaces shift toward the interior, which becomes more pronounced for the thin films. This indicates that the region near the film surface has a low concentration of I molecules. Moreover, the range considerably increases with decreasing film thickness. We define the quantity η as the ratio of the depth of the I surface phase to the film thickness. Figure 5b shows that the values of η exponentially increase with decreasing film size and approach 0.5. Because the number of I molecules partially measure the nucleation probability, we speculate that the nucleation rate in the near-surface regions is less than that in the interior. Furthermore, the nucleation rates of the films with small thickness decrease due to the large values of η, which is consistent with the calculated nucleation rates as well as the distributions of nucleation events discussed above. On the basis of these results, we proposed a nucleation mechanism of the size dependence. For large systems, the nearsurface region associated with a low nucleation rate is small relative to the total volume, so its influence on the total nucleation rate is small. In this case, nucleation is characterized as volume dominated. With decreasing film thickness, suppression of nucleation originating from the near-surface region becomes more prominent owing to large η, resulting in the decrease of the total nucleation rate. In the following section, we provide a theoretical model to account for the nucleation behavior from thermodynamics and dynamics. 3.3. Theoretical Model for Nucleation of Nanoscale Water Film. Both the distributions of critical nuclei and I molecules suggest that the near-surface regions have a low probability of nucleation. To describe the nucleation behavior, we divide the water film into two regions: the middle and the near-surface regions. We assume that nucleation in both the two regions follows the bulk behavior described by CNT but differ in their nucleation rate. Therefore, the probability of nuclei per sample in unit time, R, is written as the sum of two terms, R = JV V1 + JS V2

Z=

(8)

In our simulations, the size of the critical nucleus (n*) does not depend on the system size. This suggests that the free energy landscape near n* does not significantly change, which is also supported by the approximately constant values of Z for different systems (Figure S6, Supporting Information). Thus, the difference of the partial Z between the middle and the nearsurface regions can be ignored. The attachment rate is related to the self-diffusion coefficient (Ds), the size of the critical nucleus (n*), and the jump distance of a molecule to the critical nucleus (λ) by f = 24DSn*2/3/λ2. The self-diffusion coefficients near the surface usually exhibit different behavior from that in the bulk, and the jump distance is difficult to directly measure. A feasible route to determine f V and f S is using eq 2. To do this, we distinguish the samples that nucleate preferentially in the near-surface region from those in the middle regions for all of the trial simulations and then trace the growth of their embryos into the critical nuclei. The attachment rates f V and f S are calculated as shown in Figure 6. The values of f S are smaller than those of f V by 10% to 50%. Another quantity that influences the attachment rate is the number density of the liquid. We note that the density profile near the two surfaces exhibits layering behaviors as shown in Figure 5a. The density oscillation rapidly dissipates from the surface and can be detected three layers below the surface up to a depth of about 8.8 Å. These feature characteristics do not

(3)

where JV and JS are the nucleation rates in the middle and the near-surface regions. V1 and V2 are the volumes of the two regions, where V2 sums over the two near-surface regions of the film, and V = V1 + V2. The probability of forming a critical nucleus depends exponentially on the free-energy barrier, and the nucleation rate is given by ⎛ ΔGV ,S ⎞ JV ,S = κV ,S exp⎜ − ⎟ ⎝ kBT ⎠

|ΔG″(n*)| 2πkBT

(4)

where the subscripts V and S indicate the middle and the nearsurface regions and κ is the kinetic factor. Combining eq 3 and eq 4 gives R = JV [1 + (C − 1)x] (5) V where ⎛ ΔU ⎞ C = α exp⎜ − ⎟ ⎝ kBT ⎠

α=

(6)

κS κV

(7)

ΔU = ΔGS − ΔGV and x = where is the thickness of the near-surface region. Evaluation of nucleation rate using eq 5 2LSZ/LZ,

LSZ

Figure 6. Attachment rates of the critical nuclei in the middle and the near-surface regions for different sizes of the near-surface region. 10245



Article

surface is negligible. With decreasing film thickness, the influence of the near-surface region on nucleation becomes important. In contrast to the description of Tabazadeh et al., the near-surface region is found to have a low nucleation rate. The total nucleation rate of water films includes two contributions: from the middle and the near-surface regions. The width of the near-surface region does not change with decreasing film thickness. This means that the volume fraction of the near-surface region increases with decreasing film thickness, and its effect on nucleation becomes more important, resulting in a decrease of the total nucleation rate when the thickness of films is decreased. The nucleation process in the near-surface region still follows the description of CNT, and the low nucleation rate is essentially induced by a smaller attachment rate and a higher nucleation barrier, with the latter being dominant. 3.4. Nucleation Mechanism in the Near-Surface Region. From the discussion above, the decrease of attachment rate and the increase of nucleation barrier associated with the near-surface region jointly lower the total nucleation rate as the film thickness is decreased. The attachment rate is related to the self-diffusion coefficient and the jump distance from CNT. Numerical results for L−J liquids show the same order for the jump distance and the interparticle spacing.19 Suppose this conclusion also holds for mW water, then the difference between the near-surface region and the middle region is small due to the small density difference between them. The decrease of the attachment rate, therefore, mainly arises from the decrease of self-diffusion. We calculated the mean squared displacement (MSD), ⟨|r(t) − r(0)|2⟩, of the water molecules in the films. To obtain the MSD of the near-surface region, we sampled the water molecules that are detained in the nearsurface regions in the time interval of 30 ps. In general, the diffusion would display anisotropic behaviors in the vicinity of liquid−vapor surface. For the water films in this work, the anisotropic diffusion holds in the entire samples. Figure 8 shows the lateral and the normal MSD in the near-surface and the middle regions. The motion of molecules in the direction parallel to the surface is faster than that in the normal direction

depend on the film thickness. The oscillation magnitude of the outermost layer is about 4−5% greater than that of the bulk density, and the average density in the near-surface region is only ∼1−2% greater than the bulk value. Therefore, the influence of the number density can be ignored. Combining these results, the value of α is between 0.5 and 0.9. The next quantity to be determined is the difference in nucleation barrier (ΔU). The nucleation barrier shown in Figure 3 has a negative dependence on the film thickness. The increase of nucleation barrier essentially relates to the increasing fraction of the near-surface region. Therefore, it is reasonable to assume that nucleation barrier depends on η. Plotting the nucleation barrier as a function of η shows an approximately linear relationship. From the definition of η, no near-surface region forms when η = 0. We attempted to further increase the surface-to-volume ratio to 0.58 for nucleation. Unfortunately, we did not observe nucleation, even when the simulation time exceeded 500 ns. In this regard, it is thought that nucleation has been dominated by the near-surface region when η = 0.5. We linearly extrapolated ΔG to η = 0 and 0.5 to approximate ΔGV and ΔGS, which results in a ΔU of 3.23 kBT. Then, the calculated value of C is between 0.02 and 0.03, which is much smaller than 1, and is thus ignored in eq 5. As a result, eq 5 can be simplified as R /V ≈ JV (1 − η)

(9)

Before estimating the average nucleation rate using the above equation, the nucleation rate without the effect of the nearsurface region JV is required. Following a similar procedure to calculating ΔGV, we replotted R/V as a function of η and extrapolated it to η = 0. The corresponding value (2.70 × 1032 m−3s−1) is approximately JV. To test the extrapolation, we calculated the nucleation rate of bulk water at 205 K. The system contains 4598 mW molecules, and the simulation details are the same as the film cases. The obtained value ((2.73 ± 0.20) × 1032 m−3s−1) is consistent with the extrapolated result. Figure 7 shows the comparison of the results obtained directly from the MFPT method and using eq 9. The theoretical prediction shows good agreement with the MFPT results, which validates our assumption for the nucleation of ice in nanoscale water films. Here, we give a brief overview of its nucleation route. For large-sized water films, nucleation occurs randomly in films as with the bulk ice, as the influence of

Figure 8. Lateral and normal mean-squared displacements (MSD) of atoms in the near-surface (open circles) and the middle (solid lines) regions. (a) Lz = 4.72. (b) Lz = 20.72 nm. The dashed lines are linear fits to the MSD curves.

Figure 7. Volume nucleation rates as a function of η, where η is the normalized total thickness of the near-surface regions. The solid line is the theoretical prediction from the model proposed in this work. 10246



Article

throughout the films. In both the lateral and the normal directions, the molecules in the near-surface region have lower mobility than those in the middle region, which qualitatively agrees with the speculation from the attachment rate. Recent MD simulations reported a decrease of diffusion in confined liquids, in which diffusion was described as an activated process.38 It was shown that the confinement decreases the potential energy per molecule and leads to an increase of the activation barrier for diffusion. However, this mechanism cannot describe our results because the average potential energy per molecule actually increases as the film thickness is decreased. In fact, for droplets or liquid slabs, a negative lateral stress (pT) region is expected in the vicinity of the liquid−vapor surface due to surface tension. Figure 9 shows the lateral and

middle regions, and no nuclei are able to develop in the vicinity of surfaces. We note recent MD simulations that showed the similar decrease of nucleation rates with decreasing size in supercooled water nanodroplets.48 The suppression of ice crystallization was attributed to the Laplace pressure that depends on the curvature of droplets. From the thermodynamic model proposed by Li et al., in the presence of a pressure p, the difference in chemical potential between liquid and solid (Δμsl) should involve the contribution of the pressure and the molar volume difference between liquid and solid (ρΔν).48 This modification leads to an exponential variation of the nucleation rate with ρΔν. For tetrahedral liquids such as water and liquid silicon, Δν < 0 because of a negative slope of the liquid-ice coexistence curves. Therefore, the Laplace pressure (ρ > 0) is expected to decrease the nucleation rate. For the planar surface, the Laplace pressure becomes zero and its effect on nucleation rate vanishes. It is inferred that the negative lateral pressure near the planar surface induced by surface tension would increase ice nucleation rates according to Li’s model.48 However, our simulations show that the crystallization is not favored in the near-surface region of water films. One explanation for this contradiction is that the volume difference may be positive in our case. The MD simulations by Li et al. have shown that the density of ice is higher than that of liquid water at 205 K for the mW water model.48 Therefore, the negative lateral pressure is expected to decrease the nucleation rate because ρΔν < 0. Another reason is that the surface density layering is against the effect of negative pressure. The density layering near the surface was found to induce the fluctuations of the potential energy and the force on atoms.47 The surface layering actually exerts an additional normal pressure on the films, which is analogous to the Laplace pressure and is expected to suppress nucleation. We then evaluated the effect of the additional pressure on the nucleation barrier as follows. On the basis of CNT, the free energy of a spherical nucleus of radius r has the following form:

Figure 9. Lateral and normal pressure tensor profiles.

the normal stress (pN) distributions in the water films. The calculated lateral stress near the two surfaces reaches about −230 MPa. For simple liquids, the increase of pressure would lead to a decrease of the diffusion coefficient and fails to describe the diffusion behavior of liquid water in the low temperature regime, which has negative expansivity. Within the low temperature regime, the hydrogen bonded network distorts with increasing pressure, resulting in an increase of the diffusion coefficient until reaching a maximum.39,40 This suggests that the negative lateral pressure near the surface is responsible for the low lateral diffusion of the near-surface region. The normal stress shows fluctuation across the surfaces, which does not correctly describe the realistic system because of the violation of mechanical equilibrium, and results from the computational method used in this work (see Supporting Information).41,42 The normal stress distribution may be more complex than the expectation because our results indicate that the near-surface region is characterized by significant density layering. A similar phenomenon has been observed in several low melting point metal and alloy drops, including Hg, Ga, Sn, and Au−Si.43−46 MD simulations have shown that the motion perpendicular to the surface of these metallic drops displays a similar layering.47 The mean residence time in the outermost layer is much longer than the values in the inner ones. The directional motion toward the inner layers is restrained by the density layering, so that the atoms diffuse in a quasi two-dimensional mode, which is comparable to our results presented in Figure 8. Therefore, the slow lateral and normal diffusion jointly contribute to the decrease of the attachment rate in the near-surface region, which would negatively affect the embryo growth. As evidence, the location of nuclei formed in the near-surface region is found to be close to the boundary between the near-surface and the

4 ΔG = − πr 3ρs Δμsl + 4πr 2σsl 3

(10)

where ρs is the number density of the bulk solid and σsl is the liquid−solid interfacial tension. The interfacial tension of the near-surface and the middle regions are approximately the same, because of only a slight difference in density between the two regions. The value of Δμsl can be roughly estimated by Δμsl = ΔHmΔT/Tm, where ΔHm is the latent heat of freezing and ΔT = Tm − T is the supercooling. From Figure 5a, the surface layering actually increases the average density of the nearsurface region, which is equivalent to an increase of pressure. For a 4% increase of density that corresponds to the maximum density of the surface layering, the equivalent compression of bulk water is up to 130 MPa. The increase of the local stress associated with surface layering will change the equilibrium melting temperature. Liquid water has a negative slope of the liquid-ice coexistence curve in the region p < ∼200 MPa. Under a pressure of 130 MPa, the melting temperature of bulk water is decreased by about 20 K. The latent heat of freezing at ambient pressure is 6.0 kJ/mol and decreases to 4.2−4.5 kJ/mol at 130 MPa. The interfacial tension ranges from 40 to 50 mJ/m2, which is referred to the published value of mW water at 220 K.24 With these parameters, we evaluated that the maximum ΔG (nucleation barrier) in the near-surface region is to be 3.4 times greater than that in the middle region, which is 10247



Article

with increasing nucleation barrier. Structural analysis shows that the near-surface region has a small number of I molecules. Thermodynamically, the surface layering characterized by the density oscillation normal to the surface compels the nearsurface region to high stress, which is shown to effectively increase the local nucleation barrier, lowering the nucleation rate. The nonhomogeneous distribution of stress in the nearsurface region hinders the mobility of the molecules near the surface, which also contributes to the decrease of nucleation rate. The present work suggests that decreasing the sample size is a feasible route to obtain high supercooling of water. It is of vital importance for investigating the complex phase transitions of supercooled water in “no man’s land”.

comparable to the results from the MFPT. This confirms that the surface layering is able to effectively raise the nucleation barrier, which is the main reason for the decrease of nucleation rate with film size. The nucleation barrier of water films increases as its thickness decreases, whereas the size of the critical nuclei remains almost the same. It cannot be interpreted within the framework of CNT. In fact, our analysis shows that the nucleation rate of the middle region is much larger than that of the near-surface region. Therefore, the nucleation process of films is dominated by the middle region, in which the critical nucleus preferentially forms. The total nucleation rate includes the surface and the volume terms. As we discussed above, its decrease merely results from the increased fraction of the nearsurface region that has a low nucleation rate. The partial nucleation rate of the middle region does not vary with the film size; thus, the critical nucleus size does not show dependence on the film thickness. From the discussion above, surface layering plays a crucial role in influencing nucleation rate. It was believed that surface density layering is beneficial for the formation of locally ordered clusters in water drops.49 In contrast, the present structure analysis shows that the local ordering near the surface decreases when the density oscillation arises. A plausible explanation is that the slopes of liquid−solid coexistence curves have opposite sign, and the locally ordered clusters are difficult to form in the dense surface layer. The mechanism of ice nucleation in the presence of the free surface displays a dependence of the interaction model of water in MD simulations. Using the sixsite interaction potential, the ice nuclei were observed to form preferably in the subsurface at 250 K, which was explained as a result of a large electric field induced by the orientational ordering of water molecules in the subsurface.50 This conclusion did not take the effect of the surface density layering into account. In fact, the density layering has been confirmed to be an intrinsic structure of the water surface, which is more pronounced for the sample with small size at low temperature.51 For nanosized water films or drops in a highly supercooled region, the surface density layering cannot be neglected in the discussion of nucleation mechanism. Recent experimental studies have achieved frozen water between 202 and 215 K that was well within “no man’s land”.52 The radii of droplets in the measurements ranged from 3.2 to 5.8 nm, which are close to the film thickness in our work. We speculate that the decreased nucleation rate may be the reason for the high supercooling.

■

ASSOCIATED CONTENT

S Supporting Information *

Details of the MFPT method, description of the CHILL algorithm, and additional results in the calculation of the nucleation rate. This material is available free of charge via the Internet at http://pubs.acs.org.

■

AUTHOR INFORMATION

Corresponding Authors

*E-mail: [email protected]. Tel:+86-10-68914027. *E-mail: [email protected]. Tel: +86-10-62782709. Author Contributions

Y.J.L. designed the study and wrote the paper. X.X.Z. performed the simulations, and M.C. guided the study and revised the manuscript. Y.J.L. and X.X.Z. contributed equally Notes

The authors declare no competing financial interest.

■

ACKNOWLEDGMENTS This work was supported by the National Natural Science Foundation of China under contract Nos. 51171027 and 50876049. We thank Shanghai Supercomputer Center for providing the computing resource.

■

REFERENCES

(1) Mishima, O.; Stanley, H. E. Decompression-Induced Melting of Ice IV and the Liquid−Liquid Transition in Water. Nature 1998, 392, 164−168. (2) Harrington, S.; Zhang, R.; Poole, P. H.; Sciortino, F.; Stanley, H. E. Liquid−Liquid Phase Transition: Evidence from Simulations. Phys. Rev. Lett. 1997, 78, 2409−2412. (3) Mishima, O.; Stanley, H. E. The Relationship between Liquid, Supercooled and Glassy Water. Nature 1998, 396, 329−335. (4) Debenedetti, P. G. Supercooled and Glassy Water. J. Phys.: Condens. Matter 2003, 15, R1669−R1726. (5) Swenson, J.; Teixeira, J. The Glass Transition and Relaxation Behavior of Bulk Water and a Possible Relation to Confined Water. J. Chem. Phys. 2010, 132, 014508. (6) Krämer, B.; Hübner, O.; Vortisch, H.; Wöste, L.; Leisner, T.; Schwell, M.; Rühl, E.; Baumgärtel, H. Homogeneous Nucleation Rates of Supercooled Water Measured in Single Levitated Microdroplets. J. Chem. Phys. 1999, 111, 6521−6527. (7) Hare, D. E.; Sorensen, C. M. The Density of Supercooled Water. II. Bulk Samples Cooled to the Homogeneous Nucleation Limit. J. Chem. Phys. 1987, 87, 4840−4845. (8) Sassen, K.; Liou, K. N.; Kinne, S.; Griffin, M. Highly Supercooled Cirrus Cloud Water: Confirmation and Climatic Implications. Science 1985, 227, 411−413.

4. CONCLUSION The aims of this study were to determine whether the nucleation rate of ice depends on the system size and the role of the surface on nucleation. We employed MD simulations to calculate the nucleation rate of ice in water films with different thicknesses. The results indicate that the nucleation rate decreases when the film thickness is decreased to the nanoscale. On the basis of CNT, we propose a theoretical model that well describes the nucleation process of ice in nanosized water films. The film has two different regions: the near-surface region, which has a low nucleation rate, and the middle region, which has a high nucleation rate. The nucleation thermodynamics of the two regions follow the description of CNT and do not significantly depend on the film size, whereas the volume fraction of the near-surface region tends to increase as the film size decreases. As a result, the total nucleation rate decreases 10248



Article

(9) Irie, H.; Kondo, Y. Evidence for the Nucleation of Polar Stratospheric Clouds inside Liquid Particles. Geophys. Res. Lett. 2003, 30, 1189−1192. (10) Volmer, M.; Weber, A. Keimbildung in Ubersattigten Gebilden. Z. Phys. Chem. 1926, 119, 277−301. (11) Oxtoby, D. W. Homogeneous Nucleation: Theory and Experiment. J. Phys.: Condens. Matter 1992, 4, 7627−7650. (12) Tabazadeh, A.; Djilaev, Y. S.; Reiss, H. Surface Crystallization of Supercooled Water in Clouds. Proc. Natl. Acad. Sci. U.S.A. 2002, 99, 15873−15878. (13) Kuhn, T.; Earle, M. E.; Khalizov, A. F.; Sloan, J. J. Size Dependence of Volume and Surface Nucleation Rates for Homogeneous Freezing of Supercooled Water Droplets. Atmos. Chem. Phys. Discuss. 2009, 9, 22929−22953. (14) Tabazadeh, A.; Djikaev, Y. S.; Hamill, P.; Reiss, H. Laboratory Evidence for Surface Nucleation of Solid Polar Stratospheric Cloud Particles. J. Phys. Chem. A 2002, 106, 10238−10246. (15) Duft, D.; Leisner, T. Laboratory Evidence for VolumeDominated Nucleation of Ice in Supercooled Water Microdroplets. Atmos. Chem. Phys. Discuss. 2004, 4, 3077−3088. (16) Earle, M. E.; Kuhn, T.; Khalizov, A. F.; Sloan, J. J. Volume Nucleation Rates for Homogeneous Freezing in Supercooled Water Microdroplets: Results from a Combined Experimental and Modelling Approach. Atmos. Chem. Phys. 2010, 10, 7945−7961. (17) ten Wolde, P. R.; Ruiz-Montero, M. J.; Frenkel, D. Numerical Calculation of the Rate of Crystal Nucleation in a Lennard-Jones System at Moderate Undercooling. J. Chem. Phys. 1996, 104, 9932− 9937. (18) Auer, S.; Frenkel, D. Prediction of Absolute Crystal-Nucleation Rate in Hard-Sphere Colloids. Nature 2001, 409, 1020−1023. (19) Auer, S.; Frenkel, D. Numerical Prediction of Absolute Crystallization Rates in Hard-Sphere Colloids. J. Chem. Phys. 2004, 120, 3015−3029. (20) Bolhuis, P. G.; Chandler, D.; Dellago, C.; Geissler, P. L. Transition Path Sampling: Throwing Ropes Over Rough Mountain Passes, in the Dark. Annu. Rev. Phys. Chem. 2002, 53, 291−318. (21) Wedekind, J.; Strey, R.; Reguera, D. New Method to Analyze Simulations of Activated Processes. J. Chem. Phys. 2007, 126, 134103. (22) Wedekind, J.; Reguera, D. Kinetic Reconstruction of the FreeEnergy Landscape. J. Phys. Chem. B 2008, 112, 11060−11063. (23) Matsumoto, M.; Saito, S.; Ohmine, I. Molecular Dynamics Simulation of the Ice Nucleation and Growth Process Leading to Water Freezing. Nature 2002, 416, 409−413. (24) Molinero, V.; Moore, E. B. Water Modeled as an Intermediate Element between Carbon and Silicon. J. Phys. Chem. B 2009, 113, 4008−4016. (25) Moore, E. B.; Molinero, V. Structural Transformation in Supercooled Water Controls the Crystallization Rate of Ice. Nature 2011, 479, 506−509. (26) Nosé, S. A Unified Formulation of the Constant Temperature Molecular Dynamics Methods. J. Chem. Phys. 1984, 81, 511−519. (27) Hoover, W. G. Canonical Dynamics: Equilibrium Phase-Space Distributions. Phys. Rev. A 1985, 31, 1695−1697. (28) Moore, E. B.; de la Llave, E.; Welke, E.; Scherlis, D. A.; Molinero, V. Freezing, Melting and Structure of Ice in a Hydrophilic Nanopore. Phys. Chem. Chem. Phys. 2010, 12, 4124−4134. (29) Lundrigan, S. E. M.; Saika-Voivod, I. Test of Classical Nucleation Theory and Mean First-Passage Time Formalism on Crystallization in the Lennard-Jones Liquid. J. Chem. Phys. 2009, 131, 104503. (30) Li, T.; Donadio, D.; Russo, G.; Galli, G. Homogeneous Ice Nucleation from Supercooled Water. Phys. Chem. Chem. Phys. 2011, 13, 19807−19813. (31) Malkin, T. L.; Murray, B. J.; Brukhno, A. V.; Anwar, J.; Salzmann, C. G. Structure of Ice Crystallized from Supercooled Water. Proc. Natl. Acad. Sci. U.S.A. 2012, 109, 1041−1045. (32) Kuhs, F. W.; Sippel, C.; Falenty, A.; Hansen, T. C. Extent and Relevance of Stacking Disorder in “Ice Ic”. Proc. Natl. Acad. Sci. U.S.A. 2012, 109, 21259−21264.

(33) Yethiraj, A.; van Blaaderen, A. A Colloidal Model System with an Interaction Tunable from Hard Sphere to Soft and Dipolar. Nature 2003, 421, 513−517. (34) ten Wolde, P. R.; Ruiz-Montero, M. J.; Frenkel, D. Numerical Evidence for bcc Ordering at the Surface of a Critical fcc Nucleus. Phys. Rev. Lett. 1995, 75, 2714−2717. (35) Lü, Y. J.; Chen, M.; Yang, H.; Yu, D. Q. Nucleation of Ni−Fe Alloy near the Spinodal. Acta Mater. 2008, 56, 4022−4027. (36) Desgranges, C.; Delhommelle, J. Role of Liquid Polymorphism during the Crystallization of Silicon. J. Am. Chem. Soc. 2011, 133, 2872−2874. (37) Ostwald, W. Studien über die Bildung und Umwandlung fester Körper. Z. Phys. Chem. 1897, 22, 289−330. (38) Matsubara, H.; Pichierri, F. Mechanism of Diffusion Slowdown in Confined Liquids. Phys. Rev. Lett. 2012, 109, 197801. (39) Krynicki, K.; Green, C. D.; Sawyer, D. W. Pressure and Temperature Dependence of Self-diffusion in Water. Faraday Discuss. Chem. Soc. 1978, 66, 199−208. (40) Harris, K. R.; Woolf, L. A. Pressure and Temperature Dependence of the Self Diffusion Coefficient of Water and Oxygen18 Water. J.C.S. Faraday I 1980, 76, 377−385. (41) Todd, B. D.; Evans, D. J.; Daivis, P. J. Pressure Tensor for Inhomogeneous Fluids. Phys. Rev. E 1995, 52, 1627. (42) Ikeshoji, T.; Hafskjold, B.; Furuholt, H. Molecular-level Calculation Scheme for Pressure in Inhomogeneous Systems of Flat and Spherical Layers. Mol. Simul. 2003, 29, 101−109. (43) Magnussen, O. M.; Ocko, B. M.; Regan, M. J.; Penanen, K.; Pershan, P. S.; Deutsch, M. X-Ray Reflectivity Measurements of Surface Layering in Liquid Mercury. Phys. Rev. Lett. 1995, 74, 4444− 4447. (44) Regan, M. J.; Kawampto, E. H.; Lee, S.; Pershan, P. S.; Maskil, N.; Deutsch, M.; Magnussen, O. M.; Ocko, B. M.; Berman, L. E. Surface Layering in Liquid Gallium: An X-Ray Reflectivity Study. Phys. Rev. Lett. 1995, 75, 2498−2502. (45) Shpyrko, O. G.; Grigoriev, A. Y.; Steimer, C.; Pershan, P. S.; Lin, B.; Meron, M.; Graber, T.; Gerbhardt, J.; Ocko, B.; Deutsch, M. Anomalous Layering at the Liquid Sn Surface. Phys. Rev. B 2004, 70, 224206. (46) Shpyrko, O. G.; Streitel, R.; Balagurusamy, V. S. K.; Grigoriev, A. Y.; Deutsch, M.; Ocko, B. M.; Meron, M.; Lin, B.; Pershan, P. S. Surface Crystallization in a Liquid AuSi Alloy. Science 2006, 313, 77− 80. (47) Gonzalez, L. E.; Gonzalez, D. J. Ab initio Study of the Atomic Motion in Liquid Metal Surfaces: Comparison with Lennard-Jones Systems. J. Phys.: Condens. Matter 2006, 18, 11021−11030. (48) Li, T.; Donadio, D.; Galli, G. Ice Nucleation at the Nanoscale Probes No Man’s Land of Water. Nat. Commun. 2013, 4, 1887. (49) Zasetsky, A. Y.; Remorov, R.; Svishchev, I. M. Evidence of Enhanced Local Order and Clustering in Supercooled Water near Liquid−Vapor Interface: Molecular Dynamic Simulations. Chem. Phys. Lett. 2007, 435, 50−53. (50) Vrbka, L.; Jungwirth, P. Homogeneous Freezing of Water Starts in the Subsurface. J. Phys. Chem. B 2006, 110, 18126−18129. (51) Chacón, E.; Tarazona, P.; Alejandre, J. The Intrinsic Structure of the Water Surface. J. Chem. Phys. 2006, 125, 014709. (52) Manka, A.; Parshad, H.; Tanimura, S.; Wölk, J.; Strey, R.; Wyslouzil, B. E. Freezing Water in No-man’s Land. Phys. Chem. Chem. Phys. 2012, 14, 4505−4516.

10249


Size Effect on Nucleation Rate for Homogeneous Crystallization of

Recommend Documents