Ice is born in low-mobility regions of supercooled liquid water.

When an ice crystal is born from liquid water, two key changes occur: (i) The molecules order and (ii) the mobility of the molecules drops as they adopt their lattice positions. Most research on ice nucleation (and crystallization in general) has focused on understanding the former with less attention paid to the latter. However, supercooled water exhibits fascinating and complex dynamical behavior, most notably dynamical heterogeneity (DH), a phenomenon where spatially separated domains of relatively mobile and immobile particles coexist. Strikingly, the microscopic connection between the DH of water and the nucleation of ice has yet to be unraveled directly at the molecular level. Here we tackle this issue via computer simulations which reveal that (i) ice nucleation occurs in low-mobility regions of the liquid, (ii) there is a dynamical incubation period in which the mobility of the molecules drops before any ice-like ordering, and (iii) ice-like clusters cause arrested dynamics in surrounding water molecules. With this we establish a clear connection between dynamics and nucleation. We anticipate that our findings will pave the way for the examination of the role of dynamical heterogeneities in heterogeneous and solution-based nucleation.

The freezing of water is one of the most ubiquitous phase transitions on Earth, shaping many processes from intracellular freezing (1) to cloud formation (2), yet there are major gaps in our molecular-level understanding of ice nucleation (3). Considerable effort has gone into understanding supercooled water and the microscopic details of ice nucleation (4–7). In particular, computer simulations have proved to be of great importance in understanding ice nucleation recently. They have for instance shown that the nucleation of hexagonal ice proceeds through stacking-disordered ice (8–11) and that this polytype is entropically stabilized at cluster sizes relevant to nucleation (12). While this has greatly advanced our understanding of the structural transformation during the nucleation process, nucleation is likely also influenced by the particle mobility, i.e., the dynamics of the liquid. The liquid dynamics during or shortly before ice nucleation have, however, not been directly studied or assessed in experiment. This is particularly striking since water exhibits dynamical heterogeneity (DH) (13), the phenomenon of coexisting spatially extended domains of very different mobility.

At the molecular scale, there are several studies that found structural differences between very mobile and very immobile liquid regions. Sciortino et al. (14) have shown that increasing local coordination leads to higher diffusion of molecules. Molinero and coworkers (8, 15) have shown that ice originates from four coordinated regions both in the pure liquid (8) and in salt-water solutions (15), as well as a growing correlation length of four coordinated water patches (16). Regarding bond-orientational order, Tanaka and coworkers (17) have shown that less mobile regions have a higher degree of tetrahedrality, consistent with other studies (18, 19), and they also identified that 5-membered rings could act as locally favored structures (20). Recently they have also shown that the emergence of DH can be described by a two-state picture of ordered and disordered regions (21).

The preordering in less mobile domains could be seen as an indication that nucleation will be preferred in these regions (22). On the other hand, Mazza et al. (23) found that rotational and translational heterogeneities correlate in water, and thus, the transformative motions of the rearrangements necessary for crystallization could be hindered in strongly immobile regions. Indeed, it has been argued for metallic liquids (24) that the enhanced mobility near the surface of nanowires can explain surface-induced crystallization rather than the traditional view of a heterogeneous reduction of the nucleation barrier. Mobile molecules would have the potential for a more effective exploration of phase space and therefore the ability to more readily perform collective self-assembly. Overall, one can find arguments for both immobile and mobile regions being preferential domains for nucleation. But whether there is any preference at all is not established to date for water or indeed any other supercooled liquid.

Here, we fill this gap by performing computer simulations of supercooled liquid water and ice nucleation which show that (i) ice nucleation occurs in low-mobility regions of the liquid, (ii) there is a dynamical incubation period in which the mobility of the molecules drops before ice-like ordering, and (iii) ice-like clusters cause arrested dynamics in surrounding water molecules.

Uncovering DH

