Signatures of sluggish dynamics and local structural ordering during ice nucleation

Ice nucleation is a stochastic, fluctuation driven process of paramount importance to areas ranging from biology and atmospheric science to aeronautics and space research and energy infrastructure engineering.^1–19 Nonetheless, the mechanisms and kinetics of homogeneous ice nucleation are challenging to characterize experimentally due to the limited spatiotemporal resolution of existing techniques.^20,21 Computer simulations, by contrast, can provide the molecular-level insights necessary to understand crystallization. However, the rare and stochastic nature of ice nucleation frustrates study of this process using standard computational methods such as molecular dynamics (MD).^20,21 Consequently, few computational studies have successfully investigated ice nucleation with standard MD.^22–24 

Recent efforts have been devoted to the development of advanced computational techniques to enhance sampling of rare fluctuations that lead to nucleation.^20,25–36 The nucleation trajectories generated using these methods, however, are typically either biased or noncontiguous, making characterization of the underlying crystallization mechanisms challenging. Interestingly, the ices observed in these computational studies are predominately cubic-rich stacking disordered crystals, consistent with those that nucleate homogeneously in the atmosphere and in vapor chamber experiments.^{4,20,22–24,29,37–42} This observation has been rationalized by invoking the Ostwald step rule^34,43,44 or geometric arguments based on arrangements of structural motifs that favor the cubic polymorph.²⁰ 

According to classical nucleation theory, the casual encounter of short-living ordered environments, each with random orientations, leads to the formation of a critical nucleus that is large enough to resist melting due to the balance between its surface and volume free-energy.⁴⁵ On the other hand, non-classical pathways have been reported for numerous systems,^46–55 including stacking-fault crystals,⁴² an atomistic model of water,⁴² phase-change materials,^56–58 proteins,^59,60 and amyloid fibrils,^61–67 among others. Hence, these studies collectively suggest that non-classical pathways might be more common than previously believed.

Here, we investigate the crystallization of the ST2 model of liquid water occurring under deeply supercooled conditions, below its liquid–liquid critical point, using unbiased MD simulations. Phase separation between the low-density and high-density liquids (LDL and HDL, respectively) occurs rapidly, as reported in Ref. 24. The LDL region is characterized by lower mobility^{22,35,68–70} and enhanced structural order. Within this region, small ordered pockets form and percolate to generate a sub-critical nucleus. In contrast with other crystallites that continuously form and disappear, the sub-critical nucleus resists melting and serves as a template for the critical nucleus. Remarkably, we observe that polymorphic selection favoring Ic over Ih occurs at the stage of formation of the sub-critical nucleus, similar to the case of hard spheres.⁷¹ Given the similar free energies of Ic and Ic,⁷² we posit that Ic is kinetically favored at this stage due to purely geometric considerations. Specifically, we suggest that the ordering of the molecules in water’s second neighbor shell into a local ice Ic-like motif is more likely due to ice Ic’s symmetry and is thus favored over the formation of Ih arrangements.

We observe that crystal growth following nucleation occurs mostly via the freezing of liquid water in contact with the critical nucleus surface and via the absorption of small Ic and Ih crystallites that remain trapped in the expanding crystal. The latter process results in the formation of a stacking fault disordered ice that is rich in the cubic polymorph, similar to those observed in previous computational and experimental studies.^{4,20,22–24,29,37–42,71} On the surface of the expanding nucleation seed, we observe the occasional formation of ice-0 sites,⁷³ which transform into Ih during growth of the crystal and hence contribute to the prevalence of stacking faults in the crystallized sample. We rationalize the presence of ice-0 by observing that the oxygen atoms of the water molecules in its second neighbor shell are arranged in an intermediate structure between the Ic and the Ih configurations. Thus, the overall picture that emerges from our study confirms and extends recent observations about the non-classical nature of crystal nucleation in atomistic water models and other systems^42,56 and suggests that polymorphic selection occurs in the liquid, as observed in hard spheres.⁷¹ 

This article is organized as follows. Section II describes the MD simulation protocol and the order parameter that is used to distinguish between molecules with local coordination environments similar to those found in liquid water, ice Ic, and ice Ih, respectively. In Sec. III A, we analyze the connection between local structural order and molecular mobility and discuss the role that they play in the formation of the nucleation seed. In Sec. III B, we characterize the nucleation process using clustering analysis and elucidate the molecular origins of the polymorphic selection favoring Ic that occurs during the formation of the nucleation seed and the subsequent formation of stacking faults in the crystallized sample. Finally, conclusions and final remarks are presented in Sec. IV.

