Surface-induced water crystallisation driven by precursors formed in negative pressure regions

Domain: earth_science
Published: 2024-07-26
Editorial outcome: major_revision

Abstract
Ice nucleation is a crucial process in nature and industries; however, the role of the free surface of water in this process remains unclear. To address this, we investigate the microscopic freezing process using brute-force molecular dynamics simulations. We discover that the free surface assists ice nucleation through an unexpected mechanism. The surface-induced negative pressure enhances the formation of local structures with a ring topology characteristic of Ice 0-like symmetry, promoting ice nucleation despite the symmetry differing from ordinary ice crystals. Unlike substrate-induced nucleation via water-solid interactions that occurs directly on the surface, this negative-pressure-induced mechanism promotes ice nucleation slightly inward the surface. Our findings provide a molecular-level understanding of the mechanism and pathway behind free-surface-induced ice formation, resolving the longstanding debate. The implications of our discoveries are of substantial importance in areas such as cloud formation, food technology, and other fields where ice nucleation plays a pivotal role.

Introduction
Water freezing phenomenon is ubiquitous in nature, science and technology, such as atmospheric clouds, biological cells, meteors and even air-conditioner applications. For instance, the radiative properties of clouds are determined by the formation of ice particles. In the real world, most water exists in forms with interfaces, especially the free surface, such as water droplets in clouds and on all kinds of objects in daily life. Therefore, understanding the role of the free surface in ice nucleation and the microscopic ice nucleation mechanism at the molecular level is crucial from both fundamental and application perspectives. However, it is still challenging to directly detect the ice nucleation events at their inception since the necessary temporal and spatial resolution still need to be improved in experiments. Tabazadeh et al. proposed surface-induced crystallisation based on classical nucleation theory. Nucleation events near the free water surface are favoured thermodynamically than in bulk if \sigma_{sv} - \sigma_{lv} > \sigma_{ls}, with \sigma_{sv}, \sigma_{lv}, and \sigma_{ls}, the solid-liquid, solid-vapour, and liquid-vapour interface tensions, respectively. Due to the challenges of measuring the interface tension, it is difficult to confirm whether the criterion works for water. Furthermore, in later years, more efforts on surface freezing provide results consistent and inconsistent with this theory, resulting in controversy until now. More seriously, whether we can use the macroscopic equilibrium properties to predict the nonequilibrium microscopic process, such as crystal nucleation, is unclear.

For addressing the problems of surface freezing in water, molecular dynamics (MD) simulations are powerful alternatives that make it possible to access microscopic information with a molecular-level resolution. However, the computational studies have also been inconclusive due to their dependence on the force fields, system sizes, and even ice molecule identification. For instance, it was reported that for a free-standing water film with a six-site force field on the nanoscale, ice crystals likely nucleate near the surface due to the lack of electrostatic neutrality. However, this result has not been reproduced in larger systems and is now speculated due to the finite-system-size effects, although this needs to be confirmed. For all-atom models incorporating charges, such as TIP4P and TIP4P/ice, exploring the statistical aspects of ice nucleation through brute-force simulations is currently a formidable task. This is primarily due to the difficulty in acquiring trajectories that are both long enough and numerous enough to provide meaningful insights.

Concerning these difficulties, the coarse-grained monoatomic water model, mW water, significantly improves the computational efficiency due to the absence of long-range Coulomb electrostatic interactions. Besides, the mW water model can also achieve remarkable success in reproducing the features of water’s energetics, thermodynamics, and structures. For example, the mW model well captures the water’s locally favoured tetrahedral arrangements and the relevant ice forms at ambient pressure, such as hexagonal Ice, cubic Ice, and Ice 0. Haji-Akbari et al. and Li et al. estimated the ice nucleation rates in free-standing films and nanodroplets of the mW water with a forward-flux sampling method. Both groups found that despite the system satisfying the criterion, \sigma_{sv} - \sigma_{lv} > \sigma_{ls}, surface freezing is unfavoured in the mW water, which is against Tanazadeh’s theory. Recently, Hayton et al. also reached the same conclusion using a seeding technique. Specifically, they found that nucleation remains bulk-like in films thicker than ~3.5 nm, but the nucleation rate starts to decrease below a thickness of ~3.5 nm at 220 K.

Controversially, the conclusion drawn from the investigation of the realistic water model, TIP4P/ice, using forward-flux sampling, starkly contradicts the conclusions for the mW water model described above. This study reported a remarkable seven orders of magnitude enhancement in ice nucleation for 4-nm films compared to the bulk. However, perplexingly, their observations indicated that freezing initiates away from the air-water interface, introducing a counter-intuitive aspect to the results. On this issue, Hayton et al. claimed that the apparent discrepancy between the mW and TIP4P/ice water models can be resolved through consistent treatment of truncated interactions between homogeneous and inhomogeneous systems when using a model with a long-range potential, such as TIP4P/ice.

Consequently, the current understanding remains far from conclusive. In this context, we highlight a potential challenge in applying seeding techniques and enhanced sampling algorithms to study crystal nucleation phenomena in a heterogeneous system.

Estimating the crystal nucleation rate through seeding techniques in a heterogeneous system can be tricky, as it is influenced by the selection of the seed’s position and shape. On the other hand, techniques such as forward flux sampling and biased Monte Carlo heavily rely on the careful selection of a relevant order parameter to accurately describe the process of crystallisation.

In previous studies mentioned above, the local 6th-order bond orientational order parameter q_6 was employed for estimating the ice nucleation rate. However, a recent study on heterogeneous ice nucleation on a crystalline solid substrate, employing a forward flux sampling method, demonstrated that the commonly adopted size-based method using the order parameter q_6, inadequately describes heterogeneous ice nucleation. Instead, the study highlighted the critical importance of considering the surface-induced geometric anisotropy of ice nuclei. Furthermore, it was reported that for other tetrahedral liquids like carbon and silicon, the local q_6 is not a suitable order parameter for enhanced sampling techniques, with q_3 being more effective. These studies demonstrated that special care is necessary in selecting an appropriate order parameter when estimating the crystal nucleation rate, particularly when employing enhanced sampling techniques.

Given the unbiased simulation approach we employ to study ice nucleation, shape anisotropy may not pose a problem. However, a suitable order parameter remains crucial for accurately detecting crystal precursors and nuclei. Therefore, our aim is to identify an order parameter applicable to all tetrahedral liquids for this purpose. While the order parameter q_6 can detect the six-fold symmetry of stable ices at ambient pressure, such as cubic and hexagonal ices, the ice formation pathway is not necessarily through a direct path but can involve Ice 0-like structures that contain five- or seven-membered rings, as manifested by the Ostwald step rule. To detect these structures, we need a coarse-grained bond orientational order parameter with 12-fold symmetry, Q_{12}, capable of identifying all types of crystal structures in tetrahedral liquids. Indeed, we have confirmed that the local order parameter q_6 tends to over identify ice molecules in the liquid phase, whereas the coarse-grained order parameter utilising the 12th-order bond orientational order parameter, Q_{12}, correctly identifies ice molecules and even outperforms q_6. This issue may not significantly affect the estimation of nucleation rates in bulk systems. However, in situations where surfaces induce ordering with symmetries other than 6-fold, as in our case, utilising q_6 for biasing could lead to incorrect conclusions.

To directly explore the microscopic kinetic pathway of ice nucleation, we conduct a series of extensive brute-force molecular dynamics (MD) simulations utilising the coarse-grained mW water model. To avoid potential curvature-induced effects, our initial focus is on studying free-standing water films featuring zero-mean-curvature vapour-liquid interfaces. Our approach involves a comprehensive analysis, employing topological ring analysis and the high-degree bond-orientational order parameter Q_{12}. The Q_{12} order parameter can pickup any crystal forms of water, including Ice 0. These tools allow us to quantitatively assess the relative arrangement of neighbouring water molecules, enabling the identification of ice-crystal precursors, colloquially referred to as preorders. These preorders exhibit a liquid state but possess crystal-like local orientational order. The utilisation of these advanced techniques helped us overcome the limitations associated with standard crystalline ice identification methods, ensuring an unbiased identification of preorders—precursors—from which ice crystals are anticipated to form. Our results furnish robust numerical evidence supporting free-surface-assisted ice nucleation, occurring near the free surface but slightly inward. Ice nucleation starts from the formation of Ice 0-like precursors near the free surface, but the precursors quickly transform more stable ordinary ices. Subsequently, we unravel an unconventional microscopic physical mechanism and kinetic pathway that facilitates ice nucleation in close proximity to the free surface through ice precursors. Finally, we demonstrate the relevance of our microscopic ice nucleation mechanism by extending our investigation to nanodroplets of the mW water model and thin films of the realistic TIP4P/ice water model.

Results
Preferential ice nucleation location in thin water films
First, we explain our strategy to attack a controversial and long-standing problem concerning whether a free water surface can facilitate ice nucleation or not. In earlier works, the ice nucleation rates of free-standing water films were estimated using enhanced sampling methods, such as forward-flux sampling and umbrella sampling methods. The cumulative transition probabilities were used to quantify the nucleation rates in both thin films and bulk water. However, the properties of the free surface are pretty complicated; for example, the surface structure is strongly distorted due to the symmetry breaking at the interface. Accordingly, the dynamics near the free surface is also different from the bulk. Hence, it has been considered reasonable if the ice nucleation rate is slower for a free-standing film than the bulk.

Instead, we ask whether a free water surface can induce nucleation and whether a preferential location of ice nucleation exists. To answer these questions, we directly detect the birth location of ice nuclei during the nucleation process by brute-force MD simulations. It shows that the ice nuclei are still small initially and fluctuate across the water film. As the largest nucleus grows and reaches the critical nucleus size, the first ice crystal is born in the sub-surface region and then increases its size linearly with time. We also show the microscopic dynamic process of ice nucleation in the water film. It is suggestive of free-surface-assisted ice nucleation. We also show the histogram of ice crystals of different sizes as a function of the distance from the film centre. We observe that ice crystals are more likely to nucleate in the sub-surface region (not directly on the surfaces) than in the middle region, based on statistics from over 100 independent trajectories. The symmetric double-peak distribution of the ice nucleation probability (predominantly cubic ices) near the two surfaces of the water film strongly supports free-surface-induced ice nucleation.

