Molecular simulations of heterogeneous ice nucleation. I. Controlling ice nucleation through surface hydrophilicity

Ice formation is one of the most common and important processes on earth and almost always occurs at the surface of a material. A basic understanding of how the physicochemical properties of a material's surface affect its ability to form ice has remained elusive. Here, we use molecular dynamics simulations to directly probe heterogeneous ice nucleation at a hexagonal surface of a nanoparticle of varying hydrophilicity. Surprisingly, we find that structurally identical surfaces can both inhibit and promote ice formation and analogous to a chemical catalyst, it is found that an optimal interaction between the surface and the water exists for promoting ice nucleation. We use our microscopic understanding of the mechanism to design a modified surface in silico with enhanced ice nucleating ability.

Ice formation is one of the most common and important processes on Earth and almost always occurs at the surface of a material. A basic understanding of how the physiochemical properties of a material's surface affects its ability to form ice has remained elusive. Here we use molecular dynamics simulations to directly probe heterogeneous ice nucleation at an hexagonal surface of a nanoparticle of varying hydrophilicity. Surprisingly, we find that structurally identical surfaces can both inhibit and promote ice formation and analogous to a chemical catalyst, it is found that an optimal interaction between the surface and the water exists for promoting ice nucleation. We use our microscopic understanding of the mechanism to design a modified surface in silico with enhanced ice nucleating ability.
Upon cooling, liquid water crystallizes into solid ice. Due to the presence of a free energy barrier separating the liquid and crystalline states, however, it is possible for liquid water to remain in a metastable 'supercooled' state to temperatures far below the equilibrium melting temperature. Heterogeneous ice nucleation, that is, ice nucleation in the presence of impurity particles such as mineral dust, soot or certain types of bacteria, generally increases the rate of ice nucleation and is the dominant process by which ice forms in nature. 1 Empirically, a large variance in the propensity of different materials to nucleate ice is observed, and due to the importance of ice formation in e.g. the climate sciences, much effort has been expended in identifying and cataloging the effectiveness of different materials to nucleate ice. 1 Despite the vast amount of research into heterogeneous ice nucleation, major gaps in our knowledge still exist, especially with regard to our understanding of the underlying chemical physics; this is reflected in our inability to accurately predict a material's ice nucleating efficiency and to answer seemingly simple questions such as how does hydrophilicity affect the ice nucleation rate? Not only is an understanding of the chemical physics of heterogeneous ice nucleation needed to predict the ice nucleating efficiency of existing materials, 2,3 but it is also paramount for the rational design of new materials to either promote or inhibit ice nucleation. Controlling ice formation is desirable in a variety of fields, for example, in the cryopreservation of cells, tissues and organs, 4 the food and transport industries and even as a potential means for climate control. 5,6 In contrast to fields such as chemical catalysis 7 and materials design, 8 there is currently no comprehensive set of design principles in terms of molecular 'descriptors' for making new substances to control ice formation. Put more simply, we do not know which are the a) Electronic mail: angelos.michaelides@ucl.ac.uk relevant microscopic properties of a material that determine its macroscopic ice nucleating efficiency. Often, the so-called 'requirements' for a good ice nucleating agent (INA) have been discussed, such as the requirement for a good crystallographic match to ice and the ability of water to chemically bond to the surface of the particle (i.e. hydrophilicity). 9 Although properties such as a good crystallographic match are important in heterogeneous nucleation of some systems, 10,11 such criteria have neither served as a full set of guidelines to identify good INAs, 1 nor have they aided the systematic improvement of ice nucleation inhibitors or promoters. Experimentally there is disagreement regarding the role of hydrophilicity. For example, Alizadeh et al. 12 have found ice nucleation to be slower on superhydrophobic surfaces, which they attribute not only to a lower contact area between the water and the surface, but also to a larger free energy barrier to nucleation. In contrast, Li et al. 13,14 found ice nucleation to be enhanced at hydrophobic modified silicon wafers relative to their unmodified hydrophilic counterparts, which was attributed to a faster dynamics of water at the hydrophobic interface. Also, it is found on kaolinite (a known hydrophilic INA) and platinum 15 that the most stable water overlayer can inhibit the growth of subsequent water layers. 16,17 Furthermore, in the case of requiring a good crystallographic match, evidence for icelike structures at surfaces is in general lacking 18 and there are also instances where materials with a good crystallographic match to ice are poor INAs. 19,20 We are therefore either faced with the prospect of relying on experiments to determine the efficacy of INAs on a case-by-case basis, or we can try and rationalize their behavior by elucidating the underlying molecular processes that control heterogeneous ice nucleation.
Here we present results from molecular dynamics (MD) simulations where we directly probe heterogeneous ice nucleation in the presence of a face centered cubic (FCC) nanoparticle (NP). The NP exposes its hexagonal (111) surface as its principal facet and can therefore act as a template for the hexagonal basal face of ice. With this NP completely immersed in water, shown in Fig. 1(a), we perform a series of studies in which we systematically explore the dependence of the nucleation rate on the hydrophilicity of the NP. By comparing these results to reference simulations of homogeneous nucleation we find a very interesting dependence of the nucleation rate on NP hydrophilicity; the NP can both promote and inhibit ice nucleation and exhibits a maximum nucleation rate at intermediate interaction strengths with the water. By examining the molecular level details of the nucleation processes at different hydrophilicities of the NP, we find that the structure in the immediate vicinity of the interface couples strongly with the nucleation rate. We then use this understanding of the underlying chemical physics to design an improved INA.
We have used the single site mW potential to model the interactions between water molecules, 21 which allows us to investigate length-and time-scales inaccessible to ice nucleation simulations that employ more traditional empirical potentials. 22 The NP was modeled as an FCC crystal with a lattice constant of 0.392 nm, consisting of 380 atoms. The FCC NP was hemispherical and exposed its (111) face as its primary facet. For the interaction between the NP atoms and the water molecules, the Lennard-Jones potential U (r) = 4 (σ/r) 12 − (σ/r) 6 was used, where r is the distance between an NP atom and a water molecule. The hydrophilicity of the NP was controlled by varying (a constant value of σ = 0.234 nm was used throughout). Interactions were truncated after 0.753 nm. This setup yielded contact layers at a height between 0.2-0.25 nm above the (111) surface of the NP, which is in reasonable agreement with values obtained from density functional theory calculations of water at metal surfaces. 23,24 All simulations were performed using the LAMMPS simulation package 25 with 2944 mW molecules in a periodic supercell. Previous simulation studies have suggested that the critical ice nucleus varies from ∼10 water molecules at 180 K 26 to ∼85-265 at 220 K 27,28 giving us confidence that our simulations should not be subject to serious finite size effects. Furthermore, simulations using a slab geometry (approximate dimensions of 66Å 2 ) with 4000 mW molecules confirm that the conclusions drawn from this work are not affected by changes to the box size and shape. For each value of the water-NP interaction energy, 16 MD simulations were performed at 205 K and 1 bar. Under these conditions, bulk liquid mW water is still metastable (as opposed to unstable) 29 but undergoes homogeneous nucleation on a timescale accessible to computer simulation such that statistically meaningful rates can be obtained. To detect 'ice-like' molecules, we have used the CHILL algorithm of Moore et al. 30 By monitoring the potential energy, we are able to determine the induction time to nucleation for each simulation and thus the probability P liq (t) that a given system remains liquid after a time t from the start of the simulation. We are able to determine the ice nucleation rate R by fitting P liq (t) = exp [−(Rt) γ ], where γ > 0 is also a fitting parameter. Further details of the fitting procedure and simulation setup are provided in the Supplemental Material. 31 In order to gauge the effectiveness of the NP as an INA, we have also studied bulk homogeneous nucleation, using identical settings.
Explicit simulations of heterogeneous ice nucleation (a) optimal (b) too weak (c) too strong have only recently started to emerge in the literature (see e.g. Refs. 22, 32-34) and to enable a systematic study, we draw conclusions from over 200 successful nucleation trajectories. Fig. 1(b) shows the dependence of the nucleation rate on the water surface interaction, the main finding of this study. Specifically, we have plotted log 10 (R/R hom ) vs E ads /∆H vap , where R hom is the bulk homogeneous rate and ∆H vap is the enthalpy of vaporization of bulk mW water (10.65 kcal/mol at 298 K). 21 A rich variety in the ice nucleating behavior is seen: the NP is seen to both promote and, surprisingly, inhibit ice formation. At low values of E ads , the heterogeneous nucleation rate is approximately two times lower than R hom . Thus when the particle is very hydrophobic, it tends to inhibit nucleation. As the water-surface interaction strength increases so too does the nucleation rate until it reaches a maximum at E ads /∆H vap ≈ 0.4 that is nearly 25 times faster than bulk homogeneous nucleation. Beyond the maximum, the rate steadily decreases until E ads /∆H vap ≈ 1.0. Further beyond this, the rate remains roughly constant and slightly below R hom .
We now try to understand this intriguing dependence of the nucleation rate on the hydrophilicity of the NP. To this end, we have examined in detail the mechanisms by which nucleation occurred on the NP for the various interaction strengths. As the (111) surface of an FCC crystal exhibits hexagonal symmetry, one possible mechanism for heterogeneous ice nucleation is a template effect whereby the molecules in the contact layer form an hexagonal structure commensurate with the surface. Fig. 2(a) confirms this, where we show a typical ice nucleation event in the presence of the NP with E ads /∆H vap ≈ 0.3 (close to the maximum rate). Here we can clearly see that the water molecules in contact with the (111) surface of the NP do indeed form an hexagonal structure commensurate with the surface that resembles the basal face of ice. We can also see that the water molecules directly above this contact layer also form a similar hexagonal structure. The surface is therefore acting to promote ice nucleation by providing an arrangement of adsorption sites that resemble the structure of ice, thereby stabilizing structural fluctuations towards ice-like arrangements in the liquid. Now that we have established that the NP acts to promote ice nucleation by acting as a template for ice, the dependence of the rate on E ads can easily be understood as a competition between water-water and water-surface interactions. In Fig. 2(b) we show the structure of water in contact with the NP for E ads /∆H vap ≈ 0.08 (the weakest E ads investigated, which inhibits ice nucleation). Clearly, such a weak water-surface interaction is unable to stabilize ice-like configurations and in fact, ice nucleation is seen to occur away from the surface in a homogeneous manner. Fig. 2(c), on the other hand, shows the structure of water in contact with the NP for E ads /∆H vap ≈ 1.2. While this NP also inhibits ice nucleation, it does so for the opposite reason: the water-surface interaction is too strong, meaning that the water molecules cannot rearrange to form an ice-like layer at the surface. It is also clear that the coverage is higher than when ice forms at the (111) surface, as shown in Fig. 2(a). For this strongly interacting scenario, we also see that ice forms away from the surface in a homogeneous manner. Although the coupling between the molecular mechanism and the ice nucleation rate as we change the surface hydrophilicity is actually rather simple, it has neither been demonstrated (color online) Surface modification to promote ice nucleation for strong water-NP interaction strengths (E ads /∆Hvap > 1.0). By introducing small adsorbates (colored yellow) at the (111) surface, the template effect can be recovered (c.f. Fig. 2(c)). The nucleation rate is increased by approximately a factor 50 compared to the bare NP for E ads /∆Hvap ≈ 1.2, as indicated in Fig. 1(b).
nor proposed previously. The advantage of such simplicity is that it can potentially lead to the control of ice nucleation, which we demonstrate next.
When an ice-like hexagonal overlayer forms at the (111) surface of the NP, such as in Fig. 2(a), it does so with sub-monolayer coverage i.e. not all of the available adsorption sites are occupied by water molecules. We refer to these unoccupied adsorption sites on the (111) terrace as "excess" sites. For E ads /∆H vap > ∼ 0.6, these excess sites are occupied for long times and for nucleation to occur, an area of decreased coverage at the surface must occur such that an hexagonal motif can form. This motif can then act as a template for the hexagonal basal face of ice (a movie showing this for E ads /∆H vap ≈ 0.9 is provided 31 ). When E ads > ∆H vap it becomes favorable for a water molecule to occupy a site on the surface, including the excess sites, rather than a position in the bulk liquid. This prevents the water molecules in the contact layer from forming the hexagonal arrangements required for ice nucleation at this NP. By this rationale, if the density of available adsorption sites was lower, then the template effect (and the enhanced nucleation rate) may be preserved at higher values of E ads . To this end, we have modified the (111) surface of the NP by adsorbing small molecules, at the excess sites, which only have a weak interaction with water 35 , and recomputed the nucleation rate with E ads /∆H vap ≈ 1.2. As seen in Fig. 3, the template effect is indeed recovered. It is also seen in Fig. 1(b) that this modified surface enhances ice nucleation by a factor of 50 compared to the unmodified surface. This is a clear demonstration of how the molecular insight into heterogeneous ice nucleation can be used to rationally design surfaces of different ice nucleating ability. Experimentally, this could be realized through adsorption of small molecules to the surface (e.g. carbon monoxide) or through surface alloying. In fact, alloying a platinum (111) surface with tin is observed to promote the formation of an hexagonal ice-like bilayer under ultrahigh vacuum conditions, 36,37 in a fashion analogous to our modified surface (the tin atoms at the surface act in part to reduce the density of adsorption sites). Although not reported, we also note that in our simulations, the nucleation rate can also be decreased by adsorbing small molecules such that the NP can no longer act as a template.
More generally, the sensitivity of the nucleation rate on surface hydrophilicity could be tested by e.g. using nanoparticles of gold or silica functionalized with organic molecules of varying hydrophobicity. In addition to using well established methods such as the droplet freezing techniques, 38 it may also be possible to exploit recent advances in femtosecond X-ray scattering techniques that have allowed real-time monitoring of homogeneous ice nucleation in micron sized water droplets. 39 Not only could such an experimental protocol be used to compare rates of ice nucleation in the presence of immersed NPs, but information regarding the impact of such NPs on the microscopic structure of the liquid should also be available.
In summary, we have used computer simulations to systematically compare heterogeneous ice nucleation rates in the presence of a nanoparticle of varying hydrophilicity. We have seen that the nanoparticle can promote ice nucleation by acting as a template for the hexagonal ice lattice, but that the ice nucleating efficiency is lost if adsorption is too strong, due to a high coverage of water molecules destroying the template effect. Modification of the surface such that the coverage of water molecules is reduced recovers this template effect and enhanced nucleation can be achieved for strongly adsorbing surfaces, clearly demonstrating how molecular level understanding of heterogeneous ice nucleation can be used to manipulate the rate of ice formation. The use of molecular descriptors to predict useful macroscopic properties of materials has been successfully used in other fields, such as chemical catalysis. 7 Designing new catalysts for reactions such as methanation (CO + 3 H 2 −→ CH 4 + H 2 O) has relied upon the establishment of a Sabatier principle based on a computationally tractable quantity (in this case the dissociation energy of CO at the surface). 40 We have seen that for the surface investigated in this study, the adsorption energy of a single water molecule can be used to describe the heterogeneous ice nucleation rate. Although a comprehensive set of rules still requires further experimental and theoretical investigation, the results presented here suggest that if the surface acts as a template for ice, then one must tune either the density of adsorption sites, or the propensity of water to adsorb to the surface. In Ref. 41, we show that the variation of the ice nucleation rate upon surface hydrophilicity is dependent upon the surface topography, demonstrating that combined effect of different surface properties needs to be considered when trying to understand what makes a good INA. Other properties such as the crystallographic match to ice and the role of surface defects are also likely to be important, and the results presented in this Communication serve as a platform upon which future studies can be conducted.
Dr. Gabriele Sosso is thanked for sharing his data from larger simulations in a slab geometry. We are grateful to the London Centre for Nanotechnology and UCL Research Computing for computation resources. S.M.K.