We performed classical MD simulations of a system of N = 2500 water molecules described by the ST2 interaction potential.⁷⁴ A previous MD simulation study of the ST2 model demonstrated that it exhibits spontaneous crystallization on ∼μs time scales under deeply supercooled conditions²⁴ near its liquid–liquid critical point.⁷⁵ Whereas the liquid–liquid transition in the ST2 model of water has been studied extensively,⁷⁵ its crystallization behavior remains incompletely understood and is thus the focus of our study. Following Ref. 24, we generate crystallization trajectories by performing MD simulations in the NVT ensemble at a fixed mass density of ρ = 0.98 g/cm³, using a periodic rectangular simulation cell with relative x, y, z-dimensions of (L_x, L_y, L_z)/L_z = (1/4, 1/4, 1). This geometric configuration minimizes the interfacial free energy between LDL and HDL,^24,76,77 hence, favoring spontaneous phase separation below the liquid–liquid critical temperature T_LLCP ≈ 245 K.^78–80 The equations of motion were numerically integrated using the velocity Verlet algorithm⁸¹ with 2 fs timestep, and the temperature was maintained using a Nosé–Hoover chain thermostat⁸² with a 2 ps time constant. van der Waals and Coulombic interactions were truncated 0.78 nm, and long-range electrostatics were handled using the reaction field method.⁸³ The system was equilibrated at 300 K for 0.02 µs and then quenched to T = 235 K < T_LLCP by instantaneously adjusting the set point of the thermostat. In analyzing the behavior of the system, the point of the quench is taken to be the time origin (i.e., t = 0 µs).

We gain insight into nucleation by analyzing the MD trajectories using a recently developed local order metric (LOM)⁸⁴ that has been successfully applied to understand a variety of complex molecular processes in systems ranging from water^85–90 to phase change materials⁵⁶ and boron nitride nanotubes.⁹¹ We briefly describe the order parameter below; full details of its numerical implementation can be found in Ref. 84.

The coordination environment of an atomic site j in a snapshot of a molecular dynamics or Monte Carlo simulation is characterized by the local pattern formed by M neighboring sites. Typically, the M sites include the first and/or the second coordinating neighbors of the site j. In this study, j indexes the oxygen atoms on the water molecules in the system, and we consider the pattern formed by the M = 16 oxygen atoms on the neighboring molecules in the first and second coordination shells. The local reference structure is the set of the same M neighboring sites in an ideal lattice of choice, the spatial scale of which is fixed by setting its nearest neighbor distance equal to d, the average equilibrium value in the system of interest. For a given orientation of the reference structure and a given permutation $P$ of the pattern indices, we define the LOM S(j) as the maximum overlap between the pattern and reference structure in the j neighborhood as

where θ, ϕ, and ψ are Euler angles, $PiPj$ and $Rij$ are the pattern and the reference position vectors in the laboratory frame of the M neighbors of site j, respectively, and A^j is an arbitrary rotation matrix about the pattern centroid. The parameter σ controls the spread of the Gaussian functions (σ = d/4 in this work, where d is the characteristic length of the local pattern). The LOM satisfies the inequalities 0 ≲ S(j) ≤ 1. The two limits correspond, respectively, to a local pattern with randomly distributed points [S(j) → 0] and to an ordered local pattern matching perfectly the reference [S(j) → 1]. We also define a global order parameter based on S(j), as the average score function S,

The LOM and S are endowed with generality and a high-resolving power that allow the investigation of a variety of systems independent of their specific atomic/molecular nature.^56,84–91

We first examine the time evolution of the global order parameter S [Eq. (2)] and the percentage of ice-like sites f_ice during crystallization (Fig. 1). The latter quantity was evaluated by using the LOM S(j) to distinguish between liquid- and ice-like environments. Following Ref. 84, the local environment of j is defined to be liquid-like if S(j) < 0.80 and ice-like otherwise. The LOM was computed using three reference structures created using the position of the oxygen atoms in the second shell of Ic, Ih, and ice-0, respectively. For each water molecule j, the LOM is evaluated as the maximum value computed for the three reference structures S(j) = max{S_Ic(j), S_Ih(j), S_I0(j)}. The population of ice-0 structures is generally small and hence the analysis primarily reflects contributions from ice Ic- and Ih-like motifs. Nonetheless, as discussed in Sec. III B, ice-0 structures play an important role in the expansion of the nucleation seed.

