The limit of macroscopic homogeneous ice nucleation at the nanoscale †

Nucleation in small volumes of water has garnered renewed interest due to the relevance of pore condensation and freezing under conditions of low partial pressures of water, such as in the upper troposphere. Molecular simulations can in principle provide insight on this process at the molecular scale that is challenging to achieve experimentally. However, there are discrepancies in the literature as to whether the rate in con ﬁ ned systems is enhanced or suppressed relative to bulk water at the same temperature and pressure. In this study, we investigate the extent to which the size of the critical nucleus and the rate at which it grows in thin ﬁ lms of water are a ﬀ ected by the thickness of the ﬁ lm. Our results suggest that nucleation remains bulk-like in ﬁ lms that are barely large enough accommodate a critical nucleus. This conclusion seems robust to the presence of physical con ﬁ ning boundaries. We also discuss the di ﬃ culties in unambiguously determining homogeneous nucleation rates in nanoscale systems, owing to the challenges in de ﬁ ning the volume. Our results suggest any impact on a ﬁ lm ’ s thickness on the rate is largely inconsequential for present day experiments.


Introduction
4][15][16] For example, many properties of clouds are affected by the relative compositions of ice and water, 3,[17][18][19] and consequently, the accuracy of climate models rely heavily on parameterizations to predict ice nucleation in the atmosphere.
Broadly speaking, when ice forms, it can do so either heterogeneously, where the surface of, say, a solid particle facilitates nucleation, or homogeneously, in the absence of such surfaces.Despite heterogeneous nucleation being by far more common, homogeneous nucleation is important when temperatures approach approx.1][22][23][24] Yet, even in the absence of surfaces presented by solid particles, the finite volume occupied by the liquid means an interface (e.g., with vapor or oil) nonetheless remains.A longstanding issue for homogeneous nucleation has therefore been to establish whether nucleation is enhanced close to the liquid-vapor interface, or suppressed.Owing to the small length and fast times scales involved, however, it is experimentally challenging to establish whether homogeneous nucleation occurs near the interface, or in the bulk of the fluid.6][27][28][29][30] While pioneering simulation studies from Jungwirth 25,26 and co-workers suggested that nucleation is enhanced at the liquid-vapor interface, their results are likely affected by the finite size of the simulation cell. 31The broad consensus from multiple simulation studies is that ice formation occurs away from the interface, in regions of the fluid that are bulk-like.These conclusions are also supported by thermodynamic arguments made by Qiu and Molinero. 32 ice forms preferentially in the bulk of the liquid, what does this mean for the observed rate of ice nucleation?In 2013, using forward flux sampling, Li et al. 27 found, upon decreasing the droplet radius from approx.4.9 nm to approx.2.4 nm, that nucleation rates at 230 K were decreased relative to that of bulk water by eight orders of magnitude; for radii ≳ 5 nm, nucleation rates were virtually indistinguishable from that of bulk.With the exception of the smallest droplet investigated, this significant decrease of the rate upon decreasing the radius appeared to be explained well by the associated change in Laplace pressure.These findings were consistent with a previous study by Johnston and Molinero. 28However, subsequent simulation studies that investigated ice nucleation in thin water films, 29,33 whose planar interfaces correspond to zero Laplace pressure, found that nucleation rates were noticeably suppressed relative to that of bulk, even for film thicknesses ≳ 5 nm.Owing to the computationally demanding nature of simulating ice formation, these studies used a coarse-grained representation of water's interactions, the mW model. 34In this model, instead of explicitly representing the hydrogen bond network, the local tetrahedrality is enforced by a three body contribution to the potential energy function.][36][37] To add further complication to the above apparent discrepancy between droplets and films, tour-de-force simulations employing the TIP4P/ice model, 38 which unlike the coarse grained mW model, accounts explicitly for electrostatic interactions and water's hydrogen bond network, found that nucleation in thin water films was increased relative to that of bulk water by approx.seven orders of magnitude. 30This increase was attributed to an increase in average crystalline order and cage-like structure for TIP4P/ice (in the liquid state) in a thin-film geometry, which persisted to film thicknesses even as large as 9 nm.These pronounced structural changes were not observed with the coarse grained mW model.While the origin of such a discrepancy might be attributed to a lack of essential physics in the coarse grained model, this may raise a question mark over the use of mW to investigate ice formation in confined geometries, despite its success, 37,39,40 e.g., in explaining ice nucleation from vapor via a pore condensation and freezing mechanism. 41n this article, we first aim to resolve this apparent discrepancy between the coarse grained mW model and the all-atom TIP4P/ice model.Then, in Sec. 3, using the computationally more efficient mW model, we directly assess how the simplest confining geometry of all-a thin film of water in contact with its vaporimpacts ice formation.In Sec. 4 we extend our investigation to a more realistic case where water is confined between two solid surfaces.We conclude our findings in Sec. 5. A broad overview of the methodology is described throughout the manuscript, with full details provided in Sec.6 and the Electronic Supporting Information (ESI †).