Moreover, we directly measure the nucleation rates of both water films and bulk water. Nucleation is inherently a stochastic process; therefore, even for the same system, we must wait for varying time intervals to observe nucleation events. In this context, we determine the birth time of the first crystal, referred to as the induction time t_{ind}, for each trajectory by monitoring the time evolution of potential energy. Utilising the distribution of induction times t_{ind}, we estimate the probability for a single simulation trajectory to remain in the liquid state at time t, denoted as P_{liq}(t). The results are shown in Supplementary Fig. 2b (see “Methods” for the definition). We quantify the nucleation rate R by fitting a stretched exponential function to P_{liq}(t) (see “Methods” for the details).

In Fig. 1f, we illustrate the comparison of nucleation rates per volume, R/V (nm^{-3}s^{-1}), between bulk water and water films. Our calculated nucleation rate for bulk water aligns well with the results reported in refs. 25, 38. Notably, the nucleation rate in the water film surpasses that of the bulk, exhibiting an increasing trend as the water film thickness decreases. This observation aligns with our visual evidence that ice crystals tend to form predominantly in the sub-surface region of a thin water film. Here, it is worth noting that in thin films, crystal nucleation is influenced by dynamic capillary waves, a factor absent in quiescent bulk states. We plan to investigate the impact of this effect on crystal nucleation and growth in future studies.

Our findings diverge qualitatively from earlier results obtained using the same mW water model. In 2003, Lü et al. demonstrated that the total nucleation rate of water films increases and approaches the bulk one when the thickness of films increases. They reported that the ice nucleation rate decreases by only about a factor of two for the thinnest film (4.72 nm) compared to the bulk. On the contrary, Haji-Akbari et al. reported a significantly more pronounced decrease in the nucleation rate by two to three orders of magnitude in 5-nm-thick films and a decline of seven orders of magnitude for 2.5-nm-thick films. It is important to note that Lü et al. and Haji-Akbari et al. employed different enhanced sampling methods, namely, the mean first-passage time vs. forward-flux sampling, and different order parameters to identify ice molecules, i.e., the CHILL algorithm vs. q_6, respectively. Additionally, a recent study by Hayton et al., utilising a seeding technique, further supported the notion that the nucleation rate is lower for thin films compared to the bulk.

Now, we consider the reasons for the disparity in conclusions between our study and those mentioned above. As discussed earlier, the sampling methods are highly sensitive to the choice of order parameters associated with crystals. The previous studies mentioned above utilised the CHILL algorithm or q_6 order parameters, assuming the relevant ice type to be ice I, i.e., cubic or hexagonal ices. However, as we will demonstrate later, this assumption might not be valid. For tetrahedral liquids, such as water, the most relevant order parameter is Q_{12}, involving its 16 neighbours, as it can detect any crystals favoured by local tetrahedral symmetry, including Ice 0 and Clathrate. We employ this order parameter to identify ice-like preorders and ice crystals.

Moreover, it is crucial to note that while seeding techniques may prove valuable for examining crystal nucleation behaviours, their efficacy in investigating the influence of the free surface on ice nucleation could be questionable. Ice nucleation near a free surface may follow a distinct kinetic pathway involving Ice 0-like precursors, adding complexity to the process. Unlike in bulk, the selection of the seed’s position and shape may also play a crucial role in estimating the nucleation rate in the presence of symmetry-breaking surface fields.

Taking these factors into consideration, we assert that a direct approach—employing brute-force simulations—provides a more thorough and convincing estimation of the crystal nucleation rate in free-standing films, compared to enhanced sampling algorithms or seeding methods, even though this method demands more computation time.

Surface-induced nucleation mechanism
Due to its unique structural and mechanical properties, the water surface plays a crucial role in ice nucleation. For instance, Odendahl et al. recently reported that water and ice share similar characteristic structural features, implying the presence of nanometre-scale ice-like domains at the air-water interface. In our work, we also observe a layering phenomenon and ordered structure near the free water surface in Fig. 2a, suggesting a structured liquid at the surface. Upon closer examination of the first few layers of the water surface, we identify an ordered local structure (red region) in the 2D plane of the well-defined outermost two molecular layers, composed of five-, six-, and seven-membered rings. A recent experimental study of 2D ice crystals formed on a solid substrate using scanning tunnelling microscopy (STM) indicates that the preordered structure of 2D hexagonal ice consists of abundant five-, six-, and seven-membered rings. This implies that the water surface possesses a 2D hexagonal ice-like local structure. Moreover, by comparing the angular distribution of a water molecule with its neighbours in the 2D plane parallel to the water surface, we also observe pronounced peaks in the angular distribution of water molecules in the first two layers, s1 and s2, similar to the 2D hexagonal ice. The peak position in the three-body-angle distribution is at 120°. Therefore, both previous results and our observations seem to indicate that the ice-like ordered water surface may aid in facilitating crystallisation, contrary to the notion of ‘bad’ ordering suppressing nucleation reported in previous works.

Consequently, assuming the presence of surface-induced precursors on the water surface appears to be physically reasonable. We use the 2D orientational order parameter \lambda_1 (see “Methods” for the definition) to quantify the surface structural order and its correlation with ice nucleation in the sub-surface region seen in Fig. 1. As shown in Supplementary Fig. 3a and b, we also find highly ordered surface molecules (indicated by green beads), which strongly correlate with the ice precursors and ice clusters formed in the interior yet close to the surface. This result is consistent with the scenario of surface-induced ice nucleation in free-standing films. However, it is essential to note that ice nuclei are formed predominantly in the sub-surface regions, not directly on the free surfaces, as depicted in Fig. 1d and e. Hence, the critical question arises: whether the 2D ice-like order of the surfaces is responsible for ice nucleation, or if some other source plays a more dominant role in inducing ice nucleation in the sub-surface regions. Below, we will show that the latter is the case.

To address this question, analysing water structures in the sub-surface regions is crucial. In this context, we opt for ring analysis as an alternative to the conventional bond-orientational order parameter. This choice is motivated by the surface’s potential to distort the symmetry of locally ordered structures in its vicinity, a phenomenon absent in the bulk. A striking revelation emerges from our analysis: 5-membered rings are significantly more abundant at the surface than in the bulk and extend deep into the interior of the water film, reaching a characteristic length of 1.5 nm, as illustrated in Fig. 2d. The 5-membered ring exhibits a preference in energy over the 4-membered ring and a greater favorability in entropy than the 6-membered ring. This competition, coupled with the surface constraint, serves as the origin for the enrichment of 5-membered rings at the surface.

Russo et al. reported the spontaneous creation of structural orders associated with a metastable ice form known as Ice 0, containing 5-, 6-, and 7-membered rings, in a supercooled state. These structures act as ice precursors by lowering the free-energy barrier of ice nucleation. Ice 0 structures exhibit two distinct local structural motifs: (5,5,5,5,6,6) and (5,5,5,7,6,6), as depicted in Fig. 2e, where the numbers 5, 6, and 7 indicate the count of water molecules in rings including the central molecule as a member. It was shown that Ice 0-like structures (see Fig. 2e) play a crucial role in the ice nucleation of bulk water. These structures are more readily formed in supercooled water than Ice I-type orders (cubic and hexagonal arrangements) since their ring distribution closely resembles that of liquid water. A molecule is classified as Ice 0-like if at least five of its six rings closely match the Ice 0 ring structures illustrated in Fig. 2e. These Ice 0-like preorders play a pivotal role in fostering the nucleation of ice crystals by substantially diminishing the crystal-liquid interfacial energy cost. Ice 0-like structures are not stable, and as the order develops, they undergo rapid transformation into ordinary ices, predominantly in cubic ices, without forming Ice 0.

Interestingly, we observe an abundance of Ice 0-like preordered structures in the sub-surface regions. In Fig. 2f, the density of Ice 0-like precursors, \rho_{I0-like}, in the pink-shaded region exhibits a prominent maximum near the position adjacent to the 2D ice-like surface layers (see the blue-shaded region). Notably, the characteristic thickness of this sub-surface region is ~1.5 nm, as observed in Fig. 2f. The pressure within the water film, both in the normal and lateral directions, is negative in these sub-surface regions (the pink-shaded regions). It is known that for the bulk case, the fraction of Ice 0-like molecules increases as pressure decreases (see Supplementary Fig. 4). In connection with this, it is worth noting that the homogeneous ice nucleation rate is enhanced under negative pressures. Hence, we deduce that the abundant Ice 0-like environments stem from the surface-induced negative pressure in the sub-surface regions. Ice nuclei emerge in these preordered liquid environments, specifically in highly concentrated Ice 0-like environments, as illustrated in Fig. 2g. This surface-induced preordering in liquid water provides a free-energy-based rationale for surface-induced ice nucleation. Interestingly, a similar phenomenon of surface-induced nucleation in a negative pressure region was reported for silicon, suggesting a general mechanism of surface-induced nucleation for tetrahedral liquids showing a density decrease upon solidification.

The above finding indicates that the primary role in surface-induced ice nucleation is played by Ice 0-like preordering in the sub-surface regions (the pink-shaded region in Fig. 2f) rather than 2D ice-like order on the surfaces (the blue-shaded region in Fig. 2f). Symmetry breaking at the surfaces not only induces layering, i.e., short-range translational order, leading to 2D ice-like order, but also stretches the surfaces, resulting in negative pressure in the sub-surface regions. While the former contributes to the induction of 2D ice-like order, the latter is more influential in enhancing Ice 0-like preorders. Our study suggests that the latter mechanism plays a more substantial role than the former in the formation of ice nuclei in water films. This observation may be attributed to the fact that, unlike the case of a solid substrate, the 2D ice order at the surface is not as robust due to surface fluctuations, i.e., capillary waves, which strongly perturb it.

Besides thermodynamic factors, dynamics is another crucial factor controlling crystal nucleation. The dynamics is also heterogeneously distributed in water films. Specifically, the surface mobility of the mW water is slower than that of the interior (see Supplementary Fig. 5). For bulk water, Martin et al. proposed that ice nucleation in supercooled water occurs in less mobile liquid regions. This phenomenon can be explained through two steps. Firstly, preordered regions exhibit lower dynamics in a supercooled liquid, and secondly, preorders facilitate crystal nucleation because of the angular symmetry matching that lowers the energy cost of the crystal-liquid interface. The slower diffusion near the water surface suggested ordered surface structures, as shown in Supplementary Fig. 5b, providing additional support for the surface-induced ice nucleation mechanism described above.