We begin by characterizing the extent of DH in supercooled liquid water represented by the atomistic Transferable Intermolecular Potential, 4-point, Ice (TIP4P/Ice) (25) model. To this end we perform molecular dynamics (MD) simulations in a homogeneous water system in the temperature range 230–273 K, using iso-configurational analysis (ISOCA) (26, 27). This technique allows us to obtain spatially resolved maps of DH. We quantify the tendency of each molecule to move with a dynamical propensity (DP),where is the position vector of molecule at time , is the time of maximum heterogeneity (definition in ), and is the mean-square displacement of all oxygen atoms. In this approach we average the outcome over many trajectories that start from the same initial configuration but with a different set of random velocities, indicated by the notation . By doing this we are able to evaluate the effect of structure alone on DH. In Fig. 1 we show two snapshots of initial configurations used for the ISOCA at 230 K and 273 K, where each oxygen atom is colored according to its . This choice of temperatures is to illustrate the maximum difference in the extent of DH as well as the relevance of strong supercoolings like 230 K in homogeneous nucleation of ice (11, 12, 28). It can be seen that spatially localized domains of relatively immobile (blue) and mobile (red) particles emerge and that the spatial aggregation of the domains differs drastically. As reported in Fig. 1, the probability density distribution of the at 230 K is rather broad with a factor of 30 between the mobilities of particles at opposite tails. From the distribution of the DP we select the top and bottom 5% and label the corresponding molecules as most mobile (MM) and most immobile (MI) regions for further analysis. We show in that the choice of this threshold has no major influence on our results.

Fig. 1.

Dynamical heterogeneity in supercooled liquid water with the TIP4P/Ice model. (A) Spatial distribution of the dynamical propensity (DP) at 230 K and 273 K. Molecules (only oxygens shown) are colored according to the scale at Left. (B) Probability density distribution of the DP at 230 K. Blue and red shaded regions highlight the 5% of water molecules labeled as most immobile (MI) and most mobile (MM).

Dynamics of Precritical Fluctuations

The simulation of nucleation with atomistic water models currently remains a challenge, coming at enormous computational cost (11). Hence, as a first step to understand the connection between DH and ice nucleation, we focus on precritical clusters, i.e., the ice-like clusters that form via frequent thermal fluctuations and thus are readily probed by unbiased MD (29). The key results of our simulations with the TIP4P/Ice model are reported in Fig. 2, where we find a strong tendency for the precritical ice nuclei to form within MI domains, rather than within MM domains. To quantify their preference to form in the immobile regions we have split the whole range of (sorted) DP values into 20 equal sections (i.e., each corresponding to 5% of the whole range). This means that the molecules in the first (last) DP section are the same molecules as the ones in the MI (MM) region. In Fig. 2 we plot the average overlap of molecules in the biggest ice-like cluster with these DP sections. The expected overlap in the absence of any correlation with the DP would be 5%. We clearly see that there is a strong preference for precritical ice clusters to belong to the lower DP sections (i.e., more immobile molecules). In addition, we find that the formation of precritical clusters is suppressed in the MM domains as the overlap values there are below the baseline of 5%. Considering the fact that nucleation is stochastic in nature, this is strong evidence for the connection between immobile and ice-forming regions.

Fig. 2.

Connection between DH and precritical cluster formation in the TIP4P/Ice model. (A) Representative snapshot of a spontaneously formed cluster (red bonds and spheres) immersed in the MI region (transparent blue surface representation). (B) Average overlap between the molecules in the largest ice-like cluster and molecules in the relative DP range. Each bar corresponds to a 5% fraction of (sorted) DP values; i.e., the first (last) bar corresponds to the MI (MM) region. The expected overlap if clusters were uncorrelated with the DP would be 5% (indicated by the dashed line). (C) Average evolution of the mobility (DP) and crystallinity () for molecules in a cluster before its first time of assembly (taken to be ). Dashed lines indicate the mean values of DP and of the liquid. is the structural relaxation time. Shaded regions indicate 95% confidence intervals. All data were obtained with the TIP4P/Ice water model at 240 K.