Structural properties of thin water films converge to their bulk values on a microscopic length scale
As discussed above, previous work from Haji-Akbari and Debenedetti, 30 using the TIP4P/ice water model at 230 K has reported that the nucleation rate J(W ) in a film of thickness W ≈ 4 nm is roughly seven orders of magnitude larger than in bulk water, whose rate we denote J(∞).This significant enhancement of the rate has been attributed to a structure of liquid water "deep" in the interior of thin films that is distinctly different from that found for homogeneous bulk water.In particular, measures of local order, as prescribed, e.g., by a Steinhardt order parameter 42 for the i th water molecule, indicate that the average structure of liquid water is more icelike for thin films when compared to homogeneous bulk water under the same conditions.‡ The prime indicates that the sum only includes the ν (i) nearest neighbors of molecule i (see ESI †), where ri j is the unit vector pointing from molecule i to molecule j, and Y 6m is the m th component of a sixth-rank spherical harmonic.
Similarly, the number of cage-like structures was also found to be increased in thin water films.In contrast, when using the mW potential, such measures of local structure converge to their bulk values within approx. 1 nm of the interface.
In a recent study, Atherton et al. 46 discussed the sensitivity of the melting point of simple point charge models such as TIP4P/ice to the choice of cutoff, r * , for the non-electrostatic interactions between water molecules.There it was shown that decreasing r * led to a systematic increase in the melting temperature.More importantly, it was noticed that choices of r * typically used in molecular simulations effectively correspond to negative pressures of a few hundred bar, when compared to homogeneous bulk systems that use a mean-field treatment to correct for truncated interactions.This observation was argued to be particularly relevant when comparing nucleation in thin water films (where standard mean-field treatments of truncated interactions have no impact) to homogeneous bulk systems.In Ref. 46, rough arguments based on results from Bianco et al. 47 were used to suggest that this effective negative pressure was the root of faster ice nucleation in thin films, but firm evidence was lacking.
For a thin film of water in coexistence with its vapor, it can readily be argued that the pressure far from the interfaces is approximately 0 bar. 48Therefore, to obtain a suitable reference for the local structure in the homogeneous bulk fluid, we have performed simulations in the isothermal-isobaric ensemble at a temperature T = 230 K, and pressure p = 0 bar (see ESI † for full simulation details), and computed where ρ is the density of bulk water, N is the number of water molecules in the simulation, and angled brackets indicate an ensemble average.An important detail is that, when using the TIP4P/ice potential, we have truncated and shifted nonelectrostatic interactions at r * = 8.5 Å; following Ref.46 we indicate this molecular model as TIP4P/ice (8.5) .We find ⟨ q6 ⟩/ ρ ≈ 0.266.We now turn our attention to a thin film, as shown in Fig. 1a, also at T = 230 K.The specific system we consider comprises 3072 TIP4P/ice (8.5) water molecules, and the lateral dimensions of the simulation cell are chosen such that W ≈ 4.3 nm.Throughout this paper, the width of a film is defined as W = N/A ρ, where A is Fig. 1 Structural properties of TIP4P/ice (8.5) and mW converge to their bulk values within approx.1-1.5 nm of the liquid-vapor interface.(a) Snapshot of a TIP4P/ice (8.5) film of thickness W ≈ 4.3 nm at T = 230 K.The z axis (normal to the liquid-vapor interface) lies along the horizontal.Local structure away from the interface soon becomes bulk-like, as shown by ⟨ q6 (z)⟩/⟨ρ(z)⟩ in (b) for TIP4P/ice (8.5) (blue) and mW (orange).Dashed lines show ⟨ q6 ⟩/ ρ obtained from simulations of the bulk fluid at T = 230 K and p = 0 bar.When truncated interactions are accounted for in a mean-field fashion, a discrepancy in local structure in the center of the film and homogeneous bulk water emerges, as indicated by the dashed green line.The vertical dotted lines indicate the location of the liquid-vapor interface at z = ±W /2. the cross-sectional area.W is then varied by changing A. This simple measure is roughly consistent with the separation of Gibbs dividing surfaces.
To analyze the spatial variation in the structure of this film, we compute where z is the coordinate normal to the average liquid-vapor interface, z i is the z coordinate of i th water molecule's oxygen atom, and ⟨ρ(z)⟩ = 1 is the number density profile.The result of this analysis is shown in Fig. 1b; While some interesting structure is observed close to the interface on the liquid side, the important point is that ⟨ q6 (z)⟩/⟨ρ(z)⟩ ≈ 0.266 converges to its bulk value within approx.1-1.5 nm of the interface, as indicated by the dashed blue line in Fig. 1b.In addition, the dashed green line shows ⟨ q6 ⟩/⟨ρ(z)⟩ ≈ 0.260 obtained from a simulation of the bulk fluid at T = 230 K and p = 0 bar in which a mean-field treatment for truncated interactions has been employed; we denote this molecular model TIP4P/ice (8.5→∞) .We see that the result obtained with TIP4P/ice (8.5→∞) lies, on average, below ⟨ q6 (z)⟩/⟨ρ(z)⟩ obtained with TIP4P/ice (8.5) .Also shown in Fig. 1b is a similar analysis for the mW model.While we see quantitative differences with the TIP4P/ice (8.5) result, the important point is that ⟨ q6 ⟩/⟨ρ(z)⟩ ≈ 0.226 converges to its bulk value in the center of the film, on a similar length scale to TIP4P/ice (8.5) .This result demonstrates that the structural differences between mW and TIP4P/ice that lead to apparently qualitative differences in nucleation rates can be resolved through a consistent treatment of truncated interactions.Combined with previous observations that the mechanism of ice formation in thin films is similar between the two models, 27,29,30 the results in this section strongly suggest that the coarse-grained mW model likely contains the essential physics to describe ice nucleation in these systems.
In the following section, we will exploit the computational efficiency of the mW model to explore how the nucleation rate in thin water films depends upon W . Specifically, using the seeding technique, we will investigate how both the size of the critical nucleus and its growth depend upon W . Preliminary results for TIP4P/ice (8.5) with W = 4.3 nm are given in the ESI †. in describing the structure of thin water films, we will now investigate the extent to which W impacts ice nucleation.2][53][54][55][56][57] As the seeding approach has been detailed previously by others, 51,52,54,55 here we will only briefly sketch an outline of the procedure.