Based on the behavior of S and f_ice, we identify three distinct windows, or regions, in time. In time window R₁ $(t=0–0.8μs)$⁠, the system is predominately liquid and ice-like crystallites spontaneously form and melt, constituting $≈10%$ of the sample at any given time. Interestingly, at t ≈ 0.62 µs (r₁ in Fig. 1), we observe slight increases in S and f_ice. These increases last for $≈0.15μs$ and result from attempts of the liquid phase to generate the nucleation seed. As discussed in Sec. III B, structural reorganization of small, transient crystallites begins at their surfaces and propagates inward to their cores, causing them to lose stability and melt. The attempts to create a stable nucleation seed become successful at t ≈ 0.80 µs (time window R₂, $t=0.8–1.0μs$⁠), and the formation of a stable seed is signified by the marked jumps in both S and f_ice. The nucleation seed consists of $≈300$ water molecules (vide infra), corresponding to roughly 12% of the sample. The plateaus in both S and f_ice in time window R₂ suggests that the nucleus remains stable for $≈0.1μs$⁠. During this period of time, the inner core of the seed is no longer affected by the perturbations caused by structural reorganizations occurring at its surface. At t ≈ 0.9 µs, the nucleation seed rapidly grows, as indicated by the increases in both S and f_ice, resulting in nearly complete crystallization at t ≈ 1.00 µs (time window R₃, t > 1.0 µs).

To elucidate any possible link between molecular mobility and structural ordering, we divide the simulation box into four cubic sub-boxes of equal volume (Fig. 2). In each sub-box, we compute the order parameter S and the number of switching events (or exchanges) in which a neighboring water molecule transitions from the first to the second coordination shell of a central molecule. The number of switching events per molecule is a measure of short-range molecular mobility, and we compute it by constructing a list of the four nearest neighbors of each molecule at time t and compare it with the list at earlier time t − 1. If a molecule was identified as a nearest neighbor at time t − 1 but not at t, then it switched from being a first to a second nearest neighbor. The local molecular mobility was characterized by computing n_sw, the number of switching events normalized by the average number of water molecules in each sub-box, using a lag time of 1 ps between t and t − 1.

We investigate structural ordering and molecular mobility by examining S and n_sw in the four sub-boxes of the simulation cell as a function of time (Fig. 2). In time window R₁, the order parameter S fluctuates around constant values in all sub-boxes [Fig. 2(a)], except in sub-box No. 1, where it begins to slightly increase at t ≈ 0.62 µs. This time is similar to where slight increases are seen in S and f_ice computed for the entire system (r₁ in Fig. 1), indicating that the attempts to generate the nucleation seed are localized inside the sub-box No. 1. The number of switches n_sw in sub-box No. 1 is also remarkably lower than in the other sub-boxes, suggesting that molecules in this region exhibit lower mobility. This observation suggests a direct correlation between enhanced structural order and reduced molecular mobility. It is consistent with the recent computational study of the mW and TIP4P/2005 water models, which demonstrated that enhanced structural ordering and ice nuclei formation are preceded by a reduction in local molecular mobility.⁹² Evidence of low molecular mobility is also observed in n_sw in sub-box No. 4, which is adjacent to sub-box No. 1, and in sub-box No. 1 at earlier times t < r₁. We posit that the reduction in molecular mobility may be a prerequisite for the development of structurally ordered regions, and that the formation of a domain with enhanced structural order requires a large region of reduced mobility where it can be fostered and accommodated. Although S in sub-box Nos. 2 and 4 does not increase over the time interval $t=0.50–0.80μ$s, large fluctuations in n_sw are observed in these sub-boxes during this period. Given that sub-box Nos. 2 and 4 are both adjacent to sub-box No. 1 via the periodic boundary conditions, this suggests that the low mobility region in sub-box No. 1 attempts to move or expand into these regions.

In time window R₂ (⁠$t=0.80–1.0μ$s), which coincides with the formation and growth of the nucleation seed, we observe marked increases and decreases in S and n_sw, respectively, in sub-box Nos. 1 and 4. This behavior further points toward a strong correlation between structural order and molecular mobility. The increases in S in sub-box Nos. 1 and 4 suggest that the nucleation seed partially occupies these two regions. The seed, however, which contains $≈300$ water molecules, is too small to completely span the two boxes. Thus, the reduction in molecular mobility in sub-box Nos. 1 and 4 suggests that the seed is surrounded by slowly moving liquid-like molecules in these two regions.

