Heterogeneous Ice Nucleation Controlled by the Coupling of Surface

Jan 5, 2016 - Heterogeneous Nucleation on Concave Rough Surfaces: Thermodynamic Analysis and Implications for Nucleation Design. Dongming Yan ...
0 downloads 0 Views 5MB Size
Download PDF
Article pubs.acs.org/JPCC

Heterogeneous Ice Nucleation Controlled by the Coupling of Surface Crystallinity and Surface Hydrophilicity Yuanfei Bi, Raffaela Cabriolu,† and Tianshu Li*

J. Phys. Chem. C 2016.120:1507-1514. Downloaded from pubs.acs.org by KAOHSIUNG MEDICAL UNIV on 09/03/18. For personal use only.

Department of Civil and Environmental Engineering, George Washington University, Washington, D.C. 20052, United States ABSTRACT: The microscopic mechanisms controlling heterogeneous ice nucleation are complex and remain poorly understood. Although good ice nucleators are generally believed to match ice lattice and to bind water, counter examples are often identified. Here we show, by advanced molecular simulations, that the heterogeneous nucleation of ice on graphitic surface is controlled by the coupling of surface crystallinity and surface hydrophilicity. Molecular level analysis reveals that the crystalline graphitic lattice with an appropriate hydrophilicity may indeed template ice basal plane by forming a strained ice layer, thus significantly enhancing its ice nucleation efficiency. Remarkably, the templating effect is found to transit from within the first contact layer of water to the second as the hydrophilicity increases, yielding an oscillating distinction between the crystalline and amorphous graphitic surfaces in their ice nucleation efficiencies. Our study sheds new light on the long-standing question of what constitutes a good ice nucleator.

■

INTRODUCTION In understanding and controlling the formation of ice, perhaps the most crucial question to answer is what makes an effective ice nucleation center? The extensive studies on ice formation in supercooled cloud droplets may provide some useful insight into this question.1 In clouds, the leading candidates for heterogeneous ice nucleation centers have been identified to be bacteria, pollen grains, mineral dust (e.g., clay), emission from aircraft (e.g., soot), and high-molecular-weight organic compounds (long-chain alcohols).2 While general distinction exists in their chemical nature, these ice-nucleating agents share common structural features. For example, the surfaces of these substances may contain chemical groups capable of forming hydrogen bond with water. As a consequence, layers of water molecules can be hydrogen-bonded to the surfaces. The surfaces may also exhibit the ordering patterns that resemble the structure of ice. Therefore, water layers bound to surfaces may be ice-like, providing template for ice to nucleate. On the basis of laboratory experiments, general criteria for a surface to be a good ice nucleator (IN) were summarized,3 and these are (1) IN should be highly water-insoluble; (2) IN should have similar hydrogen bonds available at its surface; and (3) IN should have an arrangement of atoms or molecules on its surface as close as possible to that of water molecules in some low index plane of ice. However, it should be stressed that these criteria may only be suggestive but not predictive, and no unique correlation can be established between ice nucleation threshold and these surface characteristics. For example, lattice match can be a good indicator, but it cannot account for the fact that different materials with the same lattice mismatch with ice, e.g., long-chain alcohol4 and lead iodide,5 exhibit different freezing temperatures. Amorphous or poorly ordered materials such as soot6 are also known to be effective ice nucleator, © 2016 American Chemical Society

despite their noncrystalline nature. Crystalline soluble salts such as ammonium sulfate was also found to nucleate ice.7 In addition, recent laboratory experiment using droplet freezing technique8 further showed that nonclay mineral feldspar dominate ice nucleation in mix-phase clouds. Here we report direct computational evidence that the crystallinity and the hydrophilicity of the carbon surface are strongly coupled to dominate the kinetics of ice nucleation. By systematically varying the water−carbon interaction strength over a wide range and computing the ice nucleation rates explicitly, we observe a rich spectrum of heterogeneous ice nucleation behaviors on both crystalline and amorphous carbon surfaces. Remarkably, we find that only within certain ranges of hydrophilicity does the crystallinity of the carbon surface play an active role in nucleating ice, and within these ranges, crystalline carbon surface is found to be significantly more efficient than the amorphous. In other hydrophilicities, the role of crystalline ordering becomes negligible, and no appreciable difference is observed in the directly computed ice nucleation rates between the crystalline and amorphous carbon surfaces. Our study demonstrates that the commonly adopted individual criterion for good IN alone may not be sufficient for predicting the ice nucleation capacity of a surface, and points at ways of engineering surfaces for controlling ice formation.

■

METHODS Model of Amorphous Graphene. To study the role of crystallinity of the carbon surface, we amorphize the graphene

Received: October 5, 2015 Revised: December 23, 2015 Published: January 5, 2016 1507