Investigating ice nucleation in thin water films with seeding simulations
The principal idea behind the seeding method is to initialize the system with a preformed ice nucleus (a 'seed') and observe whether the seed, on average, grows such that the system ends up as ice, or shrinks such that the system ends up as liquid water.
If the seed is smaller than the critical ice nucleus, it will tend to melt, whereas if it is larger, it will tend to grow; the size of the critical ice nucleus can be determined from seeds that have equal tendency to grow or melt.In our case, this was achieved by first equilibrating a bulk crystal of hexagonal ice comprising 16000 molecules at 220 K and 0 bar, and then carving out a spherical cluster.The cluster was then inserted into the center of a thin film of liquid water of thickness W , which itself had been equilibrated at 220 K; this was done after first removing water molecules to create a spherical cavity of an appropriate size to accommodate the ice cluster.With the molecules in the ice cluster held fixed, the surrounding fluid was then relaxed by performing a short 80 ps molecular dynamics simulation.We subsequently calculated the size of the largest cluster of ice-like molecules in the system, n cl , which we took to define the size of our initial ice seed.An example of such a seed is shown in Fig. 2.
We have investigated film thicknesses in the range 2.5 nm ≲ W ≲ 6 nm, with N ≈ 6000 throughout (i.e., different thicknesses are achieved by varying A).For each W , approximately 700-900 seeding simulations were performed.The range of initial cluster sizes depended on film width, with 60 ≲ n cl ≲ 220 overall; this was sufficient to span both pre-and post-critical cluster sizes.In Fig. 3a, we present the probability p ice (n cl ;W ) that a seed of size n cl goes on to form ice.The critical cluster size n ‡ cl is estimated from p ice (n ‡ cl ;W ) = 0.5.We observe that for W ≳ 3.5 nm, n ‡ cl ≈ 120, which compares well to the value obtained in a bulk simulation of mW.To help give this result some perspective, the radius of these critical clusters is r ‡ cl ≈ 0.96 nm.If we consider W = 3.5 nm, this gives (W −2r ‡ cl )/2 ≈ 0.8 nm, which is broadly in line with accepted thicknesses of the liquid-vapor interface. 27,29,58In other words, as soon as the films are thick enough to accommodate a bulk-like region, critical nuclei appear ambivalent to the nearby presence of the liquid-vapor interface.For W ≲ 3 nm, deviations of n ‡ cl from its bulk value are observed, with n ‡ cl ≈ 170 in the thinnest film investigated (W ≈ 2.5 nm).For all systems investigated, n ‡ cl vs. W is plotted in Fig. 3b.Table 1 Summary of results from seeding simulations.For each W , the range of critical cluster sizes (indicated in parentheses) is obtained by splitting its data set into three.The range for J(W ) is calculated by using the largest and smallest values for n ‡ cl for a given W in Eq. 6 (∆µ/k B = 122 K and ρ = 33.3774nm −3 ), and is also reported in parentheses Computation of nucleation rates with the seeding approach relies upon classical nucleation theory (CNT): [49][50][51][52] where ∆µ is the chemical potential difference between bulk ice and liquid is the Zeldovich factor, and f + is the attachment frequency.For the chemical potential difference, a range of values spanning approx.118 K ≲ |∆µ|/k B ≲ 126 K have been reported in the literature 59,60 and for simplicity, we take |∆µ|/k B = 122 K. (We have repeated the following analyses with the extremal values of this range, and find that the results are virtually indistinguishable [see ESI †].)As we have found that n ‡ cl remains roughly constant for W ≳ 3.5 nm, the remaining source for deviations of the rate from that of bulk may lie in the dynamics, as codified by f + .In Fig. 3c we present the attachment frequency obtained from as first proposed by Auer and Frenkel. 61,62In practice, , it is clear that the attachment frequency is largely insensitive to W . Bringing together the above, in Fig. 3d we present J(W ) for the different film thicknesses investigated (see also Table 1).As a check of our implementation, our result for J(∞)-obtained from a simulation of bulk water under periodic boundary conditions-agrees well with that previously reported by Sanchez-Burgos et al. 63 at T = 220 K and p = 1 bar.The biggest drawbacks of the seeding method are that it presupposes the nucleation mechanism and relies upon CNT to compute J. Previously, Lupi et al. 64 have provided strong evidence that using the size of the critical cluster is a good "reaction coordinate" for nucleation; in the ESI † we present committor analyses for W = 5 nm and W = 2.5 nm that suggests that this is also the case for the thin film systems we consider.