Having shown that precritical ice clusters strongly overlap with the immobile regions from their early stages we now study the temporal connection between immobility and clusters before their first time of assembly. In Fig. 2 we show the average value of the DP and the crystallinity parameter (9) of molecules that belong to a cluster at its first time of assembly (taken to be at ). It can be seen that the mobility drops at ps, which is significant compared with the structural relaxation time of the liquid: ps at 240 K (). In addition, this drop occurs earlier and is much less abrupt than the change in structure, which can be associated with the rise of the order parameter at about 400 ps before the assembly. This finding is crucial and it confirms that immobility on average precedes ice-cluster formation by a significant timespan. This can be thought of as a dynamical incubation period in which the dynamics of the molecules change before the structural change toward ice. While this mechanism is not necessarily orthogonal to the commonly applied reasoning of purely structural ordering, our results suggest that arguing in terms of a process that involves distinct dynamical and subsequent structural steps is a viable route for a better description of nucleation.

Connection Between Nucleation and Dynamics

The results reported above point strongly toward an interplay between structural motifs pertinent to nucleation and the dynamics of the system. Unbiased MD simulations, however, cannot directly sample nucleation events, except under extreme conditions close to the homogeneous nucleation temperature (8). In this section, we report results from transition path sampling (TPS) (30) simulations that allow us to sample many nucleation events at reasonably high temperature. Specifically, we study a system of water molecules represented by the coarse-grained monoatomic water (mW) model (31) at 235 K. In the case of the brute force simulations, we used an ISOCA to quantify mobility. However, in the case of TPS we harvest many more (7,500) reactive trajectories, making ISOCA for each frame of each trajectory impractical. To fully exploit the statistical sampling provided by TPS, we therefore use the enduring displacement (32, 33) formalism, which permits on-the-fly calculation of each particle’s mobility. Using this approach has the added benefit of allowing us to validate the robustness of our previous results (we show that there exists a correspondence between quantifying the mobility with enduring displacements and the DP method in ).

To identify regions of space as either ice-like or liquid-like and either immobile or mobile (and the boundaries separating them), we introduce the coarse-grained crystallinity field and immobility field . In brief, is obtained by smearing each ice-like molecule’s position with a positive, normalized Gaussian and each liquid-like molecule with a negative Gaussian. In regions that are predominantly ice-like, will take values close to the density of ice, whereas in regions that are predominantly liquid-like, it will take values close to the negative density of liquid water. Boundaries separating ice-like and liquid-like regions are defined by surfaces with . In , we show that performs well at identifying regions as ice-like that are consistent with our intuitive understanding, while simultaneously neglecting small fluctuations. In a similar fashion, is obtained by smearing immobile particles with a positive Gaussian and mobile particles with a negative Gaussian.

In Fig. 3 we show a snapshot of and before nucleation [only regions and are shown] from a typical trajectory harvested from TPS. While there are no regions identified as ice-like, we do see large immobile domains in the supercooled liquid. Fig. 3 shows similar snapshots after the onset of nucleation. Crucially, and consistent with our findings for the precritical nuclei, it is clear that the ice nucleus forms within an immobile domain. In addition to this trajectory, we calculated statistical quantities characterizing the whole ensemble of TPS trajectories. These are discussed in and show that the behavior identified in Fig. 3 is indeed typical for all nucleation trajectories.

Fig. 3.

Ice nucleation occurs in relatively immobile domains of supercooled water. Shown is time evolution of the coarse-grained immobility (translucent silver) and crystallinity (opaque blue) fields, from a trajectory harvested by TPS with the mW model. (A) Before nucleation we see large immobile domains and an absence of crystalline order. (B and C) During nucleation, the ice nucleus forms (B) and grows (C) within the immobile domain. The ice cluster in snapshots B and C comprises 83 and 296 molecules, respectively. The diameter of the ice-like region in C is ∼3.4 nm.