Ice nucleation in a water nanodroplet
A water droplet is one of the most common forms of liquid water in nature and is crucial for atmospheric science, such as clouds and rains. It is also of technological importance. A water droplet, characterised by a large surface-to-volume ratio, is widely employed to investigate the significance of the surface on water’s properties, including its impact on dynamic behaviours and ice nucleation.

To further confirm the effects of the free surface on nucleation, we directly measure the spatial distribution of Ice 0-like and ice molecules and track their evolutions in water droplets. Remarkably, nucleation events in water droplets are not random; instead, ice nuclei exhibit a higher probability of occurring in positions close to the surface in nanodroplets, as illustrated in Fig. 3a (see also Supplementary Movie 2). This behaviour is evident when projecting the positions of molecules belonging to ice precursors and crystals onto a 2D cross-section map, as shown in Fig. 3b. The distribution of ice molecules along the radial direction, depicted in Fig. 3c, aligns with the 2D spatial map of ice molecules (see the inset of Fig. 3c). We also provide a microscopic view of the dynamic ice nucleation process in the water droplet at 150 K in Supplementary Movie 2. This result strongly supports the surface-assisted nucleation mechanism for water droplets, a topic that has been debated for decades. Based on Fig. 3a and b, we can conclude that ice nucleation in water droplets occurs heterogeneously in space, facilitated by the curved surface. Our structural analysis further reveals the formation of Ice 0-like precursors near the surface, as evidenced by the distributions of ring structures and Ice 0-like structures along the droplet’s radial direction, as shown in Fig. 3d and e, respectively.

In water droplets, Ice 0-like structures are abundant in the sub-surface region and gradually decrease into the interior. Similar to water films, we observe a negative tangential pressure P_T near the water droplet surface due to surface tension, as reported by Li et al. We conclude that this surface-induced negative pressure is the physical reason for the rich Ice 0-like environments close to the droplet surface. As well-known, the surface curvature of a water droplet creates extra pressure called the Laplace pressure, P_L = P_N - P_T, which increases the pressure inside the droplet. As the droplet size decreases, the internal pressure increases (P_L), which tends to suppress the probability of a nucleation event in the droplet, as shown in Li et al.’s work. However, due to the relatively low tangential pressure (P_T) close to the surface, ice nucleation still prefers to occur near the surface for the droplet sizes studied here.

Comparison with TIP4P/Ice water model
As described above, the use of a coarse-grained water model (mW water in our case) provides us with the capability to directly observe the microscopic dynamics of ice nucleation, facilitated by its computational efficiency. However, there might be concerns regarding the reliability of this model in faithfully representing water properties, given that simulation outcomes can be influenced by the choice of force field and system size.

In this regard, the TIP4P/ice model stands out as one of the most accurate classical molecular models of water, particularly in describing ice formation. Nevertheless, the direct simulation of ice nucleation events using this water model proves challenging due to its computationally demanding nature. Recently, it has been demonstrated that crystal nucleation is initiated by the spontaneous formation of crystal precursors (or preordering) within a liquid state. Consequently, our specific focus is directed towards exploring preordering rather than directly simulating nucleation events. To identify preorders, we employ ring analysis to capture the water structures near free surfaces, with the expectation that these structures would exhibit an Ice 0-like configuration, as observed in the results obtained with the mW water model. The ring and density distributions of Ice 0-like preordered structures are illustrated in Fig. 4a and b, respectively. Consistent with the mW water results (see Fig. 2), the fraction of the five-membered rings increases towards the surfaces, and the concentration of Ice 0-like structures, formed by hydrogen bonds, exhibits distinct maxima in the sub-surface regions.

It is reasonable to anticipate that ice crystals would more readily be nucleated in regions with a high degree of Ice 0-like preorder, as observed in the mW water case presented in this study and the bulk case. Thus, our TIP4P/ice model results provide robust support for surface-induced freezing in real water. The notion of ‘bad’ interfacial ordering in prior works might be a misguided conclusion stemming from the usage of stable bulk ice orders (Ice c and h) as the critical structural order parameters to describe ice nucleation. These structures, associated with stable ices, differ significantly from those of liquid water (whether in the bulk or near the surface). Due to the high interfacial energy cost arising from the substantial structural differences between the crystal and liquid phases, direct nucleation of ordinary ices is not favoured. Instead, Ice 0-like precursors, exhibiting a ring distribution similar to that of liquid water, are demonstrated to be preferentially formed near the free water surfaces, acting as triggers for ice nucleation. This two-step-like ordering pathway aligns with Ostwald’s step rule of phase ordering.

Methods
We perform molecular dynamics simulations using a large-scale atomic/molecular massively parallel simulator (LAMMPS) for both mW and TIP4P/ice systems. Bulk water and films & droplets are conducted in the NPT and NVT ensembles, respectively, with periodic boundary conditions. The Nose-Hoover thermostat and the Parrimello-Rahman barostat are applied to control the temperature and pressure, respectively. The dimensions of the simulation boxes along x, y, and z directions are 9.44 nm, 9.44 nm, and 21 nm, respectively. The timestep is 1 fs. 3D local structural order is analysed using the coarse-grained Steinhardt bond orientational order parameter Q_{12} and the CHILL algorithm, whereas 2D structural order on the free surface is analysed using quasi-2D order parameter \lambda_1.

Ice 0 has two different local environments around the central molecule: (5,5,5,5,6,6) and (5,5,5,7,6,6), as shown in Fig. 2e. A molecule is classified as Ice 0-like if at least 5 of the six rings around it with 16 neighbours within cutoff r_c = 0.488 nm have the correct number of the Ice 0 ring members. f_n indicates the fraction of the n-membered ring around a molecule.

Q_{12} and connections
To characterise bond orientational order of ice-like structures, we use the coarse-grained order parameter, Q_{12}, defined as follows,

q_{12,m}(i) = \frac{1}{N_b(i)} \sum_{j=1}^{N_b(i)} Y_{12,m}(r_{bij}), (1)

and then the averaged form is,

Q_{12,m}(i) = \frac{1}{N_{gb}(i)} \sum_{k=1}^{N_{gb}(i)} q_{12,m}(k): (2)

The scalar product between Q_{12,m} of two molecules is defined as S_{ij} = \sum_{m=-12}^{12} Q_{12,m}(i)Q^*_{12,m}(j). A pair of molecules i and j is identified as connected if S_{ij} is larger than 0.75. If the number of connections that the molecules have is above 75% N_b(i), the molecule is identified as crystalline.

The CHILL algorithm
The local order around each water molecule i is defined by a local orientational bond order parameter vector q_l with 2l+1 complex components,

q_{lm}(i) = \frac{1}{4} \sum_{j=1}^{4} Y_{l,m}(r_{bij}), (3)

where q_{l,m}(i) projects the bond orientational order of the four closest neighbours of a molecule based on spherical harmonics. The alignment of the orientation of the local structures is measured by the normalised dot product of q_l between each pair of neighbour molecules:

c_{ij} = \frac{\sum_{m=-l}^{l} q_{lm}(i)q^*_{lm}(j)}{\sqrt{\sum_{m=-l}^{l} q_{lm}(i)q^*_{lm}(i)} \sqrt{\sum_{m=-l}^{l} q_{lm}(j)q^*_{lm}(j)}} (4)

The quasi-2D ice nucleation
We chose the order parameter \lambda_1(i) to identify the crystal particles. However, the 2nd layer is not an exact 2D but a quasi-2D with bilayer substructures.

\lambda(i) = \frac{1}{N_b(i)} \sum_{j=1}^{N_b(i)} \sum_{m=-6}^{6} q_{6,m}(i)q^*_{6,m}(j) (5)

\lambda_1(i) = \frac{1}{1+N_b(i)} (\frac{1}{1+N_b(i)} \lambda(i) + \sum_{j=1}^{N_b(i)} \lambda(j) + \sum_{k=1}^{N_b(j)} \lambda(k)) (6)

According to the distribution of \lambda_1(i) of interfacial liquid and ice, we can identify the ice molecules in the contact layer as follows.

\lambda_1(i) \le 0.23 liquid
0.23 < \lambda_1(i) \le 0.4 ice-preorder
\lambda_1(i) > 0.4 ice

Then, the ice molecules in the second layer are identified as

\lambda_1(i) \le 0.2 liquid
0.2 < \lambda_1(i) \le 0.26 ice-preorder
\lambda_1(i) > 0.26 ice

The estimation of the nucleation rate R
Here, we determine the waiting time, denoted as the induction time t_{ind}, for each trajectory through the time evolution of potential energy (refer to supplementary Fig. 2). Utilising the distribution of induction times for nucleation, we estimate the probability for a single simulation trajectory to remain in the liquid state at time t using the

P_{liq}(t) = 1 - \frac{1}{N_{sim}} \sum_{i=1}^{N_{sim}} \Theta(t - t_{ind,i}), (7)

where N_{sim} represents the total number of simulations conducted for each system (in this case, 50), t_{ind,i} is the induction time determined for the i-th simulation, and \Theta(x) is the Heaviside step function. Subsequently, applying Eq. (7), we derive P_{liq}(t). To assess the nucleation rate from the shape of P_{liq}(t), we employ a stretched exponential function:

P_{liq}(t) = exp[-(Rt)^{\gamma}], (8)

where R denotes the nucleation rate. The values of R are acquired by fitting P_{liq}(t).

Data availability
All study data are included within the article and Supplementary Information. Source data are provided with this paper.

Code availability
The simulation codes used in this study are available from the corresponding authors upon request.