Inspection of the time evolution of the density profile along the major axis (z-axis) of the simulation cell reveals that the system is in a phase separated state in time windows R₁ and R₂ (Fig. 2). Indeed, as reported in Ref. 24, spontaneous HDL–LDL phase separation occurs rapidly after quenching the system below T_LLCP at t = 0 µs, on a time scale approximately at least an order of magnitude smaller than that relevant to crystallization. Moreover, the time-dependent density profile indicates that the LDL phase primarily occupies sub-box Nos. 1 and 4 in time windows R₁ and R₂. Thus, the analysis suggests that the formation and growth of the nucleation seed occur within the LDL phase.

Next, we examine the aggregation of ice-like molecules and polymorph selection. To characterize the aggregation of ice-like molecules, we computed the pair correlation function g(r) between Ic–Ic and Ih–Ih species in different time windows (Fig. 3). In calculating the correlation functions, we label ice-like molecule j as Ic if S_Ic(j) > S_Ih(j) and as Ih otherwise; molecules with ice-0 motifs are neglected in this analysis because their extremely low prevalence frustrates proper statistical sampling of correlation functions. For t ≤ 0.6 µs, the Ic–Ic and Ih–Ih g(r) functions mostly overlap, indicating competition between Ic- and Ih-like motifs in the early stage of the simulation [Figs. 3(a)–3(c)]. This competition can be understood by examining the positions of the oxygens in the second neighbor shells of Ic and Ih (Fig. 4). As discussed in Refs. 84 and 86, the 12 water oxygens in the second shell of Ic sit at the vertices of a cuboctahedron, whereas in Ih, they form a anticuboctahedron. The two geometries are found in fcc and hcp crystals,⁸⁴ respectively, and differ by a rotation of one of the two base triangles by 60°.⁸⁴ By contrast, the 12 oxygens in the second neighbor shell of supercooled liquid water are arranged in open cages, similar to those in amorphous ices,⁸⁶ which can close to form either Ic or Ih, generating competition between the two polymorphs. Here, we observe that the fcc and hcp structures of Ic and Ih, respectively, are connected by a smooth geometrical transformation, in which the icosahedron structure of ice-0⁷³ naturally arises as an intermediate state (Fig. 4), suggesting that icosahedral-like structures may be present in liquid water.⁹³ 

In the time window t = [0.60–0.65] μs, which encompasses point r₁ (t ≈ 0.62 µs) marking the initial slight increase in S [Fig. 1(a)], differences begin to appear in the Ic–Ic and Ih–Ih g(r) functions [Fig. 3(d)]. The Ic–Ic g(r) function subsequently develops a weak shoulder on the second peak and hints of third and fourth peaks at larger distances {Fig. 3(e), t = [0.65–0.70] μs}. The appearance of a weak shoulder on the second peak of g(r) is a general precursor to freezing transitions.^94,95 Furthermore, the emergence of these features in the Ic–Ic g(r) indicates that the increase in S in this time window is not only related to the increase in f_ice [Fig. 1(b)] but also related to the development of long-range correlations that temporally coincide with the onset of reduced molecular mobility, prior to the appearance of the nucleation seed. As time progresses, the Ic–Ic g(r) function develops signatures of increased structuring {Fig. 3(f), t = [0.70–0.75] μs}, and the nucleation seed eventually forms {Fig. 3(g), t = [0.75–0.80] μs}. By contrast, during this time, the Ih–Ih g(r) function does not develop a shoulder on the second peak or peaks at larger distances. The Ih–Ih g(r) function also remains largely unchanged as the nucleation seed begins to stabilize {Fig. 3(h), t = [0.80–0.85] μs}. These observations suggest that ice Ic is the dominant polymorph involved in the formation of the nucleation seed. Indeed, it is not until the later stages of seed stabilization that a weak shoulder emerges on the second peak of Ih–Ih g(r) {Fig. 3(i), t = [0.85–0.90] μs}, and the shoulder does not fully develop until the initial stage of seed growth {Fig. 3(j), t = [0.90–0.95] μs}. During expansion and growth of the seed [Figs. 3(k) and 3(l), respectively; t > 0.95 µs], the ice Ic–Ic and Ih–Ih g(r) functions again become nearly indistinguishable, suggesting that the crystallized portion of the sample is a mixture of the two polymorphs. It also indicates that growth proceeds by following the symmetry of the adsorbed crystallites, as interconversion between Ic and Ih would require rotation of the lower tetrahedron (Fig. 4), which would involve an extensive rearrangement of the hydrogen bond network and, therefore, is unfavorable.