Overall, the results obtained from TPS support and extend our observations from the unbiased simulations and strengthen our conclusion since we now also have insight (i) for larger clusters, (ii) for a large ensemble of nucleation trajectories, and (iii) with a different method of classifying immobility. In Fig. 3 it can also be seen that the surrounding immobile region is larger than the ice-like region, which suggests that ice-clusters further slow down their surroundings. Indeed, we also find for the atomistic model () that water molecules within ∼2.5 hydration shells have arrested dynamics as their DP values are substantially lower than in the bulk. This is relevant to crystal growth and theoretical modeling as the liquid molecules in the direct vicinity are significantly less mobile than those in the bulk.

Structural Hallmarks of Nucleating Regions

Given the presented evidence that shows that ice nucleation occurs in relatively immobile regions rather than in mobile ones, it is interesting to investigate the structural differences between these two domains. As mentioned in the Introduction, there have been a number of works already highlighting different kinds of preordering, i.e., regarding tetrahedrality (14, 19) or coordination number (8, 14). We add to this by analyzing the distribution of primitive rings (i.e., not divisible into smaller ones) in the regions of extreme mobilities (MM and MI as defined in Fig. 1) for the atomistic water model. As can be seen from Fig. 4 the MI regions have a rings distribution strongly peaked around members while the distribution for the MM domains is very broad. Moreover, the amount of entirely H-bonded rings is substantially higher in the MI region. Indeed, if we were to regard hydrogen bonds between members as a necessary criterion for being a ring [as is done in other literature (11, 34)], the MM region would be almost free of rings. In particular, an abundance of 6 -membered hydrogen-bonded rings can be regarded as the key structural characteristic of the MI domains in the liquid.

Fig. 4.

Structural differences in regions of adverse mobility in TIP4/Ice water. Shown is the number of -membered primitive rings within the respective domain at 230 K. The dashed portions of the bars represent the fraction of those rings fully connected by H bonds. Top Insets show an example of a fully and nonfully H-bonded 5-membered ring, where solid lines between oxygens are a guide to the eye and do not imply H bonds.

These results add to the findings of Haji-Akbari and Debenedetti (11), who showed that the nucleating ice nucleus exhibits a similar rings distribution, and Pirzadeh et al. (34), who highlight the presence of -membered rings near growing ice surfaces. We stress here, however, that in the liquid snapshots we investigated for the rings analysis, we find only negligible amounts of actual ice (according to different criteria; for details see ). Since the majority of 6-membered rings in the MI domains can be seen as ice-like if regarded in isolation, this means that it is the relative orientation between rings that is different from the crystal and thus the missing ingredient in forming ice. Because this happens in the MI region, which has a reduced diffusivity compared with other regions, we can speculate that the mechanism giving rise to the initial formation of ice-like clusters in the liquid is collective in nature [a similar argument based on density changes was made by Errington et al. (18)]. This would be consistent with a picture of reorienting rings rather than a picture of single-particle attachments via diffusive motion.

Discussion and Conclusions

Our results have established a clear link between DH and ice nucleation, suggesting that the complex nature of dynamics should not be overlooked in theoretical descriptions of nucleation. We have shown that liquid molecules in the vicinity of a nucleus are slowed down significantly, which implies that their diffusivity (connected to the attachment rate) is reduced compared with bulk molecules, an aspect neglected by classical nucleation theory.

We verified that our results on the rings distribution also hold for the coarse-grained mW model (31) even though for mW the extent of DH is much smaller (). This means that the characteristic features identified in our study are not sensitive to the specific hydrogen bond parameterization, but rather caused by the tetrahedral order inherent in the modeled material. Thus, our findings may be of relevance not just to water but to a much broader range of materials, evidently ones with tetrahedral order (such as group IV elements and silica). More broadly, it remains to be seen whether the immobility of precrystalline structures is connected to nucleation in nontetrahedral liquids in the same manner since the population of, e.g., rings is material specific. However, based on our findings we can suggest that it is the correspondence between immobile and crystalline topological features (such as rings) that connects immobile regions with nucleation.