Comparing, and reconciling, results from seeding with previous studies
Despite good agreement of J(∞) with previous literature values, our results for W ∼ 5 nm are at odds with previous work.For example, using forward flux sampling, for mW under the same conditions, Haji-Akbari et al. 29 found J(5 nm)/J(∞) ≈ 0.006.Similarly, by computing mean first passage times at 205 K, Lü et al. 33 found that ice nucleation rates in thin water films were suppressed compared to that of bulk.How can we rationalize this apparent discrepancy of our seeding simulations, which find that J(W )/J(∞) ≈ 1 once the film is thick enough to accommodate a critical nucleus?We, in fact, argue that this discrepancy is largely superficial.In the case of Lü et al., the reported impact of W on the rate was relatively modest, with J(5 nm)/J(∞) ≈ 0.5, within the range of uncertainty of our result.But there is also a more subtle aspect at play.Methods such as forward flux sampling and calculating the mean first passage time obtain a rate first by computing the probability to undergo a nucleation event per unit time in a particular sample, and then normalize by the volume of the sample.For homogeneous systems, the definition of this volume is unambiguous.In contrast, for inhomogeneous systems such as thin water films, what to take for the volume is less clear-cut, and J(W ) becomes increasingly sensitive to this normalization procedure as W decreases.The seeding method, on the other hand, provides an estimate of the rate directly (see Eq. 6), albeit conditioned, in our study, on the nucleus forming in the bulk-like region.The advantage of the seeding method is that it clearly demonstrates that n ‡ cl and f + remain roughly constant for W ≥ 3.5 nm.We note that in the supporting information of Li et al., 27 who used forward flux sampling, the nucleation rate in a film comprising 4096 mW water molecules was found to be indistinguishable from that of bulk; it appears that Li et al. account for the surface region in their normalization procedure.
To help further understand the differences between J(W ) computed with seeding and previous works, it is instructive to introduce the following simple model, which is similar to existing models introduced by others. 27,33,65We suppose that the thin film can be separated into two interfacial regions of thickness ℓ i that sandwich a central bulk-like region of thickness W − 2ℓ i .Accounting for the volume occupied by the critical nucleus, the volume accessible for nucleation in the bulk-like region is A(W − 2ℓ i − 2r ‡ cl ) (see Fig. 4).In the remaining volume, 2A(ℓ i +r ‡ cl ), we assume that nucleation occurs with a rate J i .For a sample of total volume AW , the number of nuclei that form per unit time is determined by an effective nucleation rate J eff : In (b) we present the fraction of frozen samples for W /nm = 3.5 and 6, and bulk water.Each system comprises 6000 mW molecules and is cooled at a rate of 0.2 K/ns from an initial temperature of 220 K.There is a modest shift to lower temperatures with the films compared to bulk water.The freezing rate, R, for W = 3.5 nm, as shown in (c), differs at most by a factor five compared to the bulk result.
Rearranging for J eff , and assuming that J i ≈ 0, we find Using r ‡ cl ≈ 0.96 nm obtained from our seeding simulations, a suppression of the nucleation rate of J eff (5 nm)/J(∞) ≈ 0.006, as obtained from forward flux sampling, 29 requires ℓ i ≈ 1.53 nm.While this order of magnitude is reasonable for an interfacial thickness, it is inconsistent with our finding that the size of critical nucleus remains constant for W ≳ 3.5 nm.As we did above, we can instead estimate ℓ i = (3.5 nm − 2r ‡ cl )/2 ≈ 0.8 nm.For W = 5 nm, we then obtain J eff (5 nm)/J(∞) ≈ 0.3.This more modest suppression in the rate appears comparable to that obtained at 205 K by Lü et al. 33