Peer review rounds
[{"round": 1, "reviewer_comments": {"reviewer_1": "In this work, the authors carried out a molecular dynamics study of ice nucleation within a\nfreestanding water film. The study was motivated by the controversy regarding the role of water-air\nsurface on ice nucleation. Despite the proposed surface-induced crystallization, there is a clear lack of\nexperimental evidence, and results from molecular modeling are also inclusive. The authors employed\nthe mW water model, which has been used in prior modeling studies that did not find evidence of\nsurface freezing, to conduct MD simulation and reported the preference of formation of Ice 0 structure\nwithin the subsurface of water film. The preference was attempted to be rationalized by the authors\nthrough the abundance of topological fragments resembling ice 0 near surface induced by negative\npressure due to surface tension.\nDespite the significance of the topic itself, there are a few major scientific issues that make the paper\nnot suitable for publishing on Nature Communications. The most critical issue is that there is NO clear\nevidence of surface freezing in this study. The only “evidence” is the formation of ice 0 structure near\nsurface, which is far from being sufficient. Even in the case where water film yields a rate lower than\nor close to homogeneous rate, there was still a tendency for ice nucleation to occur near surfaces. In\nparticular, the thickness of water film in the study 8 nm is barely enough to differentiate bulk region\nfrom surface region. Finite size effect has been previously demonstrated to play an important role in\nwater films. How would the picture change with an increasing film thickness?\nMore relevant questions should be addressed to claim surface induced ice nucleation: Did water film\nyield a higher nucleation rate than bulk water of the same volume, or did water film crystallize at a\nstatistically higher temperature? Only a positive answer to either question can constitute a more\nmeaningful evidence of surface crystallization. The authors have already obtained 50 independent\ntrajectories which should allow them to extract unbiased ice nucleation rate (along with its statistical\nuncertainty) for water film. How does it compare with homogeneous ice nucleation at the same\ntemperature (206 K)? If the rates are indistinguishable, how would one conclude free surface\npromotes ice nucleation?\nAlternatively, they can also investigate the free energy barrier. The authors made a strong claim that\nthe order parameters used in previous studies did not capture the relevant pathway and argued the\nmost relevant order parameter should be Q12. Although this reviewer does not necessarily agree on\nthis, the authors can at least try to make it self-consistent by showing ice nucleation in subsurface\nindeed leads to a lower free energy barrier using the order parameter based on Q12. They do have all\nthe data to allow them to do so.\nOther than these major issues, the overall writing of the paper can be improved. Certain part of the\npaper is difficult to understand. For example, page 9-10 where the authors discussed the role of\nLaplace pressure and commented on previous studies. For those who may not know exactly these\nreferences, it’s rather difficult to understand the main points through these sentences because they\nare mixed with authors’ opinion. By the way, when claiming others’ conclusion “is not correct” multiple\ntimes in the paper, in what position do the authors believe they have to make such a strong\nstatement? Lastly, when trying to generalize their conclusion to other materials like tetrahedral liquids,\nthe authors should be aware of early studies that already attempted to tackle this problem, for\nexample, Nat. Mater. 8, 726 (2009); JCP 18, 4102 (2016).", "reviewer_2": "This study reports molecular simulations showing that ice tends to nucleate with greater probability at\na certain distance from the ‘free surface’ of supercooled water, and attributes this to surface induced\nnegative pressure favouring so-called Ice-0 like symmetry. It’s an interesting study however I have a\nfew major concerns and would not suggest publication of the manuscript in its current form. I have\ndetailed these concerns below.\nMajor concerns:\nA central claim of the paper is that ice nucleation is favoured some distance, roughly 2 nm, from the\nwater surface. I have a few issues with this. Firstly, it is not obvious to me that what is shown in Fig. 1\n(b) is necessarily statistically significant. I would suggest the authors find some way of comparing\ntheir results to the null hypothesis that nucleation takes place uniformly across the box of water. The\nasymmetry of the peaks at roughly 2 and -2 nm, suggests that the finding may not be robust.\nI find the idea of negative pressure in a droplet or film a bit slippery. Do you not need to consider the\ndirection of pressure? This said, wouldn’t you expect the negative pressure due to surface tension to\nbe at a maximum closer to the surface than 2 nm? See Sega et al. (J. Chem. Theory Comput. 2016,\n12, 9, 4509–4515). Similarly, if the claim is that enhancement in nucleation rate is due to negative\npressure wouldn’t it make sense to work out a pressure profile for the system and compare it to the\ncluster distribution profile across the system?\nI have a concern about the experimental testability and real-world relevance of the study. The\nintroduction motivates the work of the basis of relevance to various applied problems where ice\nformation matters yet there is no discussion of how these findings might help us understand real\nsystems. I think he question in the end is ‘do we expect a small water droplet, with a water-vapour\ninterface to freeze at a warmer or colder temperature than would be anticipated for an equivalent\nvolume of ‘bulk’ water?’. The work here would presumably predict that smaller droplets, containing a\ngreater proportion of ‘interfacial’ water would tend to have a higher nucleation rate per volume than\nlarger droplets? To the best of my knowledge, there is no experimental evidence for surface\nenhancement of apparently homogeneous ice nucleation (Atkinson et al. J. Phys. Chem. A 2016, 120,\n33, 6513–6520 and Duft and Leisner Atmos. Chem. Phys., 2004, 4, 1997–2000). I would note here\nthat some of the work cited (Shaw et al. J. Phys. Chem. B 2005109, 9865–9868) where nucleation is\nenhanced near a surface refers to heterogeneous ice nucleation and would seem to me to be of limited\nrelevance to the present study.\nIt is worth noting that Hayton et al. (DOI: 10.1039/D3FD00099K (Paper) Faraday Discuss., 2023)\nhave very recently attempted to assess how we might expect the thinness of water films to impact\ntheir ice nucleation rate compared to a more naively calculated homogeneous ice nucleation rate and\nfound that the difference is minimal. An approach of this type might help clarify the relevance of the\npresent study.\nFinally, I do not think that that the importance of identifying ice-0 like precursors is sufficiently\nwell articulated at present. I am all for fundamental science, but I think it important to lay out how\nknowing the structure of critical cluster in a given situation might have wider scientific impact.", "reviewer_3": "This manuscript investigates ice nucleation by means of molecular dynamics simulations. From a\nmethodological standpoint, the authors study ice nucleation in water films and droplets, make use of\ntwo water models, namely mW and TIP4P/ice, and also put the spotlight on the relevance of chosing\nsuitable order parameters.\nThe main result of the work consists in revealing a novel nucleation mechanism in which the free\nsurface-induced negative pressure promotes nucleation via the emergence of Ice 0-like precursors\nwith ring topologies containing 5 and 6-membered rings. Since these precursos are expected to be\nmore friendly to liquid water as comapred to standard (bulk) ice-like topologies, they are expected to\nlower the interfacial energy cost for ice nucleation. This finding would rectify previous statements\nregarding free surface order as “bad” for ice nucleation. Additionally, unlike substrate-induced ice\nnucleation that occurs directly at the outermost layer, here the process occurs slightly inwards (at the\nsub-surface layer, but very close to the surface). Hence, the authors not only reproduce the surface-\ninduced crystallisation observed in experiments, but they also unravel its physical molecular\nunderpinnings.\nFrom the above-expounded considerations I find this work to tackle a problem of utmost relevance,\nboth from fundamental and practical perspectives, and to present a solid, comprehensive and carefully\ndescribed approach that yields highly significant results. More importantly, the study provides a new\nconceptualization that proposes a solution to the core physical and molecular basis of the problem.\nThus, I am glad to recommend publication.\nThere is only one point I would like the authors to comment. The local orientational bond order,\nq_lm(i), is calculated by considering the four closest neighbours of a molecule. For the molecules at\nthe first (outermost) layer, where important 3-fold hydrogen bond coordination is expected (with the\npresence of dangling bonds), the fourth neighbor would lie at the second coordination shell instead of\nbelonging to the first-shell. Thus, unlike the quasi-2D lambda index or the study of five or six-\nmembered rings, this index might produce artifacts in such case. Please comment on the application\nof this metric in your study."}, "author_response": "Replies to the comments of Reviewer #1:\nIn this work, the authors carried out a molecular dynamics study of ice nucleation within a freestanding water film. The study was motivated by the controversy regarding the role of water-air surface on ice nucleation. Despite the proposed surface-induced crystallization, there is a clear lack of experimental evidence, and results from molecular modeling are also inclusive. The authors employed the mW water model, which has been used in prior modeling studies that did not find evidence of surface freezing, to conduct MD simulation and reported the preference of formation of Ice 0 structure within the subsurface of water film. The preference was attempted to be rationalized by the authors through the abundance of topological fragments resembling ice 0 near surface induced by negative pressure due to surface tension.\nFirst, we thank the reviewer for carefully reading our manuscript and providing comments. The feedback and guidance provided by the reviewer have been invaluable in enhancing the quality and credibility of our research.\nBefore going to the specific details, we would like to emphasise that our simulations represent brute-force investigations of direct ice nucleation processes. Our findings consistently demonstrate that ice formations always occur in proximity to the free surface, regardless of film thickness or geometry, encompassing thin films and nanodroplets. It is important to note that we do not employ any enhanced sampling techniques or seeding methods in our simulations. As a result, our conclusion that ice formations tend to occur near the surface rather than in the central region is based on empirical observations and does not depend on any underlying assumptions, including the presence of ice 0. We believe that this strengthens the validity and reliability of our findings.\nWe also directly measure the nucleation rates of water films and bulk. Nucleation is intrinsically a stochastic process; therefore, even for the same system, we must wait for different time intervals to observe nucleation. Thus, to estimate the nucleation rate R, we performed more than 50 simulations for each system under the same condition. Subsequently, we determined the crystal birth time, i.e., the induction time for crystal nucleation (t ), for each trajectory by the time evolution of the potential energy (see Fig. R1a). Based on the distribution of the induction times of the nucleation, we could estimate the probability for a single simulation trajectory to stay in the liquid state at time t:\n, (1)\nwhere N is the total number of simulations performed for each system, t i is the induction time determined for the i-th simulation, and Θ(x) is the Heaviside step function. Then, following Eq. (1), we calculated P (t).\nliq\nBy examining the shape of P (t), we can quantify the nucleation rate by fitting the stretched exponential function to P (t):\nliq\n(2)\nwhere R is the nucleation rate. The fitting results are displayed in Fig. R1b. To compare the nucleation rates of bulk water and water films, we illustrate the nucleation rate per volume, denoted as R/V (nm-3ns-1), in Fig. R1c. It is evident from the results that the nucleation rate of the water film is higher than that of the bulk at different temperatures, and furthermore, it increases as the water film becomes thinner.\nFigure R1. Estimation of the nucleation behaviour. The nucleation rates were directly measured by tracking the time evolution of potential energy for both bulk water and water films with thicknesses of L=8 nm, 6 nm and 4 nm at T=200 K and 206 K. a, The determination of the induction time (t ) by monitoring the time evolution of the potential energy. b, The probability of a system remaining in a liquid state as a function of time, P (t), for different conditions. c, The comparison of the nucleation rate per volume R/V between bulk water and water films with different thicknesses.\nIn the following, we provide individual responses to each of the reviewer's comments.\nDespite the significance of the topic itself, there are a few major scientific issues that make the paper not suitable for publishing on Nature Communications. The most critical issue is that there is NO clear evidence of surface freezing in this study. The only “evidence” is the formation of ice 0 structure near surface, which is far from being sufficient. Even in the case where water film yields a rate lower than or close to homogeneous rate, there was still a tendency for ice nucleation to occur near surfaces. In particular, the thickness of water film in the study 8 nm is barely enough to differentiate bulk region from surface region. Finite size effect has been previously demonstrated to play an important role in water films. How would the picture change with an increasing film thickness?\nWe respectfully disagree with the reviewer's assessment of the points raised. We are afraid that there might be some misunderstanding. The reviewer specifically noted our observation of “ice 0”. However, it is crucial to understand that metastable “ice 0”-like structures only exist transiently in the early stage of crystal nucleation. Ice 0-like structures formed near the free surface quickly transform into standard ice crystals, such as cubic ice. Thus, our observation primarily pertains to ice nucleation near the surface, a fact substantiated by the clear peak evident in the probability distribution function (see new Fig. 1b), which indicates a heightened likelihood of ice nucleation occurring near the surface.\nIt is also important to highlight that the identification of ice nuclei is directly carried out using the bond orientational order parameter Q , rather than “Ice 0”-like structures. Additionally, the direct measurement of the nucleation rate for both bulk and water films, as illustrated in Fig. R1 above, serves to confirm our observation regarding nucleation induced by the free surface.\nWe want to emphasise that our conclusion is drawn from direct observations in brute force simulations, eliminating reliance on assumptions. Our findings are rooted in thorough and extensive computational investigations, ensuring the reliability and robustness of our results.\nWhile it is true that our system size may be comparatively small, it is important to acknowledge that this size is widely accepted as standard practice within the field of computer simulations focused on ice nucleation (see, e.g., references [2,3] below). Notably, even renowned research groups led by world experts such as Professor Pablo Debenedetti and Professors Angelos Michaelides and Stephen J. Cox employed similar simulation setups in their investigations, with film thicknesses of 4 nm in the former [3] and ranging from 2.5 to 6 nm in the latter [2]. Moreover, we confirmed our conclusion by using a droplet with a diameter of 12 nm. It is worth highlighting that our film thickness, at 8 nm, exceeds those utilised in these seminal works.\nThe findings of both groups supported the notion that ice nucleation is more likely to occur in the central region of a thin water film. We attribute the discrepancy between the conclusion of Professor Debenedetti's group and ours to a critical factor - the selection of bond-orientational order parameters. We have elaborated on this aspect in detail in our manuscript, explaining how the choice of these parameters may have led to differing conclusions between our study and Professor Debenedetti's work.\nOn the other hand, the differences in the conclusions between our study and the research conducted by Professors Michaelides and Cox can be attributed to the utilisation of a seeding technique. Specifically, they employed a method that involved seeding a spherical nucleus in the middle of the film [3], a choice that may not be appropriate.\nHere, we briefly discuss a potential problem with the seeding method. First, the interfacial energy cost to form the initial ice nucleus is not associated with standard ice structures, such as cubic or hexagonal ices, but rather with ice 0-like structures at the interface, as discussed in our previous paper [4] (we note that the existence and metastability of the ice 0 crystal was confirmed through density functional theory and quantum Monte Carlo simulations in Ref. [5]). When introducing a seed into a liquid, the surface structures of the seed should rapidly equilibrate, spontaneously forming ice 0-like structures on its surface if these structures are relevant. Consequently, the nucleation frequency in bulk should be appropriately estimated using the seeding method for bulk in this context.\nHowever, for thin films, it is crucial to recognise that identifying the likely locations of the formation of ice nuclei and understanding their shapes represents a distinct and more complex aspect of the investigation. For example, our study indicates that ice nuclei tend to form in the subsurface regions and not in the middle region. These considerations indicate that the selection of the seed's position and shape may play a crucial role in estimating the nucleation rate in the presence of symmetry-breaking surface fields, unlike in bulk.\n[1] Emily B. Moore and Valeria Molinero, Structural transformation in supercooled water controls the crystallization rate of ice. Nature 479, 506 (2011).\n[2] John A. Hayton, Michael B. Davies, Thomas F. Whale, Angelos Michaelides, and Stephen J. Cox,The limit of macroscopic homogeneous ice nucleation at the nanoscale, Faraday Discuss. (in press).\n[3] Amir Haji-Akbari and Pablo G. Debenedetti, Perspective: Surface freezing in water: A nexus of experiments and simulations, J. Chem. Phys. 147, 060901 (2017).\n[4] John Russo, Flavio Romano & Hajime Tanaka, New metastable form of ice and its role in the homogeneous crystallization of water, Nat. Mater. 13, 733–739 (2014).\n[5] D. Quigley, D. Alfè, and B. Slater, On the stability of ice 0, ice i, and Ih, J. Chem. Phys. Rapid Communications 141, 161102 (2014).\nConcerning the finite size effect, we do agree with the reviewer's suggestion that we should examine the finite size effect. We have indeed investigated this aspect by changing the film thickness from 4 nm to 8 nm and can confidently affirm that our conclusions remain unaffected by the finite-size effect. We would like to note using a film thicker than 8 nm becomes computationally impractical due to the substantial computational demands associated with the need for a large number of trajectories. Regardless of the film thickness or the system's geometry (be it flat films or nanodroplets), the consistent outcome is the nucleation of ice crystals at a distance of approximately 2 nanometers from the free surface. It is worth noting that our results, obtained using the mW water model that employs a short-range potential, are not so strongly influenced by finite-size effects. This is in contrast to systems with long-range electrostatic interactions, where such effects can be significant.\nWe believe that the provided explanations effectively address any concerns the reviewer may have had regarding the reliability of our conclusion.\nMore relevant questions should be addressed to claim surface induced ice nucleation: Did water film yield a higher nucleation rate than bulk water of the same volume, or did water film crystallize at a statistically higher temperature? Only a positive answer to either question can constitute a more meaningful evidence of surface crystallization. The authors have already obtained 50 independent trajectories which should allow them to extract unbiased ice nucleation rate (along with its statistical uncertainty) for water film. How does it compare with homogeneous ice nucleation at the same temperature (206 K)? If the rates are indistinguishable, how would one conclude free surface promotes ice nucleation?\nWe would like to express our appreciation to the reviewer for providing valuable comments.\nPrevious computational investigations in this field have largely approached the problem by examining the locations where ice crystals tend to form. However, we acknowledge the merit in the reviewer's suggestion to incorporate a direct comparison of nucleation frequencies between the thin water film and the bulk phase, as it can offer valuable insights into the phenomenon.\nWe conducted direct measurements of the nucleation rate per unit volume for both water films and bulk (refer to Fig. R1 above). The results indicate that at the same temperature, specifically for T=200 K and 206 K, the nucleation rate of the water film surpasses that of the bulk. This strong support for the promotion of ice nucleation by the free surface is evident.\nAlternatively, they can also investigate the free energy barrier. The authors made a strong claim that the order parameters used in previous studies did not capture the relevant pathway and argued the most relevant order parameter should be Q12. Although this reviewer does not necessarily agree on this, the authors can at least try to make it self-consistent by showing ice nucleation in subsurface indeed leads to a lower free energy barrier using the order parameter based on Q12. They do have all the data to allow them to do so.\nWe thank the reviewer for the valuable comment. Indeed, estimating the free energy barrier for crystal nucleation is an excellent approach to further validate our argument. However, as previously mentioned, we have already conducted the method suggested by the reviewer in the earlier comment and found that it supports the validity of our argument. Thus, we would like to leave this for future investigation.\nRegarding the use of Q12 order parameters, we have conducted extensive research that demonstrates their effectiveness in detecting various crystal structures, including ice 0. We previously made the discovery of a new metastable ice, named “ice 0”, and documented it in reference [4]. Subsequently, its existence and stability were confirmed through density functional theory and quantum Monte Carlo simulations in reference [5]. In our paper [4], we unveiled that ice-0-like structures serve as precursors or catalysts for the nucleation of cubic ice in bulk. In our current study, our findings demonstrate that these ice-0-like structures are preferentially generated near the free surface and play a pivotal role in promoting ice nucleation in proximity to that surface.\nIn reference [4], we also emphasised the critical importance of utilising the Q12 bond-orientational order (BOO) parameter when studying tetrahedral liquids. It is worth noting that the conventional Q6 parameter is unable to detect certain significant crystal structures, including metastable ice 0. For example, by employing the Q12 parameter, we not only uncovered the existence of the ice 0 crystal in water but also made a discovery of a new crystalline form in silicon. This form was later identified as the β-Sn crystal structure through a straightforward atom displacement and a reduction in cell symmetry, as described in reference [7]. These findings underscore the paramount importance of incorporating the Q12 order parameter into investigations of crystallisation in tetrahedral liquids.\nFurthermore, it is essential to highlight that using the Q6 parameter in biased Monte Carlo or forward-flux sampling simulations can potentially lead to the oversight of ice nucleation processes involving ice 0-like structures. We believe that this oversight may have contributed to the contrasting conclusions reached in previous computer simulations compared to our own research.\n[6] Flavio Romano, John Russo, and Hajime Tanaka, Novel stable crystalline phase for the Stillinger-Weber potential, Phys. Rev. B 90, 014204 (2014).\n[7] Domagoj Fijan, Mark Wilson, The characterisation of the “X” crystal structure in the Stillinger-Weber potential, Chem. Phys. Lett. 685, 316–321 (2017).\nOther than these major issues, the overall writing of the paper can be improved. Certain part of the paper is difficult to understand. For example, page 9-10 where the authors discussed the role of Laplace pressure and commented on previous studies. For those who may not know exactly these references, it’s rather difficult to understand the main points through these sentences because they are mixed with authors’ opinion. By the way, when claiming others’ conclusion “is not correct” multiple times in the paper, in what position do the authors believe they have to make such a strong statement? Lastly, when trying to generalize their conclusion to other materials like tetrahedral liquids, the authors should be aware of early studies that already attempted to tackle this problem, for example, Nat. Mater. 8, 726 (2009); JCP 18, 4102 (2016).\nWe thank the reviewer for his/her kind and valuable comments. Following suggestions, we have reorganised and improved the sections. Furthermore, we have incorporated the references mentioned by the reviewer into our revised manuscript to enhance the credibility and completeness of our work. We trust that our revised manuscript represents a substantial improvement, with our conclusions now being more robustly supported. We sincerely hope that the reviewer finds our revised work well-suited for publication in Nature Communications.\nReplies to the comments of Reviewer #2:\nThis study reports molecular simulations showing that ice tends to nucleate with greater probability at a certain distance from the ‘free surface’ of supercooled water, and attributes this to surface induced negative pressure favouring so-called Ice-0 like symmetry. It’s an interesting study however I have a few major concerns and would not suggest publication of the manuscript in its current form. I have detailed these concerns below.\nFirst, we thank the reviewer for carefully reading our manuscript and providing comments. We are pleased to learn that the reviewer found our study interesting. The feedback and guidance offered by the reviewer have been instrumental in elevating the quality and credibility of our research.\nBefore going to the specific details, we would like to emphasise that our simulations represent brute-force investigations of direct ice nucleation processes. Our findings consistently demonstrate that ice formations always occur in proximity to the free surface, regardless of film thickness or geometry, encompassing thin films and nanodroplets. It is important to note that we do not employ any enhanced sampling techniques or seeding methods in our simulations. As a result, our conclusion that ice formations tend to occur near the surface rather than in the central region is based on empirical observations and does not depend on any underlying assumptions, including the presence of ice 0-like structures. We believe that this strengthens the validity and reliability of our findings.\nIn the following, we provide individual responses to each of the reviewer's comments.\nMajor concerns:\nA central claim of the paper is that ice nucleation is favoured some distance, roughly 2 nm, from the water surface. I have a few issues with this. Firstly, it is not obvious to me that what is shown in Fig. 1 (b) is necessarily statistically significant. I would suggest the authors find some way of comparing their results to the null hypothesis that nucleation takes place uniformly across the box of water. The asymmetry of the peaks at roughly 2 and -2 nm, suggests that the finding may not be robust.\nWe appreciate the reviewer's valuable, constructive comments. We fully agree with the reviewer's observation that the double peak shown in Fig. 1(b) should be symmetric if our proposed mechanism is correct. However, we would like to clarify that the asymmetry in the distribution is primarily attributed to limited sampling. To address this issue, we have significantly increased the number of samples, resulting in a more symmetric distribution. We believe this improvement will remove the reviewer's concerns regarding the validity of our conclusion.\nFurthermore, we have undertaken a comparison of ice nucleation frequencies between thin films and the bulk, as suggested by Reviewer #1. The results are presented in Fig. R1 above in this response letter (new Fig. 1c in the manuscript). Our findings unequivocally confirm that the ice nucleation frequency is indeed higher in thin films than in the bulk.\nWe believe that these results will effectively address any concerns raised by the reviewer and substantiate the validity of our conclusion.\nI find the idea of negative pressure in a droplet or film a bit slippery. Do you not need to consider the direction of pressure? This said, wouldn’t you expect the negative pressure due to surface tension to be at a maximum closer to the surface than 2 nm? See Sega et al. (J. Chem. Theory Comput. 2016, 12, 9, 4509–4515). Similarly, if the claim is that enhancement in nucleation rate is due to negative pressure wouldn’t it make sense to work out a pressure profile for the system and compare it to the cluster distribution profile across the system?\nWe appreciate the thoughtful comment from the reviewer. Firstly, we need to consider the direction of pressure. In Fig. 2 in the manuscript, it is evident that the lateral direction displays significant negative pressure instead of the normal direction. This observation aligns with the findings in the work of Amir and Pablo [1]. Additionally, we consistently observe different profiles along the radial direction for the normal pressure and tangential pressure of a water droplet in Fig. 3e and Fig. 5 in Shahraza et al.’s study [2]. Specifically, it is the tangential direction that manifests negative pressure in a water droplet, with the minimum of the negative pressure occurring approximately 0.4~0.5nm from the surface.\nAccording to our results and Sega’s findings [3], the normal pressure P remains nearly constant N throughout the system for both flat and spherical surfaces, ensuring mechanical stability. The difference between the normal (P ) and lateral (P ) components of the pressure tensor is crucial for N T interfaces, as it characterises the surface tension: γ=P -P . It is the negative lateral and tangential N T pressure P that contributes to a high surface tension in proximity to the surface. our results are in T line with Sega’s observations. Moreover, it is the negative pressure in lateral (for flat surface) and tangential (for droplet) directions that plays an important role in lowering the free energy barrier for nucleation. In other words, nucleation is facilitated in regions with negative pressure, and the minimum of the negative pressure aligns well with the peak position of the Ice 0 density profile.\nHowever, pinpointing small nuclei precisely at the position corresponding to the minimum of the negative pressure is not necessary due to the substantial dynamic fluctuations of interfaces (see Fig.1a in the manuscript), i.e., capillary waves. As the ice nuclei grow to a sufficiently large size, we can observe a pronounced preferential position for stable ice nuclei in thermodynamics (see new Fig. 1b). However, since crystal nuclei formed near the surfaces can only grow inwards, the peak positions of stable ice nuclei (Nc > 60 molecules) are inevitably observed to be shifted inwards compared to the minimum position of the negative pressure.\nFinally, in response to the reviewer's suggestion, we have conducted an in-depth examination of the correlation between the pressure profile and the profiles of pre-ordered structures, i.e., the Ice 0-like clusters. Our findings reveal a close correlation between the two profiles, as illustrated in the newly added Figure 3e in the revised manuscript. We believe that this additional analysis will help to convince the reviewer of the validity of our argument. In fact, Shahraza et al. have assessed the pressure profile of water nanodroplets and clearly show a negative tangential pressure region close to the free surface [3], which strongly supports our findings.\n[1] Amir Haji-Akbari and Pablo G. Debenedettia, Proc. Natl. Acad. Sci. 114, 3316–3321 (2017).\n[2] Shahrazad M A Malek, Francesco Sciortino, Peter H Poole and Ivan Saika-Voivod, J. Phys.: Condens. Matter 30, 144005 (2018).\n[3] Marcello Sega, Balazs Fabian and Pal Jedlovszky, J. Chem. Theory Comput. 12, 4509–4515 (2016).\nI have a concern about the experimental testability and real-world relevance of the study. The introduction motivates the work of the basis of relevance to various applied problems where ice formation matters yet there is no discussion of how these findings might help us understand real systems. I think the question in the end is ‘do we expect a small water droplet, with a water-vapour interface to freeze at a warmer or colder temperature than would be anticipated for an equivalent volume of ‘bulk’ water?’. The work here would presumably predict that smaller droplets, containing a greater proportion of ‘interfacial’ water would tend to have a higher nucleation rate per volume than larger droplets? To the best of my knowledge, there is no experimental evidence for surface enhancement of apparently homogeneous ice nucleation (Atkinson et al. J. Phys. Chem. A 2016, 120, 33, 6513–6520 and Duft and Leisner Atmos. Chem. Phys., 2004, 4, 1997–2000). I would note here that some of the work cited (Shaw et al. J. Phys. Chem. B 2005109, 9865–9868) where nucleation is enhanced near a surface refers to heterogeneous ice nucleation and would seem to me to be of limited relevance to the present study.\nWe thank the reviewer for providing valuable comments and for sharing information regarding related papers. As discussed earlier, we conducted a comparison of ice nucleation frequencies between thin water films and bulk water in Fig. R1 above, and our findings indicate that the former exhibits a higher nucleation rate than the latter.\nRegarding the experimental papers mentioned by the reviewer, the paper by Atkinson et al. studied water droplets immersed in an oil. Thus, it is not relevant to our research.\nOn the other hand, the paper by Duft and Leisner compared the freezing probability for droplets of radius 49 μm and 19 μm. They concluded that homogeneous freezing is a volume-proportional process and that surface nucleation might only be important, if at all, for much smaller droplets. As recognized by these authors, the droplet size in this study exceeds ten microns. In such relatively large droplets, the volume of the thin surface region is significantly smaller compared to the overall droplet volume. Consequently, when estimating ice nucleation frequencies per unit volume, it becomes challenging to observe the pronounced impact of free-surface-induced ice nucleation.\nNonetheless, it is of paramount importance to underscore that the discovery regarding the facilitation of ice nucleation by the free surface holds profound significance within the realm of fundamental scientific understanding. In a previous study by Duft and Leisner, it was mentioned in the discussion section that, “It cannot be ruled out, however, that surface nucleation is important for much smaller droplets below 1 μm.”\nWe hope that our findings will stimulate and encourage further experimental investigations in this direction. It is worth noting, however, that for very small droplets, we must consider the influence of Laplace pressure. As is known, the surface curvature of a water droplet creates an additional pressure called the Laplace pressure, P =P (normal pressure)-P (tangential pressure), which L N T contributes to an increase in internal pressure. As the droplet size decreases, the internal pressure increases (P ), increases, potentially suppressing the probability of nucleation events in the droplet, L as shown in Li et al.'s work [4]. We have included this discussion into the revised manuscript.\nIn fact, by using the surface tension measurements in experiments, X.Z. Wu et al. [5] discovered surface freezing in liquid normal alkanes 30 years ago. They reported that the freezing temperature of liquid alkanes’ surface is up to 3°C above the bulk freezing point. Although water’s surface is different from that of alkanes, this experimental work on the surface freezing of liquid alkanes strongly suggests that conclusions regarding the promotion or suppression of surface freezing in liquids should be approached with careful consideration, rather than relying solely on a few predictions from enhanced sampling methods. Through thorough and direct brute-force simulations of the ice nucleation processes, we robustly provide computational evidence of surface freezing in water. This discovery not only resolves the long-standing controversy in computational studies on the effects of free surfaces on ice nucleation but also holds significant value for further exploration in experiments, emphasising the importance of our work.\n[4] Tianshu Li, Davide Donadio and Giulia Galli, Ice nucleation at the nanoscale probes no man’s land of water. Nat. Commun. 4, 1887 (2013).\n[5] X.Z. Wu, B.M. Ocko, E.B. Sirota, S.K. Sinha, M. Deutsch, B.H.Cao and M.W. Kim, Surface Tension Measurements of Surface Freezing in Liquid Normal Alkanes. Science 261, 1018 (1993).\nIt is worth noting that Hayton et al. (DOI: 10.1039/D3FD00099K (Paper) Faraday Discuss., 2023) have very recently attempted to assess how we might expect the thinness of water films to impact their ice nucleation rate compared to a more naively calculated homogeneous ice nucleation rate and found that the difference is minimal. An approach of this type might help clarify the relevance of the present study.\nWe are grateful to the reviewer for bringing to our attention the noteworthy paper by Hayton et al., which we were not aware of, as it has been accepted in Faraday Discussions but has not yet been published. In this paper, the authors stated in the Introduction, “The broad consensus from multiple simulation studies is that ice formation occurs away from the interface, in regions of the fluid that are bulk-like.” This consensus that has emerged within the scientific community stands in stark contrast to the findings we have uncovered in our study.\nAs emphasised at the beginning of our response, our simulations represent brute-force investigations of direct ice nucleation processes. Our findings consistently demonstrate that ice formations consistently occur in proximity to the free surface, regardless of film thickness or geometry, encompassing thin films and nanodroplets. It is important to note that we do not employ any enhanced sampling techniques or seeding methods in our simulations. As a result, our conclusion that ice formations tend to occur near the surface rather than in the central region is based on empirical observations and does not depend on any underlying assumptions, including the presence of ice 0-like structures. This strengthens the validity and reliability of our findings.\nThe use of enhanced sampling techniques is warranted only when employing appropriate order parameters for biasing, as we have emphasised in our manuscript. On the other hand, seeding techniques may prove valuable for examining crystal nucleation behaviours. However, it is essential to recognize that ice nucleation near a free surface follows a distinct kinetic pathway involving ice 0-like precursors. Moreover, the selection of the seed's position and shape may also play a crucial role in estimating the nucleation rate in the presence of symmetry-breaking surface fields, unlike in bulk. Consequently, the utility of seeding techniques to investigate the influence of the free surface on ice nucleation is questionable.\nTo illustrate this point, let's consider the interfacial energy cost in the initial ice nuclei. This cost is not associated with standard ice structures, such as cubic or hexagonal ices, but rather with ice 0-like structures at the interface, as extensively discussed in our previous paper [6]. Additionally, we note that the existence and stability of the ice 0 crystal were further validated through density functional theory and quantum Monte Carlo simulations in reference [7].\nBased on these findings, it is important to acknowledge that when introducing a seed into a liquid, the surface structures of the seed should rapidly equilibrate, leading to the spontaneous formation of an ice-0-like interface if it is relevant. Consequently, the nucleation frequency in bulk can be appropriately estimated using the seeding method for bulk in this context. However, it is crucial to recognise that identifying the likely locations of ice nuclei and understanding their shapes represents a distinct and more complex aspect of the investigation. This helps to elucidate why the seeding technique yields a conclusion that is contrary to our own findings.\n[6] John Russo, Flavio Romano & Hajime Tanaka, New metastable form of ice and its role in the homogeneous crystallization of water, Nat. Mater. 13, 733–739 (2014).\n[7] D. Quigley, D. Alfè, and B. Slater, On the stability of ice 0, ice i, and Ih, J. Chem. Phys. Rapid Communications 141, 161102 (2014).\nFurthermore, it is worth emphasising that our observation of ice nucleation in regions enriched with five-membered rings remains unexplainable if we were to assume that the initial ice nuclei are standard ice structures. This is primarily due to the intrinsic incompatibility of standard ice structures with the presence of five-membered rings.\nWe hope that the above explanations remove the reviewer’s concerns.\nFinally, I do not think that the importance of identifying ice-0 like precursors is sufficiently well articulated at present. I am all for fundamental science, but I think it important to lay out how knowing the structure of critical cluster in a given situation might have wider scientific impact.\nThe reviewer’s first comment seems to focus on our observation of “ice 0”. However, it is crucial to understand that metastable “ice 0” only exists transiently in the early stage of crystal nucleation. Ice 0-like structures formed near the free surface quickly transform into standard ice crystals, such as cubic ices. More importantly, as emphasized earlier, our conclusion, derived from brute-force simulations, that the free surface assists ice nucleation does not rely on the presence of ice 0. Thus, our observation primarily pertains to ice nucleation near the surface.\nWe fully agree with the reviewer’s opinion that it is critical how our finding has a broad scientific impact. In the paper by Duft and Leisner, it was mentioned in the discussion section that, “It cannot be ruled out, however, that surface nucleation is important for much smaller droplets below 1 μm.”\nWe hope that our findings will stimulate and encourage further experimental investigations in this direction. It is worth noting, however, that for very small droplets, we must consider the influence of Laplace pressure.\nWe believe that our revised manuscript represents a substantial improvement, with our conclusions now being more robustly supported. This discovery also resolves the long-standing controversy in computational studies on the effects of free surfaces on ice nucleation. We also trust that our findings will have a broad impact on various fields related to ice nucleation. We sincerely hope that the reviewer finds our revised work well-suited for publication in Nature Communications.\nReplies to the comments of Reviewer #3:\nThis manuscript investigates ice nucleation by means of molecular dynamics simulations. From a methodological standpoint, the authors study ice nucleation in water films and droplets, make use of two water models, namely mW and TIP4P/ice, and also put the spotlight on the relevance of chosing suitable order parameters.\nThe main result of the work consists in revealing a novel nucleation mechanism in which the free surface-induced negative pressure promotes nucleation via the emergence of Ice 0-like precursors with ring topologies containing 5 and 6-membered rings. Since these precursors are expected to be more friendly to liquid water as compared to standard (bulk) ice-like topologies, they are expected to lower the interfacial energy cost for ice nucleation. This finding would rectify previous statements regarding free surface order as “bad” for ice nucleation. Additionally, unlike substrate-induced ice nucleation that occurs directly at the outermost layer, here the process occurs slightly inwards (at the sub-surface layer, but very close to the surface). Hence, the authors not only reproduce the surface-induced crystallisation observed in experiments, but they also unravel its physical molecular underpinnings.\nFrom the above-expounded considerations I find this work to tackle a problem of utmost relevance, both from fundamental and practical perspectives, and to present a solid, comprehensive and carefully described approach that yields highly significant results. More importantly, the study provides a new conceptualization that proposes a solution to the core physical and molecular basis of the problem. Thus, I am glad to recommend publication.\nWe thank the reviewer for carefully reading our manuscript and valuable comments. We are delighted to receive such a positive evaluation of our manuscript and are pleased to learn that the reviewer has recommended publication with some revisions.\nThere is only one point I would like the authors to comment. The local orientational bond order, q_lm(i), is calculated by considering the four closest neighbours of a molecule. For the molecules at the first (outermost) layer, where important 3-fold hydrogen bond coordination is expected (with the presence of dangling bonds), the fourth neighbor would lie at the second coordination shell instead of belonging to the first-shell. Thus, unlike the quasi-2D lambda index or the study of five or six-membered rings, this index might produce artifacts in such case. Please comment on the application of this metric in your study.\nWe appreciate the valuable comment from the reviewer. Due to the layering phenomenon in the first two layers, the consideration of neighbours or bond coordination becomes crucial. The values of q_lm(i) with the four closest neighbours for the first layer may deviate from the standard. Accordingly, we have introduced new criteria for q_lm (i) for ice (Ic/Ih) to identify ice nuclei. Additionally, to avoid the influence of the number of neighbours, we applied several different order parameters, including the CHILL algorithm, quasi-2D lambda, and ring analysis, to fortify the robustness of our conclusion.\nTo address the reviewer’s concern, we have incorporated considerations for the specificity of the free surface in our calculations of both the local bond-order parameters and their coarse-grained counterparts. Comprehensive details about these modifications have been included in the Methods section. We believe that this addition will offer compelling evidence of our meticulous analysis of liquid structures near the surface.\nWe sincerely hope that the reviewer finds our revised work well-suited for publication in Nature Communications."}, {"round": 2, "reviewer_comments": {"reviewer_1": "In this revision, the authors attempted to address all the raised questions. The most significant\naddition to the original manuscript is the calculated nucleation rates for water films based on the\nbrute-force MD trajectories of spontaneous crystallization. The obtained ice nucleation rates for the\nfilms were found to increase as film gets thinner and they were shown to be greater than\nhomogeneous ice nucleation rate (but only by a margin). This new analysis certainly supports the\nconclusion. However there are two major issues that need to be addressed.\nThe first issue is an easy one: what is the statistical uncertainty for the estimated rate constant. Is\nthe increase in the nucleation rate statistically meaningful?\nThe second issue is more serious: I find it rather difficult to reconcile the discrepancy between\nthe current work and prior studies using the same force field. All the previous studies, e.g., Ref. 35,\nRef. 14, and Ref. 26, showed that water films led to a decrease of nucleation rate, which is the\nopposite to what’s found in this work. The authors attributed the difference to the application of\nadvanced sampling methods and the choice of order parameter in the prior studies. Although this\nmay be a possibility, the current study along with its analysis isn’t convincing enough to prove the\npoint. First, the applicability of advanced sampling approaches in homogeneous ice nucleation\nbased on mW has been cross validated against each other and also against brute-force simulation\nat low temperature. Second, the current work employs brute-force simulation, which is supposed\nto be immune to any potential bias in advanced sampling method, but the enhancement due to\nfree surface appears to be very subtle at low temperature, less than a factor of two between bulk\nand a 4 nm film. How relevant is this change in a more meaningful (much higher) temperature? It\nis likely this difference will become more significant as temperature increases, but this only stays\nas a speculation without an explicit demonstration that water film yields a sizable increase\n(orders of magnitude) in nucleation rate. Third, the authors also blamed on the q_6-based order\nparameter for not capturing the most relevant nucleation event. Again, if a Q_12 order parameter\nis used in any of these advanced sampling methods, will it yield a higher nucleation rate in film\ngeometry? I think the explicit answer to this question is required to justify a strong claim that\noverrules the prior multiple studies from multiple groups. Lastly, the enhancement of ice\nnucleation was attributed to the fraction of ice-0 increasing with negative pressure. I find this\nargument to be redundant. Negative pressure has already been shown to enhance homogeneous\nice nucleation by prior studies, without the necessity to incur ice-0, see for example, Rosky et. al.\nChem Phys Lett 789, 139289 (2022). In that study it was also shown that forward flux sampling\ncorrectly captures homogeneous ice nucleation rate and the rate enhancement with the negative\npressure. Therefore it is rather counterintuitive that the same method would fail to describe rate\nenhancement in water thin film if there is a such thing and the enhancement is indeed due to\nnegative pressure. I think if the authors make a strong claim against prior findings, all these\nquestions ought to be addressed explicitly.", "reviewer_3": "I am satisfied with the author's response and, thus, recommend publication."}, "author_response": "Responses to the comments of Reviewer #1\nIn this revision, the authors attempted to address all the raised questions. The most significant addition to the original manuscript is the calculated nucleation rates for water films based on the brute-force MD trajectories of spontaneous crystallization. The obtained ice nucleation rates for the films were found to increase as film gets thinner and they were shown to be greater than homogeneous ice nucleation rate (but only by a margin). This new analysis certainly supports the conclusion. However there are two major issues that need to be addressed.\nWe would like to express our sincere appreciation to the reviewer for their invaluable feedback on our revised manuscript. Their constructive comments and suggestions from the previous review round have greatly improved the quality of our work. In response to the two issues raised, we will address them individually as follows.\nThe first issue is an easy one: what is the statistical uncertainty for the estimated rate constant. Is the increase in the nucleation rate statistically meaningful?\nIn this work, we directly measure the nucleation rate by using the distribution of the induction times of the nucleation,\nwhere N is the total number of simulations performed for each system, t i is the induction time determined for the i-th simulation, and Θ(x) is the Heavisi de step function.\nThe number of nucleation trajectories does indeed have a slight influence on the shape of P (t). However, the distribution of induction times converges when the number of independent statistical simulation trajectories exceeds 50. In our revised manuscript, we conducted 100 independent simulations for each system under identical conditions. Therefore, the statistical uncertainty regarding the estimated rate constant does not impact the conclusions drawn from our investigation into the nucleation behavior of water films.\nTo confirm the statistical uncertainty, we present the nucleation rates of water films with L=6 nm and 4 nm at T=200 K and 206 K, respectively, as a function of the number of the statistical simulation trajectories (50~200 trajectory samplings).\nFigure R1. The nucleation rate per volume, R/V, as a function of the number of statistical simulations. This figure includes two water films with different thicknesses, L=6 nm and 4 nm, at T=206 K and 200 K, respectively.\nWe observe that the statistical uncertainty for the estimated nucleation rate, determined through the distribution of induction times of nucleation, is quite small and can be considered negligible. Hence, the nucleation rate of water films obtained through direct measurement holds statistical significance.\nThe second issue is more serious: I find it rather difficult to reconcile the discrepancy between the current work and prior studies using the same force field. All the previous studies, e.g., Ref. 35, Ref. 14, and Ref. 26, showed that water films led to a decrease of nucleation rate, which is the opposite to what’s found in this work. The authors attributed the difference to the application of advanced sampling methods and the choice of order parameter in the prior studies. Although this may be a possibility, the current study along with its analysis isn’t convincing enough to prove the point.\nFirst, the applicability of advanced sampling approaches in homogeneous ice nucleation based on mW has been cross validated against each other and also against brute-force simulation at low temperature. Second, the current work employs brute-force simulation, which is supposed to be immune to any potential bias in advanced sampling method, but the enhancement due to free surface appears to be very subtle at low temperature, less than a factor of two between bulk and a 4 nm film. How relevant is this change in a more meaningful (much higher) temperature? It is likely this difference will become more significant as temperature increases, but this only stays as a speculation without an explicit demonstration that water film yields a sizable increase (orders of magnitude) in nucleation rate. Third, the authors also blamed on the q_6-based order parameter for not capturing the most relevant nucleation event. Again, if a Q_12 order parameter is used in any of these advanced sampling methods, will it yield a higher nucleation rate in film geometry? I think the explicit answer to this question is required to justify a strong claim that overrules the prior multiple studies from multiple groups. Lastly, the enhancement of ice nucleation was attributed to the fraction of ice-0 increasing with negative pressure. I find this argument to be redundant. Negative pressure has already been shown to enhance homogeneous ice nucleation by prior studies, without the necessity to incur ice-0, see for example, Rosky et. al. Chem Phys Lett 789, 139289 (2022). In that study it was also shown that forward flux sampling correctly captures homogeneous ice nucleation rate and the rate enhancement with the negative pressure. Therefore it is rather counterintuitive that the same method would fail to describe rate enhancement in water thin film if there is a such thing and the enhancement is indeed due to negative pressure. I think if the authors make a strong claim against prior findings, all these questions ought to be addressed explicitly.\nWe appreciate the reviewer's raising this issue. First, we want to clarify that our work does not challenge the validity of advanced sampling approaches in ice nucleation studies. Rather, our point in the manuscript is that these advanced sampling methods, which place intermediate states along the reaction coordinate or an order parameter, may be sensitive to the choice of order parameters.\nThe local bond orientational order parameter q , commonly used in previous studies (e.g., Ref. 6 14), fails to adequately describe water nucleation in heterogeneous systems. This could be due to the following reasons:\n(1) The local 6th-order harmonic q is not ideal for capturing the correct symmetry of crystals of 6 tetrahedral liquids, including the ice lattice [1,2]. Instead, the 3rd-order (q ) [1,2] and the 12th- 3 order harmonic (Q ) [3], coarse-grained up to neighboring molecules, are more suitable choices 12 for identifying crystal symmetry for tetrahedral liquids. This is particularly evident in the case of the coarse-grained 12th-order harmonic, Q (refer to Fig. R3c).\n12\n(2) The local order parameter q is measured using the nearest neighbors, which can pick up the 6 ordinary ice structures (hexagonal and cubic ices). In a low-temperature liquid, water molecules and their nearest neighbors tend to form tetrahedral structures to minimize energy. However, to form ice crystalline structures, these tetrahedral geometries need further adjustments to establish specific connections, which involve more neighbors, including at least the second nearest neighbors. Therefore, the local order parameter q generally tends to overidentify the ice 6 molecules in the liquid phase (refer to Fig. R3a).\nFigure R2. Illustration of local water geometry and ice structure. (a) Tetrahedral structure in water formed with its nearest neighbors. (b) The ice structure formed b y a water molecule and its 16 neighbors through a special tetrahedron connection.\nFigure R3. The ice molecules in a water film with L=8 nm at T=206 K and t=2.63 ns. (a) Ice molecules identified using local q as in Refs. 14, 26, and 35. (b) Ice molecules identified using 6 the local q [2]. (c) Ice molecules identified using the coarse-grained Q [3].\n3 12\nIn conclusion, we do not question the validity of advanced sampling methods in liquid films or the various potential models. However, it is important to note that these advanced sampling methods are sensitive to the choice of order parameter. Previous studies have demonstrated that the local 6th-order harmonic, q , is not the suitable choice for carbon [1] and silicon [2], as it 6 may overestimate the presence of crystalline molecules during nucleation events. On the other hand, the coarse-grained order parameter utilizing the 12th-order harmonic, Q , offers a clear 12 and efficient method for identifying ice molecules (refer to Fig. R3). This issue may not significantly affect the estimation of nucleation rates in bulk systems. However, in cases where surfaces induce ordering with symmetries other than 6-fold, as in our scenario, using q for 6 biasing could potentially lead to incorrect conclusions.\nWe have incorporated the above information into the revised manuscript and have also referenced the paper kindly notified by the reviewer (Rosky et al. Chem Phys Lett 789, 139289 (2022). These revisions are highlighted in blue. We have also included Fig. R3 as a part of Supplementary Fig. S1.\nWe believe that the above responses adequately address the reviewer's concerns and effectively convey the validity of our findings. We sincerely hope that the reviewer finds our revised manuscript suitable for publication in Nature Communications.\n[1] Ghiringhelli, L. M., Valeriani, C., Meijer, E. j. & Frenkel, D. Phys. Rev. Lett. 99, 055702 (2007). In this reference, the authors stated, “We used a local order parameter as introduced in Ref. [6]; rather than the 6th order harmonic, we found that a 3rd order harmonic can capture the correct symmetry for both graphite and diamond lattices.”\n[2] Li, T., Donadio, D., Ghiringhelli, L. M. & Galli, G. Nat. Mater. 8, 726-730 (2009). In this reference, the authors stated in the Methods section, “We use the local q (ref. 24, the above [1]), 3 a quantity that is highly sensitive to crystalline order, to identify Si crystalline clusters in the liquid.\n[3] Russo, J., Romano, F.& Tanaka, H. Nat. Mater. 13, 733-739 (2014).\nResponses to the comments of Reviewer #3\nI am satisfied with the author's response and, thus, recommend publication.\nWe are grateful to the reviewer for taking the time to review our revised manuscript again. We are delighted that the reviewer has recommended the publication of our paper in Nature Communications."}]