DOI: 10.1021/acs.jpcc.5b09740 J. Phys. Chem. C 2016, 120, 1507−1514

Article

The Journal of Physical Chemistry C surface by introducing the Stone−Wales (SW) defects through the Wooten−Weaire−Winer (WWW) bond-switching Monte Carlo (MC) method.9 Although the introduction of SW defect to crystalline graphene causes out-of-plane buckling,10 our main purpose in this work is to understand the role of surface crystallinity, therefore the amorphous graphene model is kept flat for simplicity. In this approach, the carbon−carbon interaction is described by the Keating-like potential:11 V=

∑ i,j

k 2 [rij − r02]2 + 2

2 ⎡ r02 ⎤ ∑ h⎢ rki⃗ · rkj⃗ + ⎥ 2⎦ i ,j,k ⎣

(1)

where r0 = 1.42 Å which denotes the equilibrium CC bond length. The first and second terms of the potential describe the energies for bond stretching and bond bending, with k = 7.5 kcal/mol and h (h/k = 0.2) being the corresponding spring constants, respectively. The MC simulation includes two types of moves: a single particle displacement and a WWW bondswitching move. The latter move allows introducing five- and 7fold rings in the graphene sheet while preserving the local 3fold coordination of carbon atoms. The schematic diagram for the bond-switching move is shown in Figure 1. By progressively introducing the bond-switching move, one may disrupt the translational symmetry of the graphene while minimizing the coordinational defect. In our procedure of amorphizing graphene, each MC step consists of one attempt of single particle displacement and one bond-switching move, both performed randomly throughout the graphene sheet. The standard Metropolis MC procedure was employed. To allow the structure to relax sufficiently, particularly once the bondswitching move is accepted, the single particle move and the bond-switching move are performed at kBT = 0.33 and 39 eV, respectively. After 109 MC steps, the structure was quenched at kBT = 0.033 eV for obtaining realistic atomic structure for amorphous graphene at room temperature. Figure 1(b−d) shows the CC bond length distribution, CCC bond angle distribution, and the radial distribution function of carbon atoms in the amorphous graphene, respectively. The topology analysis shows the amorphous graphene network is composed of 35% pentagons, 37% hexagons, 22% heptagons, and 4.9% octagons, in agreement with the ring statistics (34.5%, 38%, 24%, and 4.5%, respectively) obtained by Kapko et al.12 Systems. Our molecular dynamics (MD) simulations were carried out using the monotonic-water (mW) model.13 The carbon−water interaction is described based on the two-body term of the mW model. To mimic a wide range of hydrophilicity of carbon surface, the water−carbon interaction strength ε is varied14−16 between ε0 and 10ε0, where ε0 = 0.13 kcal/mol is the original water−carbon interaction strength17 that reproduces the experimental contact angle (86°) of liquid water on graphite. The simulations include 4096 water molecules and 1008 carbon atoms, in a nearly cubic cell combined with a periodic boundary condition (PBC). The Nośe-Hoover thermostat was employed to simulate the isobaric−isothermal canonical ensemble (NPT), with a relaxation time of 1 and 15 ps for the temperature and pressure, respectively. A time step of 5 fs is used throughout the simulations. Calculation of Ice Nucleation Rate. The direct calculation of ice nucleation rates was conducted by employing the forward flux sampling (FFS) method.18 In FFS, the nucleation rate R was obtained through the product of initial flux rate Φ̇λ0 and the growth probability P(λB|λ0) (namely, R =

Figure 1. (a) Bond interchange using the WWW algorithm to introduce the Stone−Wales defect in graphene. The SW defect is generated by rotating a carbon−carbon bond by 90°. Within the WWW framework, this is achieved by first randomly selecting a bond, e.g., AC, followed by a random selection of two carbon atoms connecting to A and C, respectively, but on the opposite sides of bond AC, e.g., atom B and D. Then the bond AB and CD are replaced by new bonds AD and BC (dashed lines). The final amorphous graphene structure yields a distribution of the carbon−carbon bond length (b), the carbon−carbon−carbon bond angle (c), and the radial distribution function of carbon atoms (d).

Φ̇λ0P(λB|λ0)), both of which can be calculated directly by sampling the nucleation trajectories in the parameter space defined by the order parameter λ. The pB histogram analysis in our previous study19 has demonstrated that a good order parameter λ is the number of water molecules contained in the largest ice nucleus for both homogeneous and heterogeneous ice nucleation on graphitic surface. The ice-like water molecule is numerically identified by the local bond-order parameter q6 with a q6 > 0.5.20 The initial flux rate Φ̇λ0 is obtained by N0/t0V, where N0 is the number of successful crossings to the interface λ0 from basin A, i.e., the spontaneous formation of ice nucleus containing λ0 water molecules, t0 is the total time of initial sampling, and V is the volume. The growth probability P(λB|λ0) 1508