Probing the effective rate in thin water films through the lens of an experimentalist
What the above analysis demonstrates is the sensitivity of the effective nucleation rate to the interfacial thickness ℓ i and the size of the critical nucleus r ‡ cl .While the latter is in principle a physical observable, there is a degree of arbitrariness in defining an interfacial thickness, which obfuscates direct interpretation of J eff .7][68][69] In such cases, one is probably less interested in the precise value of an effective rate, and instead more concerned whether ice can form, given that a pore/crack is filled with water under the conditions.
In this spirit, in Fig. 4b we present the fraction of frozen samples against temperature for bulk water, W = 6 nm and W = 3.5 nm, obtained by cooling 100 replicas of each system at 0.2 K ns −1 .Note that the number of water molecules is the same in each case (6000 mW).The size of the largest ice-like cluster was monitored, and the system was considered frozen when the largest cluster n cl (T ) ≥ 200.§ We see that the curves are shifted relative to each other, with the bulk samples freezing at the highest temperature, and the W = 3.5 nm film at the lowest temperature being separated by around 2 K.We are able to achieve statistically different freezing temperatures and rates in our studies primarily due the artificially high control we can exert on the number of molecules (N = 6000) and with it the volume undergoing cooling.
In experimental systems such fine control is impossible, and with it the deviations in freezing temperatures and freezing rates reported here become unimportant.To see this, we compare the impact of varying W by computing the freezing rate, R, for each of the data sets using an approach typically employed in the analysis of laboratory experiments. 70R is analogous to the radioactive decay constant (a first-order reaction rate constant) and is most easily determined from experiments conducted on completely iden-tical supercooled droplets under isothermal conditions.In the isothermal case, the fraction of droplets which freeze in a given time period can be found as: 71 N liq (t) where N liq (t) is the number of liquid droplets remaining after time t and N tot is the total number of droplets (liquid or solid) present in the experiment [equal to N liq (t = 0)].In order to calculate R from continuously cooled experiments, it is necessary to divide the experiment into small time intervals ∆t over which changes in temperature are treated as small.Rearranging Eqn. 10 we find where N i is the number of unfrozen droplets present at the start of the time interval and N f is the number of droplets which freeze during the time interval.We have calculated R for the three systems by treating each of their 100 replicas as a droplet.
We have used ∆t = 5 × 10 −9 s, resulting in temperature bins 1 K wide.The results of this analysis are shown in Fig. 4c, where the data points indicate the midpoint of the temperature bins.Confidence intervals were calculated using a simple Monte Carlo simulation implemented in Stata 17. 72 For physically identical droplets it is expected that the number of freezing events in a temperature interval will follow a Poisson distribution on repeat testing. 71,73For each temperature bin, the observed number of freezing events was taken as the expectation value for the bin.
Poisson distributed random numbers were then generated to create 5000 simulated repeats of each of the three data sets.The bars on Fig. 4c show the 5 th to 95 th percentile range of the simulated experiments.(Data points generated from bins containing so few freezing events that the 5 th percentile simulation had a value of R = 0 have not been included in Fig. 4c.)The freezing rates determined for the bulk system are clearly higher than those found for W = 3.5 nm.This difference is largest for the bins centered at 201.65 K, where the bulk rate is greater by approx.a factor of five.
In principle, the homogeneous nucleation rate, J, can be found as R/V where V is the volume of the water droplets undergoing freezing.Taking the volume of an mW water molecule in the liquid phase as 1/ ρ ≈ 3 × 10 −29 m 3 we find that 6000 water molecules occupy 1.8 × 10 −25 m 3 .Dividing our calculated freezing rates by this volume we obtain values of J entirely consistent with previous literature values 29,35,63,74 for mW (see ESI †).However, obtaining J in this fashion is made possible as we know the exact number of water molecules used in our simulations.Experimentally, it is infeasible to determine the volume with such precision, so the differences in rates we report are likely inconsequential in any practical setting.