Decision letter
REVIEWERS' COMMENTS
Reviewer #1 (Remarks to the Author):
The authors have conducted additional brute-force simulations to confirm the validity of their
conclusion. They also provided more analysis and reasoning to support their finding and try to
reconcile the discrepancy between the current work and prior studies through q6 and Q12.
Although this still isn’t sufficient either to prove the relevance of their finding in a more relevant
temperature or to explain the contradiction between the current and prior studies, I would say the
current work has its own merit as it shows through unbiased simulation that film geometry could
lead to a higher ice nucleation rate and it also points out the potential issue of q6. I find Fig. 3R to
be interesting and I would suggest the authors to include this figure in the paper. Recent two-path
model (JCP 158, 124501 (2023)) could help rationalize the discrepancy provided that it’s proven
that the pathway near the surface identified through q6 indeed leads to a lower or comparable rate
than homogeneous nucleation.
The authors have conducted additional brute-force simulations to confirm the validity of
their conclusion. They also provided more analysis and reasoning to support their finding
and try to reconcile the discrepancy between the current work and prior studies through
q6 and Q12. Although this still isn’t sufficient either to prove the relevance of their
finding in a more relevant temperature or to explain the contradiction between the
current and prior studies, I would say the current work has its own merit as it shows
through unbiased simulation that film geometry could lead to a higher ice nucleation
rate and it also points out the potential issue of q6. I find Fig. 3R to be interesting and I
would suggest the authors to include this figure in the paper. Recent two-path
model (JCP 158, 124501 (2023)) could help rationalize the discrepancy provided that it’s
proven that the pathway near the surface identified through q6 indeed leads to a lower or
comparable rate than homogeneous nucleation.
We first thank the reviewer for carefully reviewing our manuscript again. We are
delighted to learn that the reviewer found our explanation reasonable. Following the
reviewer's kind advice, we have cited the paper suggested by the reviewer. We have also
included Fig. R3 in panels a-c in new Fig. 1.
We trust that our revised manuscript is now suitable for publication in Nature
Communication.