DOI: 10.1021/acs.jpcc.5b09740 J. Phys. Chem. C 2016, 120, 1507−1514

Article

The Journal of Physical Chemistry C

Table 1. Calculated Initial Flux Rate Φ̇λ0, Growth Probability P(λB|λ0), and the Final Nucleation Rate R for Heterogeneous Ice Nucleation on Crystalline and Amorphous Graphene, With Different Water−Carbon Interaction Strength ε at 230 K crystalline

amorphous

water−carbon strength

Φ̇λ0 m−3 s−1

1ε0 2ε0 3ε0 4ε0 5ε0 5.5ε0 6ε0 6ε0 (stretched) 7.5ε0 7.5ε0 (stretched) 9ε0 10ε0 1ε0 2ε0 3ε0 4ε0 5ε0 5.5ε0 6ε0 7.5ε0 9ε0 10ε0

1.11 × 10 9.71 × 1033 5.43 × 1033 5.78 × 1033 4.60× 1033 2.12 × 1033 4.01 × 1033 8.34 × 1033 7.22 × 1033 1.53 × 1034 5.21 × 1033 1.45× 1034 1.04 × 1034 2.17 × 1034 9.02× 1033 5.23 × 1033 1.95× 1033 3.60 × 1033 3.12 × 1033 6.00 × 1033 7.07 × 1033 3.58 × 1033 34

P (λB/λ0) 1.03 1.23 5.27 3.00 1.25 7.28 2.28 1.25 1.74 3.66 2.34 6.80 1.21 1.17 1.67 1.55 9.08 1.60 2.41 5.44 2.19 1.24

± ± ± ± ± ± ± ± ± ± ± ± ± ± ± ± ± ± ± ± ± ±

0.39 0.36 1.36 0.67 0.30 2.63 1.26 0.59 0.86 1.36 0.96 2.83 0.48 0.33 0.52 0.51 3.80 0.84 1.64 2.67 1.06 0.61

× × × × × × × × × × × × × × × × × × × × × ×

nucleation rate m −6

10 10−3 10−3 10−2 10−2 10−4 10−11 10−8 10−10 10−5 10−6 10−7 10−6 10−3 10−4 10−4 10−6 10−8 10−11 10−10 10−9 10−8

1.14 1.19 2.86 1.73 5.74 1.54 9.13 1.04 1.27 5.61 1.22 9.85 1.26 2.54 1.50 8.10 1.77 5.76 7.52 3.27 1.55 4.44

± ± ± ± ± ± ± ± ± ± ± ± ± ± ± ± ± ± ± ± ± ±

0.43 0.35 0.74 0.39 1.72 0.57 5.05 0.49 0.62 2.08 0.50 4.10 0.50 0.72 0.47 2.68 0.74 3.02 5.12 1.60 0.75 2.18

× × × × × × × × × × × × × × × × × × × × × ×

−3

s−1

28

10 1031 1031 1032 1031 1030 1022 1026 1024 1029 1028 1027 1028 1031 1030 1029 1028 1025 1022 1024 1025 1025

Figure 2. (a) Calculated heterogeneous ice nucleation rates on crystalline and amorphous graphene at 230 and 235 K. Inset shows the atomic structures of both crystalline and amorphous graphene. (b) Calculated density profiles of water along the direction normal to the water−graphene interface, for different water−carbon interaction strengths.

is computed through P(λB|λ0) = ∏i =n 1P(λi|λi−1), where P(λi|λi−1) is the crossing probability for which a trajectory starts from interface λi−1 and ends on interface λi. By firing a large number Mi−1 of trial shootings from the interface λi−1 and collecting Ni successful crossings at the interface λi, one estimates P(λi|λi−1) = Ni/Mi−1. The statistical uncertainty of P(λB|λ0) consists of both the variance of binomial distributions of Ni and the landscape variance of the configurations collected at the previous interface λi−1.21 More details of ice nucleation rate calculations by FFS can be found in refs 19 and 20. The computed heterogeneous ice nucleation rates for all the conditions are listed in Table 1.

graphene. At 230 K, the computed heterogeneous ice nucleation rate (∼1028 m−3s−1) is significantly higher than the homogeneous ice nucleation rates (∼1013 m−3s−1 for the mW model20 and ∼105 m−3s−1 for the TIP4P/Ice model22) obtained at the same thermodynamic condition, highlighting the importance of heterogeneous nucleation. Interestingly, at both temperatures investigated, no distinguishable difference is found on the computed ice nucleation rate between the two surfaces. The result may appear surprising in its first glimpse but can be understood by the fact that crystalline graphene neither binds water nor templates ice. In fact, the heterogeneous ice nucleation of graphene was recently attributed to the induced layering of water density perpendicular to water−carbon interface:14,16,17 Since the oscillation of water density near the water−carbon interface matches well the density profile of ice, and the motion of water molecules in the contact layers is well restricted within the plane, the space that these water molecules can explore, hence the entropic barrier of