Assessing the impact of physical boundaries
The thin water films we have investigated so far are able to provide insight into the fundamental question of how much liquid water is needed to recover bulk-like nucleation.But, the direct relevance of thin water films in coexistence with vapor to real physical scenarios is, in fact, somewhat limited; small volumes of liquid water will form droplets with high interfacial curvature, with an associated Laplace pressure.As discussed in Sec. 1, Li et al. 27 have found that differences in ice nucleation rate in small droplets relative to bulk water can largely be accounted for by this Laplace pressure.Notwithstanding the issues discussed in Sec. 3 in calculating the rate in small systems, we argue that this analysis of Li et al. is in line with our finding that J(W )/J(∞) ≈ 1 once the film is thick enough to accommodate a critical nucleus.Instead of extending our studies to small droplets, then, we instead focus on a system whose direct comparison to thin water films is more apparent: water confined between two 'inactive' walls.
A snapshot of the system we simulate is shown in Fig. 5a.To introduce confining walls, a slab comprising 9600 atoms arranged on a FCC lattice, and exposing its (111) face, is placed such that its lowermost plane of atoms is situated at z = H/2.An identical slab is then placed such that its uppermost plane of atoms is located at z = −H/2.The atoms belonging to these slabs interact with mW water molecules by a Lennard-Jones potential.The lattice constant of the slab and the parameters for the Lennard-Jones potential (see ESI †) were chosen such that the surface only promotes ice formation at temperatures below 201 K. 75 As we are investigating nucleation at T = 220 K, we can consider these surfaces to be inactive to ice nucleation, and we will assume that nucleation occurs homogeneously in the interior 'bulk-like' region of the fluid.Between these slabs, we introduce mW water molecules such that a 'squished cylinder' forms that spans one of the lateral dimensions, shown in Fig. 5a.In the orthogonal lateral dimension, a slightly convex water-vapor interface forms.The presence of a curved interface implies that, unlike the thin films considered above, the pressure inside the fluid is nonzero.In the films considered, the density far from the interfaces is virtually indistinguishable from the case of the thin water film (see ESI †).Note that the amount of water included in the simulation depends on H such that the curvature of the liquid-vapor interface is approximately constant (see ESI †).
Following the same procedure as we did in Sec. 3, we performed seeding simulations for H/nm = 4.0, 4.5, 5.0 and 6.0.In Fig. 5, we present p ice (n cl ; H), which shows that for H ≥ 4.5 nm, the size of the critical nucleus is insensitive to H and, within statistical uncertainty, indiscernible from n ‡ cl found in the thicker unsupported thin films.Also similar to the case of the unsupported films, for the smallest value of H investigated, n ‡ cl is seen to increase.Providing a rigorous relationship between H and W is beyond the scope of this work.Instead, in Fig. 6 we present density profiles of water for both the H = 4.5 nm confined system, and the W = 3.5 nm unsupported film (i.e., the smallest value of H and W for which n ‡ cl remains bulk-like).By eye, the extent of the bulk-like region in these two systems is comparable.Superposed on these plots is a snapshot (roughly to scale) of a critical nucleus for the W = 3.5 nm system, which shows that the diameter of the critical nucleus occupies the full extent of this bulk-like region.The insight gleaned from our investigation on thin films therefore also appears to hold in the more realistic case that water is confined between surfaces that do not promote ice formation.Fig. 6 The smallest values of W (unsupported film) and H (confined between walls) for which n ‡ cl remains bulk-like both exhibit similar sized bulk-like regions, as seen in ⟨ρ(z)⟩.For the system confined between walls, ⟨ρ(z)⟩ is calculated such that we only consider a region far from the convex liquid/vapor interface.Overlayed on these plots is a snapshot of a critical nucleus, roughly to scale, which shows that it spans essentially all of the central bulk-like region.