The connection of DH and nucleation in water could be probed experimentally by investigating heavier water molecules as their mobility might be different. It is for instance established that liquid O has a higher nucleation rate than O (35). However, this cannot be taken as direct evidence in support of our observation as the change in the hydrogen bonding induced by nuclear quantum effects (36) potentially influences the nucleation rate too. A more rigorous experimental validation of our findings would be the nucleation rate comparison for H218O and ordinary water, which to the best of our knowledge has not been achieved. If one of the two liquids is more (less) mobile (diffusive), our results suggest a decreased (enhanced) nucleation rate.

Our findings are likely of relevance to heterogeneous nucleation and nucleation from solution. Generally, previous work has focused on investigating the structural or templating role of nucleating agents (37–40). However, an impurity or substrate is bound to impact the dynamics of the supercooled liquid in its vicinity, possibly leading to a novel mechanism of heterogeneous nucleation. For the example of water freezing, ice nucleation on hydrophobic surfaces (basically incapable of structuring the water network to a major extent) has been reported (37, 38) as well as intriguing alternating hydrophilic–hydrophobic patterns in the prominent ice-nucleating bacteria Pseudomonas syringae (41) and a nucleation enhancement by soluble molecules (42) that could be connected to the liquid dynamics. Moreover, for nucleation from solution it is well known that different solutes change the nucleation mechanisms, i.e., in the case of urea (43). Understanding how solutes change the dynamics and impact the formation of amorphous precursors could shed light on this issue. As such, we hope that this work will push the community to take into account the role of dynamics and particularly of DH in connection with crystal nucleation and growth.

In conclusion, we have shown that ice nuclei originate within immobile regions of the supercooled liquid and that there is a dynamical incubation period in which the mobility of particles drops before any structural change. Additionally, the presence of an ice crystallite causes arrested dynamics in water molecules that surround it and the distribution of rings can be seen as the structural hallmarks of DH in water. This connection between dynamics and structure provides another perspective on the physics of nucleation.

Supporting Information ().

We provide additional material in about (i) the role of statistics in calculating the DP, (ii) the DP threshold choice for the domains, (iii) the correspondence between DP and inherent structure displacements, (iv) the domain characterization with different order parameters, (v) DH in the coarse-grained mW model of water, (vi) the TPS simulations, (vii) arrested dynamics around ice clusters, and (vii) a visual impression of the connection between nucleation and DH which is provided in Movie S1.

Materials and Methods

Molecular Dynamics Simulations.

We mainly study the DH of a system containing 10,000 water molecules, represented by the TIP4P/Ice model (25). All our MD simulations are performed with the large-scale atomic/molecular massively parallel simulator (LAMMPS) code (44), integrating the equations of motion with a 2-fs time step and using a 10-fold Nosé-Hoover chain (45) with a relaxation time of 200 fs to control temperature. We use a cubic simulation box with 3D periodic boundary conditions and approximate volume of 68 68 Å3. Static bonds and angles have been constrained with the SHAKE algorithm (46). To avoid quenching effects upon generating starting configurations at different temperatures we performed (after 10 ns equilibration at melting temperature) a 0.5-K/ns cooling ramp in the NPT ensemble (10-fold Nosé–Hoover chain barostat with relaxation time of 2 ps). At 273 K, 260 K, 250 K, 240 K, 230 K, 220 K, and 210 K we save configurations. Those are propagated for 10 ns in the constant number of particles, volume, and temperature (NVT) ensemble at equilibrium volume to calculate dynamical properties as well as draw five snapshots for each temperature that are apart at least 1 ns to use for the ISOCA.

Work Flow to Characterize the Liquid Dynamics.

To characterize the liquid dynamics appropriate length and time scales have to be chosen, which is achieved by the following procedure:The resulting values for and for all temperatures can be found in Table 1 together with the structural relaxation time that was obtained as the -relaxation value from .

Table 1.

Overview of length () and time () scales used to characterize DH and structural relaxation time at different temperatures with the TIP4P/Ice model

Quantity	273 K	260 K	250 K	240 K	230 K	220 K	210 K
q_0, Å−1	2.01	1.96	1.91	1.84	1.80	1.77	1.76
t_0, ps	5	11	27	115	620	—	—
τ_liq, ps	2	6	14	68	356	—	—