■

RESULTS AND DISCUSSION We first carried out study of heterogeneous ice nucleation on carbon surface with the original water−carbon interaction strength (ε = ε0). Figure 2a shows the computed heterogeneous ice nucleation rates on both crystalline and amorphous 1509

DOI: 10.1021/acs.jpcc.5b09740 J. Phys. Chem. C 2016, 120, 1507−1514

Article

The Journal of Physical Chemistry C

carbon interaction strength, the crystalline graphene shows no appreciable difference from the amorphous in its ice nucleation ability, whereas it becomes much more efficient within other ranges. In particular, at ε = 5ε0 the crystalline graphene yields an ice nucleation rate almost 105 times higher than the amorphous, strongly suggesting that the crystallinity is a key factor for ice nucleation under such hydrophilicity. The results immediately suggest the existence of a coupling resulted from the two surface characteristics that dictates ice nucleation, and that the coupling strength varies nonmonotonically with the surface hydrophilicity. To understand the origin of this coupling behavior, we first examine the recently proposed layering mechanism. Figure 2b shows, however, that the density profile of water normal to the carbon−water surface is essentially insensitive to the change of surface crystallinity, within the entire range of hydrophilicity investigated in this work. Clearly, the observed difference of ice nucleation rates induced by crystallinity may not be explained by the layering mechanism. To shed light on its origin, we examine the structure of the critical ice nucleus forming on carbon surface. As heterogeneous ice nucleation aligns the basal plane of ice parallel to graphene, the ice structure in the contact layers are expected to play an active role in ice nucleation. Figure 4a,b shows the

ice nucleation, is effectively reduced. It has been further suggested16 that water molecules in the first layer experience a nearly uniform potential from graphene. In other words, water molecules in contact layers see carbon surface as an atom-less flat surface. Indeed, our calculated density of water (Figure 2b) shows that the absence of crystallinity in graphene does not alter the profiles of water density normal to the surface. This is consistent with the observed insensitivity of ice nucleation rate on surface crystallinity under the original surface hydrophilicity. The interesting results are obtained when water−carbon interaction strength ε is allowed to vary. Increasing ε makes carbon atoms bind water more strongly, and correspondingly the surface becomes more hydrophilic. As shown in Figure 3,

Figure 3. Variation of the calculated ice nucleation rates with water− carbon interaction strength ε (in the unit of the original strength ε0) at 230 K, for both crystalline and amorphous graphene. The blue diamonds indicate the calculated nucleation rate of ice forming on crystalline graphene with a stretched carbon−carbon bond length of 1.46 Å. Inset shows the ratio of the ice nucleation rates between the crystalline and amorphous graphene Rc/Ra, as a function of water− carbon strength. The red horizontal line indicates a ratio of one. Figure 4. Evolution of ice layer (blue) in the critical nucleus forming on graphene (cyan) with water−carbon interaction strength ε. Parts (a) and (b) show the ice layer forming in the first contact layer of water on the crystalline and the amorphous graphene, respectively, at the original strength ε0. With an increasing ε, the ice layer gains lattice registry with respect to the crystalline graphene, as in (c) and (d). Further hydrophilicity increase disfavors ice formation in the first contact layer, as in (e) where only dangling water molecules appear, and instead, facilitates ice nucleation in the second contact layer, as in (f) and (g). The ice layer gradually aligns registered with the crystalline graphene with the continuous increase of water−carbon strength, as in (h).

the variation of ε (hence hydrophilicity) yields a rich spectrum of heterogeneous ice nucleation behaviors on both crystalline and amorphous graphene. First, increasing surface hydrophilicity initially enhances the ice nucleation rates on both surfaces, which maximizes the ice nucleation efficiencies for the crystalline and amorphous at the water−carbon strength of 4ε0 and 2ε0, respectively. A further increase of hydrophilicity, nevertheless, starts weakening the ice nucleation capacities of both surfaces, until the minimum nucleation rates are reached at 6ε0 for both surfaces. Upon additional enhancement of hydrophilicity, both surfaces are found to regain their abilities of nucleating ice. Remarkably, in contrast to the behaviors observed at the original water−carbon strength ε0, the crystallinity of the graphene surface was found to play a significant role in facilitating ice nucleation when surface hydrophilicity is varied. As shown in Figure 3, there exist ranges of hydrophilicity where the crystalline graphene yields significantly higher ice nucleation rate than the amorphous graphene. This can be better illustrated by the inset of Figure 3 which shows the ratio of the ice nucleation rate resulted from the crystalline graphene to that of the amorphous surface, Rc/Ra, as a function of water− carbon interaction strength ε. It is evident that within the low (ε ≤ 2ε0) and midhigh (6ε0 ≤ ε ≤ 7.5 ε0) ranges of water−