Conclusions
In this article, we have investigated ice nucleation in thin films of water, both freestanding, and confined between surfaces that do not promote ice nucleation.In both cases, our results show that the critical nucleus' size is indistinguishable from that of bulk water for sample sizes that can barely accommodate its presence.At a temperature of 220 K, once the thickness is decreased below approx.3.5 nm, we found that the size of the critical nucleus increases.These results for thin films are consistent with the observation in Ref. 27 that bulk-like nucleation is recovered in nanometer sized droplets, once the Laplace pressure is taken into account, and indeed, that the rate was seen to decrease relative to that of bulk for droplets too small to accommodate a critical nucleus.Ref. 27 and this work, however, have reached this con-clusion using a coarse grained representation of water's interactions (the mW model), while previous work 30 has found a qualitative discrepancy with more a more detailed model (TIP4P/ice) when ice nucleation rates in thin films are compared to those in homogeneous bulk.So, in this work, we have also shown that this apparent discrepancy between water models can be resolved through consistent treatment of truncated interactions between the homogeneous and inhomogeneous systems, when using more detailed models such as TIP4P/ice.
In the thin film systems we have investigated, we also observed that the attachment frequencies to the critical nuclei are comparable to those in bulk water.Calculating the rate based on classical nucleation theory then gives the impression that the nucleation rate in thin films is the same as bulk water, which is seemingly at odds with some, 25,29,30,33 but not all, 27 previous works.A simple model that supposes nucleation is suppressed in the interfacial regions is able to reconcile these differences at a qualitative level, but drawing quantitative conclusions on the rate (per unit volume per unit time) is made challenging by the sensitivity of the effective rate to the definition of the interfacial thickness.One scenario where one might be interested in the kinetics of ice formation under confinement is in the presence of pores or microstructures presented by solid particles such as silicas or feldspars i.e., pore condensation and freezing. 66Describing ice formation in these complex systems may benefit from simplified models, and to this end, the results of this work suggest CNT provides a reasonable foundation.
In any case, in experiments that investigate pore condensation and freezing, one typically measures the temperature at which ice forms rather than the nucleation rate directly.To gauge the impact on readily obtainable experimental observables, we therefore compared the fraction of frozen samples containing the same number of water molecules.Our results suggest that any difference in nucleation rate in a 3.5 nm thick film of water relative to that of a macroscopic sample is slight, and likely inconsequential for any experimental measurement that can be performed in the foreseeable future.They further suggest that the ice nucleation rate for water condensed in the pores of aerosol particles under cirrus cloud conditions will not be significantly suppressed due to the confined nature of the water in those pores.