Because of the computational cost we did not consider 220 K and 210 K for the rest of the study.

From the NVT simulation of the system with 10,000 molecules at the target temperature we obtain the oxygen–oxygen radial distribution function where the sums consider oxygen atoms and their positions , is the liquid density, and the average is over all trajectory frames.

We calculate the isotropic structure factor and define as the value where has its first peak, with being a reciprocal length.

For this we calculate the quantity .

Via we obtain the self-intermediate scattering function and the dynamical susceptibility . Isotropic averages taken over 200 independent directions according to are evaluated.

The time of maximum heterogeneity is taken as the time where has its maximum, i.e., where the movements at the nearest-neighbor range are most heterogeneous.

Analyzing the Connection Between Water Structure and Dynamics.

Ice-like molecules were detected using an order parameter () according to Li et al. (9) as implemented in plugin for metadynamics 2 (PLUMED2) (47, 48). First we compute for each molecule the quantitywhere the sum goes over the neighbors of molecule , are spherical harmonics, and and are the relative orientational angles between the molecules and . We compute this quantity for all possible values of and store them in a vector with components. Finally, we calculate values according toFor the particular choice of we classify a molecule with a value of lq6 > 0.5 as ice-like and otherwise as liquid. Ice-like molecules are then grouped together if they are within 3.4 Å of each other, and we call the resulting entities ice-like clusters.

For all snapshots (five for each of the five different temperatures) we calculate the DP (Eq. ) by performing 200 independent MD runs in the NVT ensemble, starting with randomized velocities. As length of these production runs we choose the time . For our purpose it is sufficient to consider oxygens only when calculating the DP. The latter is then used to label the MM and MI molecules as the top (bottom) 5% of DP values in the respective snapshot. The structural properties of these regions are then established by calculating the number of primitive -membered rings within them [using the rigorous investigation of networks generated using simulations code (49)], as well as computing distributions of common order parameters such as or topological patterns such as cages. All results reported in the text are averages over the results for the five snapshots per temperature.

The results reported in Fig. 2 were calculated from a TIP4P/Ice trajectory at 240 K. We define as the time of first assembly of a cluster the frame that for the first time has the biggest ice-like cluster composed of molecules that have not been part of the biggest ice-like cluster in the 75 ps before that frame. This is slightly larger than the structural relaxation time of 68 ps at that temperature and larger or smaller choices did not qualitatively alter the findings.

TPS Simulations.

TPS (30) was performed using an in-house code interfacing with LAMMPS (44). The system comprised 4,000 mW molecules and a pressure of 1 atm was maintained using a barostat with a damping constant of 5 ps and a time step of 10 fs. Langevin dynamics were used to maintain a temperature of 235 K with a damping constant of 10 ps.

To define whether a trajectory was reactive or not, we used the size of the largest ice-like cluster as defined by Li et al. (9). This means we additionally include each ice-like molecule’s nearest neighbor into the cluster to effectively include a surface contribution. The system was considered to be liquid if 50 and ice if 800. An initial reactive trajectory was generated by unbiased simulation at 205 K and the velocities were rescaled by a factor . An equilibration of 1,000 TPS moves was then performed with a 2:1 ratio of shooting to shifting moves. We made use of the one-way shooting algorithm (30). The maximum length of a shifting move was 80 ps. After this equilibration a further 7,500 TPS moves were performed as a production run.

To classify a molecule as either mobile or immobile we used the enduring displacement formalism (32)where Å, ps, and is the inherent structure position of molecule at time and the Heaviside step function . To find the inherent structure positions, the fast inertial relaxation engine (FIRE) algorithm was used (50). The coarse-grained immobility () and crystallinity () fields are then defined aswhere the sum runs over all molecules, and is a normalized Gaussian with Å. and if the molecule identified with the position vector is mobile and ice-like, respectively.

Code Availability.

Our custom TPS code is available from the corresponding author upon request.