structure of the first ice layer of a critical nucleus formed on the underlying substrate with the original ε = ε0, for both crystalline and amorphous graphene. As expected, the basal plane of ice does not appear to match the underlying crystalline graphene lattice. Correspondingly, the computed ice nucleation rates do not exhibit fundamental difference when the crystallinity of graphene changes. As the water−carbon strength increases to 3ε0, the first layer of ice and the graphene lattice are found to form a nearly commensurate structure, i.e., each six member ring of water encloses a six member ring of carbon, as shown in Figure 4c. On the basis of the underlying crystalline graphene lattice, the pseudo commensurate structure can be considered as a 3 × 3 supercell, containing three graphene unit cells along 1510

DOI: 10.1021/acs.jpcc.5b09740 J. Phys. Chem. C 2016, 120, 1507−1514

Article

The Journal of Physical Chemistry C

Figure 5. In-plane distribution of water molecules in the contact layers of graphene at 230 K. For clarity, the distribution is defined as −Ln[P(x,y)],16 where P(x,y) is the probability density of finding a water molecule located at (x,y) in the plane. Parts (a) and (c) show the contour of − Ln[P(x,y)] computed for the first water layers in contact with crystalline graphene and amorphous graphene, respectively, with a carbon−water strength of ε = 3ε0. As indicated by the contour legend, a darker color represents a higher probability of distribution. Part (b) is the zoom-in of a local region in (a), which shows water tends to be adsorbed above the center of carbon hexagonal ring. To demonstrate the influence of water−carbon interaction strength ε on the distribution of water in contact layers, (d)−(i) show the distribution −Ln[P(x,y)] as a function of both x and y for the same region as in (b), with different hydrophilicity. In general, the increasing hydrophilicity of crystalline graphene patterns the first layer of water, which resembles the underlying graphene lattice. When the water−carbon strength is high enough (as in (i)), the second layer of water is also patterned, as a result of the strong localization of water in the adsorption sites in the first layer and the local tetrahedral ordering of water.

the crystalline graphene now regains its efficiency for promoting ice nucleation. The molecular analysis of ice nucleus suggests that the higher ice nucleation efficiency of the crystalline graphene within 2ε0 ≤ ε ≤ 6 ε0 and ε ≥ 7.5ε0 is closely related to its ability of forming an ice layer nearly commensurate with the underlying graphene lattice. Essentially, this can be understood on the basis of the templating effect, albeit that the templating of ice occurs over the supercell of graphene lattice. The question now is why does this templating effect only occur at certain hydrophilicity? To answer this question, we examine the in-plane density of water in the contact layers of graphene. At a low water−carbon strength (ε0 ≤ ε ≤ 2ε0), i.e., when carbon atom binds water weakly, water molecules in the first contact layer only experiences a nearly uniform potential from the underlying graphene. This can be illustrated by the in-plane distribution of water density in Figure 5d. Within this range of hydrophilicity, the heterogeneous nucleation of ice can be well explained by the layering mechanism, and an increase of surface hydrophilicity simply yields a greater out-of-plane density oscillation (i.e., a greater layering length), thus enhancing ice nucleation

each direction. We note that at this hydrophilicity, the crystalline graphene already yields an ice nucleation rate about ten times higher than that of the amorphous graphene. As shown in Figure 4d, the pseudo commensurate structure becomes even more evident at 5ε0, but interestingly, disappears upon a further increase of water−carbon strength to 6ε0. Under this hydrophilicity (6ε0), it is found that the water molecules in the first contact ice layer no longer arrange themselves into an ice like structure (Figure 4e). Instead, the basal plane of ice forms in the second layer of water, and it is further observed that the second ice layer, similar to the f irst ice layer forming at the original water−carbon strength ε0, is incommensurate with the underlying graphene lattice (Figure 4f). We also note that at this water−carbon strength (6ε0), the computed ice nucleation rates for both crystalline and amorphous graphene become equivalent again, and reach their minima (see Figure 3). More interestingly, upon a further increase in hydrophilicity, e.g., ε = 9ε0, the basal plane of ice in the second layer appears to form the 3 × 3 pseudo commensurate structure, similar to the first ice layer formed on the crystalline graphene at 3ε0∼5ε0. Through the formation of the pseudo commensurate structure, 1511

DOI: 10.1021/acs.jpcc.5b09740 J. Phys. Chem. C 2016, 120, 1507−1514

Article

The Journal of Physical Chemistry C