To gain additional insight into the nucleation process, we analyzed the size of largest Ic and Ih clusters n_max, and the number of each cluster type n_cl, as a function of time over the period $t=0.25–1.25μ$s (Fig. 5). Two molecules i and j were defined to belong to the same cluster if they were nearest neighbors and classified as the same ice polymorph type. For $t=0.25–0.60μ$s, the clusters consist of mostly monomers and dimers dispersed in the bulk liquid (Fig. 5), consistent with the absence of the second-peak shoulders and long-range order in the g(r) functions during this period [Figs. 3(a)–3(c)]. At t ≈ 0.62 µs (point r₁), we observe an increase in n_max and a slight decrease in n_cl for Ic, which coincide with the development of the structural features in Ic–Ic g(r) [Figs. 3(d)–3(f)]. The crystalline clusters formed during this period reflect attempts of the system to generate a stable nucleation seed, but they are too small to stabilize and grow. We posit that structural reorganizations occurring at the surface of these small crystallites propagate inwards, destabilizing their cores and ultimately causing them to melt.

As indicated by the analysis above, the system attempts to form a nucleation seed at t ≈ 0.62 µs (point r₁) in sub-box No. 1 of the simulation cell. We investigate this process in detail by analyzing the fraction of ice-like neighbors in the first and second neighbor shells of each molecule. For molecule j, we consider its 16 nearest neighbors and compute $fj,icenn=nj,icenn/16$⁠, where $nj,icenn$ is the number of ice-like neighbors. Because the LOM S(j), used to identify ice-like molecules, characterizes order using both first and second nearest neighbors, the quantity $fj,icenn$ is sensitive to the local structure extending up to the fourth neighbor shell. Thus, $fj,icenn=1$ indicates that molecule j is in the core of a crystallite and surrounded by ice-like neighbors up to its fourth neighbor shell.

The spatial distribution of $fj,icenn$ at times shortly after t ≈ 0.62 µs (point r₁) in sub-box No. 1 reveals the existence of two ice-rich domains that are separated by liquid regions characterized by low molecular mobility [Fig. 6(a), t ≈ 0.659 µs]. Each ice-rich domain contains several crystallites of different sizes, which are in turn separated by slowly moving liquid water. Thus, the nucleation seed formation kinetics are governed by the dynamics in low-mobility liquid domains that contain dispersed crystallites. As time progresses, one of the ice-rich domains rotates toward the other. This process does not involve rotation of the domain as a whole, but rather it occurs through melting and formation of crystallites at the peripheries of the domain [Fig. 6(b), t ≈ 0.660 µs]. The rotation reduces the distance between the two ice-rich domains such that they are able to bridge via the formation of new crystallites [Figs. 6(c) and 6(d), t ≈ 0.661 µs and t ≈ 0.662 µs, respectively]. The bridged domains serve as the incubator for the future nucleation seed. The total size of the two domains fluctuates, as do the sizes of the individual crystallites contained within, but the domains do not disappear completely because they are each stabilized by their neighbor.

In time window R₂ (⁠$t=0.80–1.0μ$s), we observe that low-mobility water molecules located between the crystallites arrange in configurations that favor percolation, leading to the formation of a stable nucleation seed that almost entirely consists of the Ic polymorph. During this period, the number of clusters n_cl of Ic decreases and the size of the largest cluster n_max of Ic increases to $≈300$ water molecules, whereas n_cl for Ih remains almost constant (Fig. 5). The size of the stable seed and the number of Ic clusters fluctuate around a constant value for $≈0.1μ$s, indicating that the seed does not grow during this time. The growth that follows, however, occurs rapidly. As growth occurs, n_cl and n_max for Ic further decrease and increase, respectively (Fig. 5). This behavior indicates that the nucleation seed grows via the formation of new Ic sites at the seed–liquid interface and via adsorption of dispersed Ic clusters in the liquid at the seed’s surface. The increases in n_cl and n_max for Ih during the growth phase suggest that the formed crystal contains stacking faults. Indeed, the presence of stacking faults is confirmed by visual inspection of the system (Fig. 7) and is in agreement with the observations presented in previous computational studies of ice nucleation.^20,38,96 Additionally, during the early stages of growth, we observe the presence of a few molecules with ice-0-like local order on the surface of the nucleation seed. Ice-0 structures are identified via the LOM using the criteria S_I0(j) > S_Ic/h(j) > 0.8. The presence of ice-0 like structures on the seed’s surface is not surprising because the second neighbor shell of ice-0 is an intermediate arrangement between Ic and Ih (Fig. 4), and it results in a lower seed–liquid interfacial free energy,⁹⁷ which stabilizes the burgeoning crystal.