Methods
Full details of the methodology are given in the ESI †.Liquid structures were initialized using the packmol software package, 76 and crystal structures with GenIce. 77,78All simulations were performed with the LAMMPS simulation package 79 using the standard velocity verlet algorithm. 80The temperature was controlled using a Nosé-Hoover thermostat 81 and, where necessary, the pressure maintained with a Parinello-Rahman barostat. 82Seeding systems were produced with the aid of the MDAnalysis software package. 83,84or TIP4P/ice, 38 simulations were performed in a similar manner, with electrostatic interactions computed with a particleparticle particle-mesh Ewald method, 85 with parameters chosen such that the root mean square error in the forces were a factor 10 5 smaller than the force between two unit charges separated by a distance of 1 Å. 86 For simulations of liquid water in contact with its vapor, we set D = 0, where D is the electric displacement field along z, using the implementation described in Refs.87 and 88 (this is formally equivalent to the commonly used Yeh and Berkowitz 89 method). Th rigid geometry of the molecules was maintained with the RATTLE algorithm.90 Parallel tempering 91,92 was employed for the TIP4P/ice (8.5) film in Sec. 2.
Steinhardt order parameters were calculated using the opensource, community-developed PLUMED library, 93 version 2.8. 94napshots were created in VMD. 95,96

Fig. 3
Fig. 3 Assessing the impact of film width, W , on ice nucleation.(a) The size of the critical cluster for a given W [as indicated in the legend in (c)] is obtained by finding the size of cluster after initial equilibration, n cl (t = 0) = n ‡ cl , for which there is equal probability to grow or shrink, p ice (n cl ;W ) = 0.5.We find that n ‡ cl only differs significantly from its value in bulk water for W ≲ 3.5 nm, as shown in (b).Errors bars indicate the range of results obtained by splitting each data set into three.In (c) we show ⟨(n cl (t) − n ‡ cl ) 2 ⟩ for each W .The attachment frequency f + is obtained by fitting a straight line after an initial transient period, as shown in bold for W = 3.5 nm.Dotted lines indicate hypothetical values of f + to give a sense of scale.Using the obtained n ‡ cl and f + , and ∆µ/k B = 122 K, the nucleation rate is obtained from Eq. 6, as shown in (d).In panels (b) and (d), the shaded region indicates the result obtained from simulations of bulk water.

5 (
f + is obtained by starting many simulations with n cl (t = 0) = n ‡ cl and fitting n cl (t) − n ‡ cl 2 to a straight line after an initial transient time.Despite significant noise in n cl (t) − n ‡ cl 2

Fig. 4
Fig. 4 Reconciling nucleation rates obtained from seeding with other approaches.(a) Schematic of the simple model described in the text.Nucleation is assumed to proceed in a bulk-like fashion in the central volume A(W − 2ℓ i − 2r ‡ cl ).In (b) we present the fraction of frozen samples for W /nm = 3.5 and 6, and bulk water.Each system comprises 6000 mW molecules and is cooled at a rate of 0.2 K/ns from an initial temperature of 220 K.There is a modest shift to lower temperatures with the films compared to bulk water.The freezing rate, R, for W = 3.5 nm, as shown in (c), differs at most by a factor five compared to the bulk result.

Fig. 5
Fig. 5 Conclusions drawn from unsupported thin films appear robust to the presence of physical boundaries.(a) Snapshot from a seeding simulation with a supported film confined between two slabs of Lennard-Jones particles separated by H = 4.5 nm.The water forms a 'squished cylinder' that spans the periodic boundaries of the simulation cell out of the plane of the page.In panel (b), we show p ice (n cl ; H), which shows the n ‡ cl remains bulk-like for H ≳ 4.5 nm.For H = 4 nm, the size of the critical nucleus increases, similar to what we observed in the unsupported films.