the in-plane g(r) has undergone a continuous yet significant change. In particular, the second peak of g(r) shifts its position from about 4.5 Å to about 3.8 Å. Noticing that these two inplane radii correspond to the distances between the second nearest water molecules in a hexagon and a pentagon, respectively, the peak shift can be understood by the change of the structure of first water layer from the predominantly hexagonal-like patchwork at the low surface hydrophilicity (Figure 6(b)), to a predominantly smaller membered ring (pentagons) at the high hydrophilicity (Figure 6(c)). The prevalence of pentagons at high hydrophilicity is in fact very similar to that of water layers confined under high pressure,23 and can be attributed to water maximizing its interaction with the surface.24 In particular, Figure 6(a) shows that the hexagonpentagon structural change is prominent as the water−carbon strength increases from 5ε0 to 6ε0. This leads to a significant decrease of ice nucleation rate (Figure 3) because the first contact water layer becomes structurally incompatible to ice I. Second, since the underlying graphene structure is rigid with a fixed carbon−carbon bond length dCC = 1.4 Å, the formation of an ideal commensurate superstructure that matches both ice and graphene would require an in-plane lattice constant of ice of 4.2 Å, i.e., three times of dCC. This implies that the ice layer in such structure must be strained (compressed in this case), given that the equilibrium in-plane lattice constant of ice is 4.5 Å. Therefore, the strain energy cost of imposing such perfect match eventually makes the formation of the strained ice layer energetically unfavorable in the contact layer. To further demonstrate the role of strain, we increase the carbon−carbon bond length in crystalline graphene to 1.46 Å, in order to better match ice lattice. The calculated ice nucleation rate at both 6ε0 and 7.5ε0 indeed show that the elongated crystalline graphene yields significant enhancement on ice nucleation rate (Figure 3) relative to those of the unstrained graphene. As the first water layer becomes inactive, the nucleation of ice consequently starts occurring in the second layer of water when ε ≥ 6ε0. In this case, as the in-plane water density in the second layer is nearly uniform (Figure 5h), the crystallinity of the underlying graphene becomes inactive again in ice nucleation, as evidenced by the synchronization of the computed ice nucleation rates within 6ε0 ≤ ε ≤ 7.5ε0. Interestingly, a further increase of hydrophilicity (ε ≥ 9ε0) is found to nearly immobilize water molecules within the first water layer in their adsorption sites. The strong localization of water in the first water layer essentially turns it into an image of the underlying crystalline graphene sheet. The second layer of water, now under the influence of the patterning from the first layer and the water−water interaction, becomes also structured (Figure 5i) and tends to facilitate the formation of the commensurate ice-like structure. In other words, the first and the second layers of water at the high hydrophilicity act almost as the graphene substrate and the first layer of water at the low hydrophilicity, respectively. In this way, the crystalline graphene sheet gains renewed ice nucleation capability, making it again superior than amorphous graphene for nucleating ice. The findings in this work have a few important implications regarding our understanding of heterogeneous ice nucleation. Perhaps the most immediate implication is that neither surface crystallinity nor surface hydrophilicity alone may be a good indicator for the ice nucleation efficiency of a surface. Instead, the combined surface characteristics may yield a complex coupling that dominates ice nucleation behaviors. This may potentially explain why materials that have similar lattice

rates. The similar behaviors have also been observed in previous studies.14,16 When carbon binds water more strongly (ε ≥ 3ε0), the in-plane density distribution of water in the first layer begins affected by the underlying atomic structures and becomes nonuniform. Consequently, the out-of-plane density profile of water (layering mechanism) is no longer sufficient for understanding the nucleation of ice, and the in-plane density distribution must be taken into account. In the case of crystalline graphene, water molecules in the contact layer are found to preferentially locate themselves above the center of the carbon hexagonal ring (see Figure 5b), making those weak adsorption sites of water. Consequently, the density of water clearly displays a pattern reminiscent to graphene lattice, as shown in Figure 5a. The ordering of water density in the first layer, which partially matches the ice basal plane of ice, thus may effectively increase the possibility of forming an ice-like fragment. In the case of amorphous graphene, the disordered, less uniform water density in the first layer (Figure 5c) in fact frustrates the crystalline ordering required for ice nucleation, which leads to a decrease of ice nucleation rate. As a consequence, the difference in the ice nucleation efficiency between a crystalline and amorphous graphene starts increasing with hydrophilicity. As water−carbon interaction becomes even stronger (e.g., 6ε0), the water molecules in the first layer are further constrained around the adsorption sites (Figure 5f). This leads to two effects that both suppress ice nucleation on crystalline graphene surface. First, because the formation of the hexagonal ice patchworks on crystalline graphene requires that only two-thirds of the adsorption sites are filled while the rest one-third are empty (see Figure 4c and d), a full coverage of water with strong binding strength in fact structurally hinders ice nucleation.16 This can be further illustrated by the evolution of the in-plane pair correlation function g(r), where r is the inplane distance, for the first contact layer with the water−carbon strength ε, as shown in Figure 6(a). As hydrophilicity increases,