Prior studies have rationalized the propensity of water to form metastable ice Ic over the thermodynamically stable Ih polymorph by invoking the Ostwald step rules^43,44 and geometric arguments.²⁰ In particular, Haji-Akbari and Debenedetti proposed that Ic is favored because cubic structural motifs, similar to those posited to play a role in the homogeneous crystallization of silicon,⁹⁸ lead to crystallites that are more compact than those with hexagonal motifs. In accord with Ref. 20, our findings support the kinetic nature of this polymorphic selection, but we invoke a different interpretation to explain the prevalence of Ic crystallites in the liquid phase. Given the small free energy difference between bulk Ih and Ic⁷² and their identical tetrahedral configurations of first nearest neighbors, we posit that the polymorph selection originates in the second coordination shell. As discussed above, and as partially reported in Refs. 84 and 86, the second neighbor shells of Ic and Ih form cuboctahedron and anticuboctahedron arrangements, respectively.⁸⁶ Whereas the anticuboctahedron in Ih has one symmetry plane, the cuboctahedron in Ic has two identical symmetry planes (Fig. 8). Thus, in a regime of reduced molecular mobility, it is more likely for the liquid to close its second neighbor shell by forming the cuboctahedron motif found in Ic. This hypothesis is consistent with the argument proposed by Russo and Tanaka,⁹⁹ who showed that ice Ic is slightly entropically favored over Ih, at least at the level of the two-body excess entropy.

We performed unbiased classical molecular dynamics simulations of the ST2 water model under deeply supercooled conditions to investigate the nucleation of ice. Our analysis reveals that the system undergoes many attempts to generate a nucleation seed before crystallization occurs. The successful nucleation seed eventually emerges in small percolated regions of enhanced structural order that are embedded within a larger low-density liquid (LDL) domain characterized by low molecular mobility.

The nucleation seed formed in this region exhibits a period of transient stability before growth begins. Growth occurs via freezing of water molecules in contact with the surface of the seed and through the adsorption of ice Ic and Ih crystallites in the liquid at the seed’s surface. On the surface of the nucleation seed, we also observed molecules with ice-0-like local structure, which subsequently transform into Ih during growth. We rationalize this finding by observing that the molecular arrangement in ice-0’s second neighbor shell is an intermediate closed structure between that of Ic and Ih, and it facilitates nucleation by reducing the seed-liquid interfacial free. The final crystallized sample that forms is a stacking fault ice that is rich in the cubic polymorph. Although the HDL–LDL phase-separated system examined here is heterogeneous, our analysis suggests that nucleation occurs in the LDL phase and results in a stacking fault ice similar to those observed in previous studies of homogeneous ice nucleation. Thus, our findings intimate that nucleation may effectively occur homogeneously in the LDL phase, but this scenario requires further scrutiny in future studies.

Finally, our results suggest that polymorphic selection in favor of cubic ice (Ic) occurs in the liquid phase prior to the formation of the nucleation seed. We explain this observation by positing that the likelihood for the liquid to form Ic is higher because of its cuboctahedron second neighbor arrangement. In contrast with the second neighbor anticuboctahedron arrangement found in hexagonal ice (Ih), which contains only a single symmetry plane, ice Ic’s cuboctahedron structure has two identical symmetry planes, making it easier to form in liquid regions with low molecular mobility and, thus, kinetically favored.

This work was supported by the Hartree National Centre for Digital Innovation, a collaboration between STFC and IBM. J.C.P. acknowledges support from the Welch Foundation (Grant No. E-1882) and the National Science Foundation (Grant No. CBET-1751173).

The authors have no conflicts to disclose.

The data that support the findings of this study are available from the corresponding author upon reasonable request.