Figure 6. (a) Evolution of the calculated in-plane radial distribution function g(r) of water in the first contact layer with surface hydrophilicity. The water molecules within 4.5 Å from the graphene surface are defined as the first contact layer (see Figure 2(b)) and are counted in the in-plane g(r). The red arrow highlights the growth of the second peak at r ≈ 3.8 Å. The typical structures of water in the first contact layer of crystalline graphene are shown for (b) ε0, and (c) 6ε0. 1512

DOI: 10.1021/acs.jpcc.5b09740 J. Phys. Chem. C 2016, 120, 1507−1514

Article

The Journal of Physical Chemistry C

into ice by self-assembled crystalline monolayers of amphiphilic alcohols at the air-water interface. J. Am. Chem. Soc. 1994, 116, 1179− 1191. (5) Bryant, G. W.; Hallett, J.; Mason, B. J. The epitaxial growth of ice on single-crystalline substrates. J. Phys. Chem. Solids 1960, 12, 189− 195. (6) Popovicheva, O.; Kireeva, E.; Persiantseva, N.; Khokhlova, T.; Shonija, N.; Tishkova, V.; Demirdjian, B. Effect of soot on immersion freezing of water and possible atmospheric implications. Atmos. Res. 2008, 90, 326−337. (7) Zuberi, B.; Bertram, A. K.; Koop, T.; Molina, L. T.; Molina, M. J. Heterogeneous Freezing of Aqueous Particles Induced by Crystallized (NH 4) 2SO 4, Ice, and Letovicite. J. Phys. Chem. A 2001, 105, 6458− 6464. (8) Atkinson, J. D.; Murray, B. J.; Woodhouse, M. T.; Whale, T. F.; Baustian, K. J.; Carslaw, K. S.; Dobbie, S.; O’Sullivan, D.; Malkin, T. L. The importance of feldspar for ice nucleation by mineral dust in mixed-phase clouds. Nature 2013, 498, 355−358. (9) Wooten, F.; Winer, K.; Weaire, D. Computer Generation of Structural Models of Amorphous Si and Ge. Phys. Rev. Lett. 1985, 54, 1392−1395. (10) Ma, J.; Alfè, D.; Michaelides, A.; Wang, E. Stone-Wales defects in graphene and other planar sp2-bonded materials. Phys. Rev. B: Condens. Matter Mater. Phys. 2009, 80, 033407−4. (11) Keating, P. N. Effect of invariance requirements on the elastic strain energy of crystals with application to the diamond structure. Phys. Rev. 1966, 145, 637. (12) Kapko, V.; Drabold, D. A.; Thorpe, M. F. Electronic structure of a realistic model of amorphous graphene. Phys. Status Solidi B 2010, 247, 1197−1200. (13) Molinero, V.; Moore, E. B. Water Modeled As an Intermediate Element between Carbon and Silicon. J. Phys. Chem. B 2009, 113, 4008−4016. (14) Lupi, L.; Molinero, V. Does hydrophilicity of carbon particles improve their ice nucleation ability? J. Phys. Chem. A 2014, 118, 7330− 7337. (15) Cox, S. J.; Kathmann, S. M.; Slater, B.; Michaelides, A. Molecular simulations of heterogeneous ice nucleation. I. Controlling ice nucleation through surface hydrophilicity. J. Chem. Phys. 2015, 142, 184704. (16) Cox, S. J.; Kathmann, S. M.; Slater, B.; Michaelides, A. Molecular simulations of heterogeneous ice nucleation. II. Peeling back the layers. J. Chem. Phys. 2015, 142, 184705. (17) Lupi, L.; Hudait, A.; Molinero, V. Heterogeneous Nucleation of Ice on Carbon Surfaces. J. Am. Chem. Soc. 2014, 136, 3156−3164. (18) Allen, R. J.; Frenkel, D.; Wolde, P. R. T. Simulating Rare Events in Equilibrium or Nonequilibrium Stochastic Systems. J. Chem. Phys. 2006, 124, 024102−024116. (19) Cabriolu, R.; Li, T. Ice nucleation on carbon surface supports the classical theory for heterogeneous nucleation. Phys. Rev. E 2015, 91, 052402. (20) Li, T.; Donadio, D.; Russo, G.; Galli, G. Homogeneous Ice Nucleation from Supercooled Water. Phys. Chem. Chem. Phys. 2011, 13, 19807−19813. (21) Allen, R. J.; Frenkel, D.; ten Wolde, P. R. Forward Flux Sampling-type Schemes for Simulating Rare Events: Efficiency Analysis. J. Chem. Phys. 2006, 124, 194111−194117. (22) Haji-Akbari, A.; Debenedetti, P. G. Direct calculation of ice homogeneous nucleation rate for a molecular model of water. Proc. Natl. Acad. Sci. U. S. A. 2015, 112, 10582−10588. (23) Johnston, J. C.; Kastelowitz, N.; Molinero, V. Liquid to quasicrystal transition in bilayer water. J. Chem. Phys. 2010, 133, 154516. (24) Cox, S. J.; Kathmann, S. M.; Purton, J. A.; Gillan, M. J.; Michaelides, A. Non-hexagonal ice at hexagonal surfaces: the role of lattice mismatch. Phys. Chem. Chem. Phys. 2012, 14, 7944−7949. (25) Ochshorn, E.; Cantrell, W. Towards understanding ice nucleation by long chain alcohols. J. Chem. Phys. 2006, 124, 054714.

mismatch with ice exhibit drastically different ice nucleation behaviors, and that no correlation has been established between ice nucleation threshold and any of the crystallographic characteristics.3 Therefore, a thorough understanding of the ice nucleation capacity for an IN should be achieved through a comprehensive study by explicitly considering all the necessary molecular details at the surface. Second, it is envisioned that the surface chemistry and surface crystallinity may also be coupled with the elasticity of the substrate. In our modeling the graphene surface is considered rigid. This is a good approximation as carbon−carbon bond is much stronger than the hydrogen bond of ice. When ice nucleates on a soft substrate that has a shear modulus lower than or comparable to that of ice, the possible local deformation of the substrate should play an active role if there exists a lattice mismatch between ice and the surface. Since a soft substrate may better template ice for a lower strain energy cost, the range of its lattice mismatch can be wider than a stiff substrate for achieving the comparable ice nucleation efficiency. This also may explain why some soft materials, e.g., self-assembled monolayers of amphiphilic alcohols,4,25 are known as the most effective INs.1 Last, the current work is carried out using a coarse-grained water model mW that does not explicitly account for hydrogen bond and electrostatic effect. In real water−surface interfaces, such effects are expected to be important for controlling initial ice crystallization. For example, the density profile of the hydrogen atom of water at the interface was found significantly different between the silver exposed and the iodide exposed AgI surfaces which have very different ice nucleation capabilities.26 However, the fact that heterogeneous ice nucleation already exhibits such a rich spectrum of behaviors even when simple models are employed14−17,19,27 simply suggests that the heterogeneous nucleation of ice, albeit ubiquitous and thermodynamically simple,19 is a complex process at the molecular level. Further studies are thus needed to elucidate the key microscopic factors controlling ice formation.

■

AUTHOR INFORMATION

Corresponding Author

*Phone: +1 (202) 994-3809. Fax: +1 (202) 994-0127. E-mail: [email protected] (T.L.). Present Address †

Laboratoire Interdisciplinaire de Physique (LiPhy), Université Joseph Fourier, Grenoble 1, BP 87, Saint Martin d’Hères, France.

Notes

The authors declare no competing financial interest.

■

ACKNOWLEDGMENTS The work is supported by NSF through award CMMI1537286. T.L. also thanks the Sloan Foundation through the Deep Carbon Observatory for supporting this work.

■

REFERENCES

(1) Murray, B. J.; O’Sullivan, D.; Atkinson, J. D.; Webb, M. E. Ice nucleation by particles immersed in supercooled cloud droplets. Chem. Soc. Rev. 2012, 41, 6519−6554. (2) Cantrell, W.; Heymsfield, A. Production of ice in tropospheric clouds: A review. Bull. Am. Meteorol. Soc. 2005, 86, 795−807. (3) Pruppacher, H.; Klett, J. Microphysics of Clouds and Precipitation; Kluwer Academic Publisher: Berlin, 2007. (4) Popovitz-Biro, R.; Wang, J. L.; Majewski, J.; Shavit, E.; Leiserowitz, L.; Lahav, M. Induced freezing of supercooled water 1513

DOI: 10.1021/acs.jpcc.5b09740 J. Phys. Chem. C 2016, 120, 1507−1514

Article

The Journal of Physical Chemistry C (26) Zielke, S. A.; Bertram, A. K.; Patey, G. N. A Molecular Mechanism of Ice Nucleation on Model AgI Surfaces. J. Phys. Chem. B 2015, 119, 9049−9055. (27) Fitzner, M.; Sosso, G. C.; Cox, S. J.; Michaelides, A. The Many Faces of Heterogeneous Ice Nucleation: Interplay Between Surface Morphology and Hydrophobicity. J. Am. Chem. Soc. 2015, 137, 13658−13669.

1514

DOI: 10.1021/acs.jpcc.5b09740 J. Phys. Chem. C 2016, 120, 1507−1514