# The limit of macroscopic homogeneous ice nucleation at the nanoscale

John A. Hayton  Yusuf Hamied Department of Chemistry, University of Cambridge,
Lensfield Road, Cambridge CB2 1EW, United Kingdom  Michael B. Davies
Department of Physics and Astronomy, University College London, London WC1E
6BT, United Kingdom  Yusuf Hamied Department of Chemistry, University of
Cambridge, Lensfield Road, Cambridge CB2 1EW, United Kingdom  Thomas F. Whale
Department of Chemistry, University of Warwick, Gibbet Hill Road, Coventry CV4
7AL, United Kingdom.  Angelos Michaelides  Yusuf Hamied Department of
Chemistry, University of Cambridge, Lensfield Road, Cambridge CB2 1EW, United
Kingdom  Stephen J. Cox  Yusuf Hamied Department of Chemistry, University of
Cambridge, Lensfield Road, Cambridge CB2 1EW, United Kingdom<EMAIL_ADDRESS>
(22 June 2023)

###### Abstract

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 confined 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 films of water are affected by
the thickness of the film. Our results suggest that nucleation remains bulk-
like in films that are barely large enough accommodate a critical nucleus.
This conclusion seems robust to the presence of physical confining boundaries.
We also discuss the difficulties in unambiguously determining homogeneous
nucleation rates in nanoscale systems, owing to the challenges in defining the
volume. Our results suggest any impact on a film’s thickness on the rate is
largely inconsequential for present day experiments.

##  I Introduction

The formation of ice from liquid water is one of the most important phase
transitions on Earth, and plays a vital role in climate science,Cantrell and
Heymsfield (2005); Satoh _et al._ (2018); DeMott _et al._ (2010); Hoose _et
al._ (2010); Gettelman _et al._ (2010, 2012); Hudait and Molinero (2016)
cryopreservation,Kiani and Sun (2011); Rall and Fahy (1985); Reed (2014)
geologyGerrard (2012); Ashton (1986) and many industrial
applications.Rykaczewski _et al._ (2013); Fukasawa _et al._ (2001, 2002); Ohji
and Fukushima (2013) For example, many properties of clouds are affected by
the relative compositions of ice and water,Baker (1997); Carslaw _et al._
(2002); DeMott _et al._ (2010); Boucher _et al._ (2013) 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. $-40$°C, e.g., in cirrus cloud formation in the
upper troposphere.Kärcher and Lohmann (2002); Jensen _et al._ (1998); DeMott
_et al._ (1998); Rogers _et al._ (1998); DeMott _et al._ (2003); Satoh _et
al._ (2018); DeMott _et al._ (2010) 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 long-standing
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. Molecular simulations have
therefore been employed by several groups to investigate where in the liquid
homogeneous nucleation occurs.Vrbka and Jungwirth (2006, 2007); Li _et al._
(2013); Johnston and Molinero (2012); Haji-Akbari _et al._ (2014); Haji-Akbari
and Debenedetti (2017) While pioneering simulation studies from JungwirthVrbka
and Jungwirth (2006, 2007) 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.Cox _et al._ (2013) 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. These conclusions are
also supported by thermodynamic arguments made by Qiu and Molinero.Qiu and
Molinero (2018)

If 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._ Li _et al._ (2013) 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 $\gtrsim 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.Johnston and
Molinero (2012) However, subsequent simulation studies that investigated ice
nucleation in thin water films,Haji-Akbari _et al._ (2014); Lü _et al._ (2013)
whose planar interfaces correspond to zero Laplace pressure, found that
nucleation rates were noticeably suppressed relative to that of bulk, even for
film thicknesses $\gtrsim 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.Molinero and Moore (2009)
In 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. Despite its simplicity, the mW model reasonably
describes the anomalies and structure of water and its phase behavior,
including the density maximum, the increase in heat capacity and
compressibility of the supercooled region and melting temperatures of both
hexagonal and cubic ice.Molinero and Moore (2009); Moore and Molinero (2011);
Limmer and Chandler (2011); Moore _et al._ (2010)

To add further complication to the above apparent discrepancy between droplets
and films, tour-de-force simulations employing the TIP4P/ice model,Abascal _et
al._ (2005) 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.Haji-Akbari and Debenedetti (2017) This
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,Bi _et al._ (2017); Rosky _et al._ (2023); Moore _et al._ (2010) e.g.,
in explaining ice nucleation from vapor via a pore condensation and freezing
mechanism.David _et al._ (2019)

In 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. III,
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
vapor—impacts ice formation. In Sec. IV we extend our investigation to a more
realistic case where water is confined between two solid surfaces. We conclude
our findings in Sec. V. A broad overview of the methodology is described
throughout the manuscript, with full details provided in Sec. VI and the
Electronic Supporting Information (ESI).

##  II 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,Haji-Akbari
and Debenedetti (2017) using the TIP4P/ice water model at 230 K has reported
that the nucleation rate $J(W)$ in a film of thickness $W\approx 4$ nm is
roughly seven orders of magnitude larger than in bulk water, whose rate we
denote $J(\infty)$. 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 parameterSteinhardt _et al._ (1983) for the $i^{\mathrm{th}}$ water
molecule,

$\bar{Q}_{6}^{(i)}=\frac{Q_{6}^{(i)}+\sideset{}{{}^{\prime}}{\sum}_{j}Q_{6}^{(j)}}{\nu^{(i)}+1},$
(1)

indicate that the average structure of liquid water is more ice-like for thin
films when compared to homogeneous bulk water under the same
conditions.111There a several similar measures of local order in the
literature,Li _et al._ (2011); Lechner and Dellago (2008); Rein ten Wolde _et
al._ (1996) but the analysis we present in this section is robust to the exact
choice.  The prime indicates that the sum only includes the $\nu^{(i)}$
nearest neighbors of molecule $i$ (see ESI),

$Q_{6}^{(i)}=\frac{1}{\nu^{(i)}}\sqrt{\sum_{m=-6}^{6}\sideset{}{{}^{\prime}}{\sum}_{j,k}Y_{6m}^{\ast}(\hat{\mathbf{r}}_{ij})Y_{6m}(\hat{\mathbf{r}}_{ik})},$
(2)

where $\hat{\mathbf{r}}_{ij}$ is the unit vector pointing from molecule $i$ to
molecule $j$, and $Y_{6m}$ is the $m^{\rm 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._ Atherton _et al._ (2022) discussed the
sensitivity of the melting point of simple point charge models such as
TIP4P/ice to the choice of cutoff, $r_{\ast}$, for the non-electrostatic
interactions between water molecules. There it was shown that decreasing
$r_{\ast}$ led to a systematic increase in the melting temperature. More
importantly, it was noticed that choices of $r_{\ast}$ 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._ Bianco _et al._ (2021) 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.Rowlinson and Widom (1982) Therefore, 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

$\langle\bar{q}_{6}\rangle/\bar{\rho}=\frac{1}{N}\big{\langle}\sum_{i}\bar{Q}_{6}^{(i)}\big{\rangle},$
(3)

where $\bar{\rho}$ 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 non-electrostatic interactions at
$r_{\ast}=8.5\,\text{\AA}$; following Ref. 46 we indicate this molecular model
as TIP4P/ice(8.5). We find $\langle\bar{q}_{6}\rangle/\bar{\rho}\approx
0.266$.

Figure 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\approx 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
$\langle\bar{q}_{6}(z)\rangle/\langle\rho(z)\rangle$ in (b) for TIP4P/ice(8.5)
(blue) and mW (orange). Dashed lines show
$\langle\bar{q}_{6}\rangle/\bar{\rho}$ 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=\pm W/2$.

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\approx 4.3$ nm. Throughout this paper, the width of a film is defined
as $W=N/A\bar{\rho}$, where $A$ is 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

$\langle\bar{q}_{6}(z)\rangle/\langle\rho(z)\rangle=\frac{1}{A\langle\rho(z)\rangle}\big{\langle}\sum_{i}\bar{Q}_{6}^{(i)}\delta(z-z_{i})\big{\rangle},$
(4)

where $z$ is the coordinate normal to the average liquid-vapor interface,
$z_{i}$ is the $z$ coordinate of $i^{\rm th}$ water molecule’s oxygen atom,
and

$\langle\rho(z)\rangle=\frac{1}{A}\big{\langle}\sum_{i}\delta(z-z_{i})\big{\rangle}$
(5)

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
$\langle\bar{q}_{6}(z)\rangle/\langle\rho(z)\rangle\approx 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
$\langle\bar{q}_{6}\rangle/\langle\rho(z)\rangle\approx 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
$\langle\bar{q}_{6}(z)\rangle/\langle\rho(z)\rangle$ 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 $\langle\bar{q}_{6}\rangle/\langle\rho(z)\rangle\approx 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,Haji-Akbari _et al._ (2014); Haji-Akbari and Debenedetti (2017); Li _et
al._ (2013) 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.

##  III Nucleation in thin water films remains bulk-like down to very small
length scales

Figure 2: A snapshot from a seeding simulation with $W=3.5$ nm after initial
equilibration at $T=220$ K, using the mW model. Molecules identified as
belonging to the largest ice-like cluster are shown by large blue spheres,
while all other molecules are shown by small red spheres. In this case, the
initial seed happens to be critical, with approximate radius $r_{\rm
cl}^{\ddagger}$, as indicated. Figure 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_{\rm cl}(t=0)=n_{\rm cl}^{\ddagger}$, for
which there is equal probability to grow or shrink, $p_{\rm ice}(n_{\rm
cl};W)=0.5$. We find that $n_{\rm cl}^{\ddagger}$ only differs significantly
from its value in bulk water for $W\lesssim 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 $\langle(n_{\rm cl}(t)-n_{\rm
cl}^{\ddagger})^{2}\rangle$ 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_{\rm cl}^{\ddagger}$
and $f^{+}$, and $\Delta\mu/k_{\rm 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.

Having established consistency between TIP4P/ice(8.5) and mW in describing the
structure of thin water films, we will now investigate the extent to which $W$
impacts ice nucleation. To do so, we will make use of the ‘seeding’ method,
first introduced by Bai and LiBai and Li (2005, 2006) in their study of
Lennard-Jones particles, and popularized for ice nucleation by Vega and co-
workers. Sanz _et al._ (2013); Espinosa _et al._ (2016); Zaragoza _et al._
(2015); Espinosa _et al._ (2014, 2016, 2016, 2018); Bianco _et al._ (2021) As
the seeding approach has been detailed previously by others, Sanz _et al._
(2013); Espinosa _et al._ (2014, 2016, 2016) here we will only briefly sketch
an outline of the procedure.

###  III.1 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_{\rm cl}$, which we
took to define the size of our initial ice seed. An example of such a seed is
shown in Fig. 2.

rowsep=4pt

$W\,\mathrm{[nm]}$ | $n_{\mathrm{cl}}^{{\ddagger}}$ | $f^{+}\,\mathrm{[ps^{-1}]}$ | $\log_{10}\left(J/(\mathrm{m^{-3}\>s^{-1}})\right)$
---|---|---|---
$\infty$ | $122.5\>(117.5-127.5)$ | $13.7$ | $25.1\>(24.4-25.7)$
6.0 | $119.5\>(113.5-126.5)$ | $11.9$ | $25.4\>(24.5-26.2)$
5.0 | $121.5\>(121.5-124.5)$ | $9.3$ | $25.1\>(24.6-25.1)$
4.0 | $121.5\>(119.5-125.5)$ | $14.2$ | $25.2\>(24.7-25.5)$
3.5 | $125.5\>(120.5-126.5)$ | $13.1$ | $24.7\>(24.5-25.4)$
3.0 | $132.5\>(129.5-134.5)$ | $9.7$ | $23.7\>(23.4-24.1)$
2.75 | $143.5\>(139.5-153.5)$ | $9.5$ | $22.4\>(21.1-22.9)$
2.5 | $172.5\>(166.5-187.5)$ | $14.2$ | $19.0\>(17.1-19.8)$
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_{\mathrm{cl}}^{{\ddagger}}$ for a given $W$
in Eq. 6 ($\Delta\mu/k_{\rm B}=122$ K and $\bar{\rho}=33.3774$ nm-3), and is
also reported in parentheses

We have investigated film thicknesses in the range $2.5\,\mathrm{nm}\lesssim
W\lesssim 6\,\mathrm{nm}$, with $N\approx 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\lesssim n_{\rm cl}\lesssim 220$ overall; this
was sufficient to span both pre- and post-critical cluster sizes. In Fig. 3a,
we present the probability $p_{\rm ice}(n_{\rm cl};W)$ that a seed of size
$n_{\rm cl}$ goes on to form ice. The critical cluster size $n_{\rm
cl}^{\ddagger}$ is estimated from $p_{\rm ice}(n_{\rm cl}^{\ddagger};W)=0.5$.
We observe that for $W\gtrsim 3.5$ nm, $n_{\rm cl}^{\ddagger}\approx 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_{\rm cl}^{\ddagger}\approx 0.96$ nm. If we consider $W=3.5$ nm, this gives
$(W-2r_{\rm cl}^{\ddagger})/2\approx 0.8$ nm, which is broadly in line with
accepted thicknesses of the liquid-vapor interface.Haji-Akbari _et al._
(2014); Li _et al._ (2013); Sega _et al._ (2016) In 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\lesssim 3$ nm, deviations of $n_{\rm cl}^{\ddagger}$ from its bulk value
are observed, with $n_{\rm cl}^{\ddagger}\approx 170$ in the thinnest film
investigated ($W\approx 2.5$ nm). For all systems investigated, $n_{\rm
cl}^{\ddagger}$ vs. $W$ is plotted in Fig. 3b.

Computation of nucleation rates with the seeding approach relies upon
classical nucleation theory (CNT):Sanz _et al._ (2013); Espinosa _et al._
(2016); Bai and Li (2005, 2006)

$J=\mathcal{Z}f^{+}\bar{\rho}\exp\big{(}-\beta|\Delta\mu|n_{\rm
cl}^{\ddagger}/2\big{)},$ (6)

where $\Delta\mu$ is the chemical potential difference between bulk ice and
liquid water, $\beta=1/k_{\rm B}T$ ($k_{\rm B}$ is Boltzmann’s constant),
$\mathcal{Z}=(\beta|\Delta\mu|/6\pi n_{\rm cl}^{\ddagger})^{1/2}$ is the
Zeldovich factor, and $f^{+}$ is the attachment frequency. For the chemical
potential difference, a range of values spanning approx.
$118\,\mathrm{K}\lesssim|\Delta\mu|/k_{\rm B}\lesssim 126\,\mathrm{K}$ have
been reported in the literature Reinhardt and Doye (2012); Jacobson _et al._
(2009) and for simplicity, we take $|\Delta\mu|/k_{\rm 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_{\rm cl}^{\ddagger}$ remains roughly constant for $W\gtrsim 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

$f^{+}=\frac{\big{\langle}\big{(}n_{\rm cl}(t)-n_{\rm
cl}^{\ddagger}\big{)}^{2}\big{\rangle}}{2t},$ (7)

as first proposed by Auer and Frenkel.Auer and Frenkel (2001, 2004) In
practice, $f^{+}$ is obtained by starting many simulations with $n_{\rm
cl}(t=0)=n_{\rm cl}^{\ddagger}$ and fitting $\big{\langle}\big{(}n_{\rm
cl}(t)-n_{\rm cl}^{\ddagger}\big{)}^{2}\big{\rangle}$ to a straight line after
an initial transient time. Despite significant noise in
$\big{\langle}\big{(}n_{\rm cl}(t)-n_{\rm
cl}^{\ddagger}\big{)}^{2}\big{\rangle}$, 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(\infty)$—obtained from a simulation of bulk water under periodic boundary
conditions—agrees well with that previously reported by Sanchez-Burgos _et
al._ Sanchez-Burgos _et al._ (2022) 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._ Lupi
_et al._ (2017) 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.

###  III.2 Comparing, and reconciling, results from seeding with previous
studies

Despite good agreement of $J(\infty)$ with previous literature values, our
results for $W\sim 5$ nm are at odds with previous work. For example, using
forward flux sampling, for mW under the same conditions, Haji-Akbari _et al._
Haji-Akbari _et al._ (2014) found $J(5\,{\rm nm})/J(\infty)\approx 0.006$.
Similarly, by computing mean first passage times at 205 K, Lü _et al._ Lü _et
al._ (2013) 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(\infty)\approx
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\,{\rm nm})/J(\infty)\approx 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_{\rm
cl}^{\ddagger}$ and $f^{+}$ remain roughly constant for $W\geq 3.5$ nm. We
note that in the supporting information of Li _et al._ ,Li _et al._ (2013) 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.

Figure 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\ell_{\rm i}-2r_{\rm cl}^{\ddagger})$. In (b) we present the fraction of
frozen samples for $W/{\rm 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.

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.Lü _et
al._ (2013); Li _et al._ (2013); Haji-Akbari and Debenedetti (2017) We suppose
that the thin film can be separated into two interfacial regions of thickness
$\ell_{\rm i}$ that sandwich a central bulk-like region of thickness
$W-2\ell_{\rm i}$. Accounting for the volume occupied by the critical nucleus,
the volume accessible for nucleation in the bulk-like region is
$A(W-2\ell_{\rm i}-2r_{\rm cl}^{\ddagger})$ (see Fig. 4). In the remaining
volume, $2A(\ell_{\rm i}+r_{\rm cl}^{\ddagger})$, we assume that nucleation
occurs with a rate $J_{\rm 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_{\rm eff}$:

$AWJ_{\rm eff}(W)=2A(\ell_{\rm i}+r_{\rm cl}^{\ddagger})J_{\rm
i}+A(W-2\ell_{\rm i}-2r_{\rm cl}^{\ddagger})J(\infty).$ (8)

Rearranging for $J_{\rm eff}$, and assuming that $J_{\rm i}\approx 0$, we find

$J_{\rm eff}(W)=\bigg{[}1-\frac{2(\ell_{\rm i}+r_{\rm
cl}^{\ddagger})}{W}\bigg{]}J(\infty).$ (9)

Using $r_{\rm cl}^{\ddagger}\approx 0.96$ nm obtained from our seeding
simulations, a suppression of the nucleation rate of $J_{\rm
eff}(5\,\mathrm{nm})/J(\infty)\approx 0.006$, as obtained from forward flux
sampling,Haji-Akbari _et al._ (2014) requires $\ell_{\rm i}\approx 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\gtrsim 3.5$ nm. As we did above, we can instead estimate
$\ell_{\rm i}=(3.5\,{\rm nm}-2r_{\rm cl}^{\ddagger})/2\approx 0.8$ nm. For
$W=5$ nm, we then obtain $J_{\rm eff}(5\,\mathrm{nm})/J(\infty)\approx 0.3$.
This more modest suppression in the rate appears comparable to that obtained
at 205 K by Lü _et al._ Lü _et al._ (2013)

###  III.3 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 $\ell_{\rm i}$ and the size of
the critical nucleus $r_{\rm cl}^{\ddagger}$. 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_{\rm
eff}$. In an attempt to gauge the significance of any suppression of $J_{\rm
eff}$ with $W$, it is useful to think of an instance when one might be
interested in ice formation in such small volumes, such as to understand how a
porous material (or a rough surface with, e.g., cracks) might promote ice
formation via a pore condensation and freezing mechanism.David _et al._
(2019); Holden _et al._ (2021); Campbell and Christenson (2018); Whale _et
al._ (2017); Pach and Verdaguer (2019) 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\,\mathrm{nm}$ and $W=3.5\,\mathrm{nm}$,
obtained by cooling 100 replicas of each system at $0.2\,\mathrm{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_{\mathrm{cl}}(T)\geq 200$. 222In
rare cases that the largest nucleus recrossed 200 molecules, the lower
temperature for which $n_{\mathrm{cl}}(T)$ passed 200 molecules was taken.
Such recrossings occurred in fewer than 1% of our simulations.  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\,\mathrm{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.Vali _et al._ (2015) $R$ is
analogous to the radioactive decay constant (a first-order reaction rate
constant) and is most easily determined from experiments conducted on
completely identical supercooled droplets under isothermal conditions. In the
isothermal case, the fraction of droplets which freeze in a given time period
can be found as:Koop _et al._ (1997)

$\frac{\mathcal{N}_{\mathrm{liq}}(t)}{\mathcal{N}_{\rm tot}}=\exp(-Rt)$ (10)

where $\mathcal{N}_{\mathrm{liq}}(t)$ is the number of liquid droplets
remaining after time $t$ and $\mathcal{N}_{\rm tot}$ is the total number of
droplets (liquid or solid) present in the experiment [equal to
$\mathcal{N}_{\mathrm{liq}}(t=0)$]. In order to calculate $R$ from
continuously cooled experiments, it is necessary to divide the experiment into
small time intervals $\Delta t$ over which changes in temperature are treated
as small. Rearranging Eqn. 10 we find

$R=\frac{-\log\left(1-\mathcal{N}_{\mathrm{f}}/\mathcal{N}_{\mathrm{i}}\right)}{\Delta
t}$ (11)

where $\mathcal{N}_{\mathrm{i}}$ is the number of unfrozen droplets present at
the start of the time interval and $\mathcal{N}_{\mathrm{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 $\Delta t=5\times 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.StataCorp  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.Koop _et al._ (1997); Blackman _et al._ (2017) For 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 5th to 95th percentile range of the simulated experiments. (Data
points generated from bins containing so few freezing events that the 5th
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\,\mathrm{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/\bar{\rho}\approx
3\times 10^{-29}\,\mathrm{m^{3}}$ we find that 6000 water molecules occupy
$1.8\times 10^{-25}\,\mathrm{m^{3}}$. Dividing our calculated freezing rates
by this volume we obtain values of $J$ entirely consistent with previous
literature values Haji-Akbari _et al._ (2014); Russo _et al._ (2014); Moore
and Molinero (2011); Sanchez-Burgos _et al._ (2022) 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.

##  IV Assessing the impact of physical boundaries

Figure 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_{\mathrm{ice}}(n_{\mathrm{cl}};H)$, which shows the $n_{\rm
cl}^{\ddagger}$ remains bulk-like for $H\gtrsim 4.5$ nm. For $H=4$ nm, the
size of the critical nucleus increases, similar to what we observed in the
unsupported films.

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. I, Li _et al._ Li _et al._ (2013) 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. III 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(\infty)\approx 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\,\mathrm{K}$.Davies _et al._ (2022) 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).

Figure 6: The smallest values of $W$ (unsupported film) and $H$ (confined
between walls) for which $n_{\rm cl}^{\ddagger}$ remains bulk-like both
exhibit similar sized bulk-like regions, as seen in $\langle\rho(z)\rangle$.
For the system confined between walls, $\langle\rho(z)\rangle$ 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.

Following the same procedure as we did in Sec. III, we performed seeding
simulations for $H/{\rm nm}=4.0,4.5,5.0$ and $6.0$. In Fig. 5, we present
$p_{\rm ice}(n_{\rm cl};H)$, which shows that for $H\geq 4.5$ nm, the size of
the critical nucleus is insensitive to $H$ and, within statistical
uncertainty, indiscernible from $n_{\rm cl}^{\ddagger}$ found in the thicker
unsupported thin films. Also similar to the case of the unsupported films, for
the smallest value of $H$ investigated, $n_{\rm cl}^{\ddagger}$ 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_{\rm cl}^{\ddagger}$
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.

##  V 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 conclusion using a coarse grained representation of water’s interactions
(the mW model), while previous workHaji-Akbari and Debenedetti (2017) 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,Lü _et al._ (2013); Haji-Akbari
_et al._ (2014); Haji-Akbari and Debenedetti (2017); Vrbka and Jungwirth
(2006) but not all,Li _et al._ (2013) 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.Holden _et al._ (2021) Describing 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.

##  VI Methods

Full details of the methodology are given in the ESI. Liquid structures were
initialized using the packmol software package,Martínez _et al._ (2009) and
crystal structures with GenIce.Matsumoto _et al._ (2017, 2021) All simulations
were performed with the LAMMPS simulation packageThompson _et al._ (2022)
using the standard velocity verlet algorithm.Swope _et al._ (1982) The
temperature was controlled using a Nosé-Hoover thermostatNosé (2006) and,
where necessary, the pressure maintained with a Parinello-Rahman
barostat.Parrinello and Rahman (1981) Seeding systems were produced with the
aid of the MDAnalysis software package.Michaud-Agrawal _et al._ (2011);
Richard J. Gowers _et al._ (2016)

For TIP4P/ice,Abascal _et al._ (2005) simulations were performed in a similar
manner, with electrostatic interactions computed with a particle-particle
particle-mesh Ewald method,Hockney and Eastwood (2021) 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
Å.Kolafa and Perram (1992) 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 BerkowitzYeh and Berkowitz (1999)
method). The rigid geometry of the molecules was maintained with the RATTLE
algorithm.Andersen (1983) Parallel temperingSwendsen and Wang (1986); Earl and
Deem (2005) was employed for the TIP4P/ice(8.5) film in Sec. II.

Steinhardt order parameters were calculated using the open-source, community-
developed PLUMED library,Bonomi _et al._ (2019) version 2.8.Tribello _et al._
(2014) Snapshots were created in VMD.Humphrey _et al._ (1996); Stone (1998)

## Data Availability Statement

The data that supports the findings of this study and input files for the
simulations, are openly available at the University of Cambridge Data
Repository, https://doi.org/10.17863/CAM.96642.

## Author Contributions

John A. Hayton: Data Curation (lead), Formal Analysis (lead), Investigation
(lead), Methodology (lead), Writing - Original Draft (lead), Writing - Review
and Editing (contributing)

Michael B. Davies: Conceptualization (contributing), Supervision
(contributing), Writing - Review and Editing (contributing)

Thomas F. Whale: Formal Analysis (contributing), Writing - Review and Editing
(contributing)

Angelos Michaelides: Conceptualization (contributing), Writing - Review and
Editing (contributing)

Stephen J. Cox: Conceptualization (lead), Supervision (lead), Investigation
(contributing) Writing - Review and Editing (lead)

## Conflicts of interest

There are no conflicts to declare.

## Acknowledgements

J.A.H. acknowledges studentship funding from the Engineering and Physical
Sciences Research Council (EP/T517847/1). T.F.W. thanks the Leverhulme Trust
and the University of Warwick for supporting an Early Career Fellowship
(ECF-2018-127). S.J.C. is a Royal Society University Research Fellow (Grant
No. URF\R1\211144) at the University of Cambridge.

## References

  * Cantrell and Heymsfield (2005) W. Cantrell and A. Heymsfield, _Bull Am Meteorol Soc_ , 2005, 86, 795–808.
  * Satoh _et al._ (2018) M. Satoh, A. T. Noda, T. Seiki, Y. W. Chen, C. Kodama, Y. Yamada, N. Kuba and Y. Sato, _Prog. Earth Planet. Sci._ , 2018, 5, 1–29.
  * DeMott _et al._ (2010) P. J. DeMott, A. J. Prenni, X. Liu, S. M. Kreidenweis, M. D. Petters, C. H. Twohy, M. S. Richardson, T. Eidhammer and D. C. Rogers, _Proc. Natl. Acad. Sci. U.S.A._ , 2010, 107, 11217–11222.
  * Hoose _et al._ (2010) C. Hoose, J. E. Kristjánsson, J. P. Chen and A. Hazra, _J. Atmos. Sci._ , 2010, 67, 2483–2503.
  * Gettelman _et al._ (2010) A. Gettelman, X. Liu, S. J. Ghan, H. Morrison, S. Park, A. J. Conley, S. A. Klein, J. Boyle, D. L. Mitchell and J. L. Li, _J. Geophys. Res. Atmos._ , 2010, 115, .
  * Gettelman _et al._ (2012) A. Gettelman, X. Liu, D. Barahona, U. Lohmann and C. Chen, _J. Geophys. Res. Atmos._ , 2012, 117, 20201.
  * Hudait and Molinero (2016) A. Hudait and V. Molinero, _J. Am. Chem. Soc._ , 2016, 138, 8958–8967.
  * Kiani and Sun (2011) H. Kiani and D. W. Sun, _Trends Food Sci. Technol._ , 2011, 22, 407–426.
  * Rall and Fahy (1985) W. F. Rall and G. M. Fahy, _Nature_ , 1985, 313, 573–575.
  * Reed (2014) B. M. Reed, II International Symposium on Plant Cryopreservation 1039, 2014, pp. 41–48.
  * Gerrard (2012) J. Gerrard, _Rocks and landforms_ , Springer Science & Business Media, 2012\.
  * Ashton (1986) G. D. Ashton, _River and lake ice engineering_ , Water Resources Publication, 1986.
  * Rykaczewski _et al._ (2013) K. Rykaczewski, S. Anand, S. B. Subramanyam and K. K. Varanasi, _Langmuir_ , 2013, 29, 5230–5238.
  * Fukasawa _et al._ (2001) T. Fukasawa, M. Ando, T. Ohji and S. Kanzaki, _J. Am. Ceram. Soc._ , 2001, 84, 230–232.
  * Fukasawa _et al._ (2002) T. Fukasawa, Z. Y. Deng, M. Ando, T. Ohji and S. Kanzaki, _J. Am. Ceram. Soc._ , 2002, 85, 2151–2155.
  * Ohji and Fukushima (2013) T. Ohji and M. Fukushima, _Int. Mater. Rev._ , 2013, 57, 115–131.
  * Baker (1997) M. Baker, _Science_ , 1997, 276, 1072–1078.
  * Carslaw _et al._ (2002) K. Carslaw, R. Harrison and J. Kirkby, _Science_ , 2002, 298, 1732–1737.
  * Boucher _et al._ (2013) O. Boucher, D. Randall, P. Artaxo, C. Bretherton, G. Feingold, P. Forster, V.-M. Kerminen, Y. Kondo, H. Liao, U. Lohmann _et al._ , _Climate change 2013: the physical science basis. Contribution of Working Group I to the Fifth Assessment Report of the Intergovernmental Panel on Climate Change_ , Cambridge University Press, 2013, pp. 571–657.
  * Kärcher and Lohmann (2002) B. Kärcher and U. Lohmann, _J. Geophys. Res. Atmos._ , 2002, 107, AAC 4–1.
  * Jensen _et al._ (1998) E. J. Jensen, O. B. Toon, A. Tabazadeh, G. W. Sachse, B. E. Anderson, K. R. Chan, C. W. Twohy, B. Gandrud, S. M. Aulenbach, A. Heymsfield, J. Hallett and B. Gary, _Geophys. Res. Lett._ , 1998, 25, 1363–1366.
  * DeMott _et al._ (1998) P. J. DeMott, D. C. Rogers, S. M. Kreidenweis, Y. Chen, C. H. Twohy, D. Baumgardner, A. J. Heymsfield and K. R. Chan, _Geophys. Res. Lett._ , 1998, 25, 1387–1390.
  * Rogers _et al._ (1998) D. C. Rogers, P. J. DeMott, S. M. Kreidenweis and Y. Chen, _Geophys. Res. Lett._ , 1998, 25, 1383–1386.
  * DeMott _et al._ (2003) P. J. DeMott, D. J. Cziczo, A. J. Prenni, D. M. Murphy, S. M. Kreidenweis, D. S. Thomson, R. Borys and D. C. Rogers, _Proc. Natl. Acad. Sci. U.S.A._ , 2003, 100, 14655–14660.
  * Vrbka and Jungwirth (2006) L. Vrbka and P. Jungwirth, _J. Phys. Chem. B_ , 2006, 110, 18126–18129.
  * Vrbka and Jungwirth (2007) L. Vrbka and P. Jungwirth, _J. Mol. Liq._ , 2007, 134, 64–70.
  * Li _et al._ (2013) T. Li, D. Donadio and G. Galli, _Nat. Commun._ , 2013, 4, 1–6.
  * Johnston and Molinero (2012) J. C. Johnston and V. Molinero, _J. Am. Chem. Soc._ , 2012, 134, 6650–6659.
  * Haji-Akbari _et al._ (2014) A. Haji-Akbari, R. S. DeFever, S. Sarupria and P. G. Debenedetti, _Phys. Chem. Chem. Phys._ , 2014, 16, 25916–25927.
  * Haji-Akbari and Debenedetti (2017) A. Haji-Akbari and P. G. Debenedetti, _Proc. Natl. Acad. Sci. U.S.A._ , 2017, 114, 3316–3321.
  * Cox _et al._ (2013) S. J. Cox, Z. Raza, S. M. Kathmann, B. Slater and A. Michaelides, _Faraday Discuss._ , 2013, 167, 389–403.
  * Qiu and Molinero (2018) Y. Qiu and V. Molinero, _J. Phys. Chem. Lett._ , 2018, 9, 5179–5182.
  * Lü _et al._ (2013) Y. Lü, X. Zhang and M. Chen, _J. Phys. Chem. B_ , 2013, 117, 10241–10249.
  * Molinero and Moore (2009) V. Molinero and E. B. Moore, _J. Phys. Chem. B_ , 2009, 113, 4008–4016.
  * Moore and Molinero (2011) E. B. Moore and V. Molinero, _Nature_ , 2011, 479, 506–508.
  * Limmer and Chandler (2011) D. T. Limmer and D. Chandler, _J. Chem. Phys._ , 2011, 135, 134503\.
  * Moore _et al._ (2010) E. B. Moore, E. De La Llave, K. Welke, D. A. Scherlis and V. Molinero, _Phys. Chem. Chem. Phys._ , 2010, 12, 4124–4134.
  * Abascal _et al._ (2005) J. L. Abascal, E. Sanz, R. G. Fernández and C. Vega, _J. Chem. Phys._ , 2005, 122, 234511.
  * Bi _et al._ (2017) Y. Bi, B. Cao and T. Li, _Nat. Commun._ , 2017, 8, 1–7.
  * Rosky _et al._ (2023) E. Rosky, W. Cantrell, T. Li, I. Nakamura and R. A. Shaw, _EGUsphere_ , 2023, 2023, 1–21.
  * David _et al._ (2019) R. O. David, C. Marcolli, J. Fahrni, Y. Qiu, Y. A. P. Sirkin, V. Molinero, F. Mahrt, D. Brühwiler, U. Lohmann and Z. A. Kanji, _Proc. Natl. Acad. Sci. U.S.A._ , 2019, 116, 8184–8189.
  * Steinhardt _et al._ (1983) P. J. Steinhardt, D. R. Nelson and M. Ronchetti, _Phys. Rev. B_ , 1983, 28, 784.
  * Li _et al._ (2011) T. Li, D. Donadio, G. Russo and G. Galli, _Phys. Chem. Chem. Phys._ , 2011, 13, 19807–19813.
  * Lechner and Dellago (2008) W. Lechner and C. Dellago, _J. Chem. Phys._ , 2008, 129, 114707\.
  * Rein ten Wolde _et al._ (1996) P. Rein ten Wolde, M. J. Ruiz-Montero and D. Frenkel, _J. Chem. Phys._ , 1996, 104, 9932–9947.
  * Atherton _et al._ (2022) D. Atherton, A. Michaelides and S. J. Cox, _J. Chem. Phys._ , 2022, 156, 164501.
  * Bianco _et al._ (2021) V. Bianco, P. M. D. Hijes, C. P. Lamas, E. Sanz and C. Vega, _Phys. Rev. Lett._ , 2021, 126, 15704.
  * Rowlinson and Widom (1982) J. S. Rowlinson and B. Widom, _Molecular theory of capillarity_ , Clarendon Press, 1982.
  * Bai and Li (2005) X. M. Bai and M. Li, _J. Chem. Phys._ , 2005, 122, 224510.
  * Bai and Li (2006) X. M. Bai and M. Li, _J. Chem. Phys._ , 2006, 124, 124707.
  * Sanz _et al._ (2013) E. Sanz, C. Vega, J. R. Espinosa, R. Caballero-Bernal, J. L. Abascal and C. Valeriani, _J. Am. Chem. Soc._ , 2013, 135, 15008–15017.
  * Espinosa _et al._ (2016) J. R. Espinosa, C. Vega, C. Valeriani and E. Sanz, _J. Chem. Phys._ , 2016, 144, 034501.
  * Zaragoza _et al._ (2015) A. Zaragoza, M. M. Conde, J. R. Espinosa, C. Valeriani, C. Vega and E. Sanz, _J. Chem. Phys._ , 2015, 143, 134504.
  * Espinosa _et al._ (2014) J. R. Espinosa, E. Sanz, C. Valeriani and C. Vega, _J. Chem. Phys._ , 2014, 141, 18C529.
  * Espinosa _et al._ (2016) J. R. Espinosa, C. Navarro, E. Sanz, C. Valeriani and C. Vega, _J. Chem. Phys._ , 2016, 145, 211922.
  * Espinosa _et al._ (2016) J. R. Espinosa, A. Zaragoza, P. Rosales-Pelaez, C. Navarro, C. Valeriani, C. Vega and E. Sanz, _Phys. Rev. Lett._ , 2016, 117, 135702\.
  * Espinosa _et al._ (2018) J. R. Espinosa, C. Vega and E. Sanz, _J. Phys. Chem. C_ , 2018, 122, 22892–22896.
  * Sega _et al._ (2016) M. Sega, B. Fábián, G. Horvai and P. Jedlovszky, _J. Phys. Chem. C_ , 2016, 120, 27468–27477.
  * Reinhardt and Doye (2012) A. Reinhardt and J. P. Doye, _J. Chem. Phys._ , 2012, 136, 054501\.
  * Jacobson _et al._ (2009) L. C. Jacobson, W. Hujo and V. Molinero, _J. Phys. Chem. B_ , 2009, 113, 10298–10307.
  * Auer and Frenkel (2001) S. Auer and D. Frenkel, _Nature_ , 2001, 409, 1020–1023.
  * Auer and Frenkel (2004) S. Auer and D. Frenkel, _J. Chem. Phys._ , 2004, 120, 3015.
  * Sanchez-Burgos _et al._ (2022) I. Sanchez-Burgos, A. R. Tejedor, C. Vega, M. M. Conde, E. Sanz, J. Ramirez and J. R. Espinosa, _J. Chem. Phys._ , 2022, 157, 094503.
  * Lupi _et al._ (2017) L. Lupi, A. Hudait, B. Peters, M. Grünwald, R. Gotchy Mullen, A. H. Nguyen and V. Molinero, _Nature_ , 2017, 551, 218–222.
  * Haji-Akbari and Debenedetti (2017) A. Haji-Akbari and P. G. Debenedetti, _J. Chem. Phys._ , 2017, 147, 060901.
  * Holden _et al._ (2021) M. A. Holden, J. M. Campbell, F. C. Meldrum, B. J. Murray and H. K. Christenson, _Proc. Natl. Acad. Sci. U.S.A._ , 2021, 118, e2022859118.
  * Campbell and Christenson (2018) J. M. Campbell and H. K. Christenson, _Phys. Rev. Lett._ , 2018, 165701\.
  * Whale _et al._ (2017) T. F. Whale, M. A. Holden, A. N. Kulak, Y. Y. Kim, F. C. Meldrum, H. K. Christenson and B. J. Murray, _Phys. Chem. Chem. Phys._ , 2017, 19, 31186–31193.
  * Pach and Verdaguer (2019) E. Pach and A. Verdaguer, _J. Phys. Chem. C_ , 2019, 123, 20998–21004.
  * Vali _et al._ (2015) G. Vali, P. DeMott, O. Möhler and T. Whale, _Atmos. Chem. Phys._ , 2015, 15, 10263–10270.
  * Koop _et al._ (1997) T. Koop, B. Luo, U. M. Biermann, P. J. Crutzen and T. Peter, _J. Phys. Chem. A_ , 1997, 101, 1117–1133.
  * (72) StataCorp, _Stata Statistical Software: Release 17. College Station, TX: StataCorp LLC._ , https://www.stata.com/.
  * Blackman _et al._ (2017) L. D. Blackman, S. Varlas, M. C. Arno, A. Fayter, M. I. Gibson and R. K. O’Reilly, _ACS Macro Lett._ , 2017, 6, 1263–1267.
  * Russo _et al._ (2014) J. Russo, F. Romano and H. Tanaka, _Nat. Matter._ , 2014, 13, 733–739.
  * Davies _et al._ (2022) M. B. Davies, M. Fitzner and A. Michaelides, _Proc. Natl. Acad. Sci. U.S.A._ , 2022, 119, e2205347119.
  * Martínez _et al._ (2009) L. Martínez, R. Andrade, E. G. Birgin and J. M. Martínez, _J Comput Chem_ , 2009, 30, 2157–2164.
  * Matsumoto _et al._ (2017) M. Matsumoto, T. Yagasaki and H. Tanaka, _J Comput Chem_ , 2017, 39, 61–64.
  * Matsumoto _et al._ (2021) M. Matsumoto, T. Yagasaki and H. Tanaka, _J Chem Inf Model_ , 2021, 61, 2542–2546.
  * Thompson _et al._ (2022) A. P. Thompson, H. M. Aktulga, R. Berger, D. S. Bolintineanu, W. M. Brown, P. S. Crozier, P. J. in ’t Veld, A. Kohlmeyer, S. G. Moore, T. D. Nguyen, R. Shan, M. J. Stevens, J. Tranchida, C. Trott and S. J. Plimpton, _Comp. Phys. Comm._ , 2022, 271, 108171.
  * Swope _et al._ (1982) W. C. Swope, H. C. Andersen, P. H. Berens, K. R. Wilson, W. C. Swope, H. C. Andersen, P. H. Berens and K. R. Wilson, _J. Chem. Phys._ , 1982, 76, 637–649.
  * Nosé (2006) S. Nosé, _Mol Phys_ , 2006, 52, 255–268.
  * Parrinello and Rahman (1981) M. Parrinello and A. Rahman, _J. Appl. Phys._ , 1981, 52, 7182–7190.
  * Michaud-Agrawal _et al._ (2011) N. Michaud-Agrawal, E. J. Denning, T. B. Woolf and O. Beckstein, _J Comput Chem_ , 2011, 32, 2319–2327.
  * Richard J. Gowers _et al._ (2016) Richard J. Gowers, Max Linke, Jonathan Barnoud, Tyler J. E. Reddy, Manuel N. Melo, Sean L. Seyler, Jan Domański, David L. Dotson, Sébastien Buchoux, Ian M. Kenney and Oliver Beckstein, Proceedings of the 15th Python in Science Conference, 2016, pp. 98 – 105.
  * Hockney and Eastwood (2021) R. W. Hockney and J. W. Eastwood, _Computer simulation using particles_ , crc Press, 2021.
  * Kolafa and Perram (1992) J. Kolafa and J. W. Perram, _Mol Simul_ , 1992, 9, 351–368.
  * Sayer and Cox (2019) T. Sayer and S. J. Cox, _Phys. Chem. Chem. Phys._ , 2019, 21, 14546–14555.
  * Cox and Sprik (2019) S. J. Cox and M. Sprik, _J. Chem. Phys._ , 2019, 151, 064506.
  * Yeh and Berkowitz (1999) I.-C. Yeh and M. L. Berkowitz, _J. Chem. Phys._ , 1999, 111, 3155–3162.
  * Andersen (1983) H. C. Andersen, _J. Comput. Phys._ , 1983, 52, 24–34.
  * Swendsen and Wang (1986) R. H. Swendsen and J.-S. Wang, _Phys. Rev. Lett._ , 1986, 57, 2607\.
  * Earl and Deem (2005) D. J. Earl and M. W. Deem, _Phys. Chem. Chem. Phys._ , 2005, 7, 3910–3916.
  * Bonomi _et al._ (2019) M. Bonomi, G. Bussi, C. Camilloni, G. A. Tribello, P. Banáš, A. Barducci, M. Bernetti, P. G. Bolhuis, S. Bottaro, D. Branduardi, R. Capelli, P. Carloni, M. Ceriotti, A. Cesari, H. Chen, W. Chen, F. Colizzi, S. De, M. D. L. Pierre, D. Donadio, V. Drobot, B. Ensing, A. L. Ferguson, M. Filizola, J. S. Fraser, H. Fu, P. Gasparotto, F. L. Gervasio, F. Giberti, A. Gil-Ley, T. Giorgino, G. T. Heller, G. M. Hocky, M. Iannuzzi, M. Invernizzi, K. E. Jelfs, A. Jussupow, E. Kirilin, A. Laio, V. Limongelli, K. Lindorff-Larsen, T. Löhr, F. Marinelli, L. Martin-Samos, M. Masetti, R. Meyer, A. Michaelides, C. Molteni, T. Morishita, M. Nava, C. Paissoni, E. Papaleo, M. Parrinello, J. Pfaendtner, P. Piaggi, G. M. Piccini, A. Pietropaolo, F. Pietrucci, S. Pipolo, D. Provasi, D. Quigley, P. Raiteri, S. Raniolo, J. Rydzewski, M. Salvalaglio, G. C. Sosso, V. Spiwok, J. Šponer, D. W. Swenson, P. Tiwary, O. Valsson, M. Vendruscolo, G. A. Voth and A. White, _Nat. Methods_ , 2019, 16, 670–673.
  * Tribello _et al._ (2014) G. A. Tribello, M. Bonomi, D. Branduardi, C. Camilloni and G. Bussi, _Comput Phys Commun_ , 2014, 185, 604–613.
  * Humphrey _et al._ (1996) W. Humphrey, A. Dalke and K. Schulten, _J. Mol. Graph. Model._ , 1996, 14, 33–38.
  * Stone (1998) J. Stone, _MSc thesis_ , Computer Science Department, University of Missouri-Rolla, 1998.
  * Haji-Akbari (2018) A. Haji-Akbari, _J. Chem. Phys._ , 2018, 149, 072303.

Electronic Supporting Information for The limit of macroscopic homogeneous ice
nucleation at the nanoscale

This document includes: methodology for all simulations in the main paper; a
sensitivity analysis of the nucleation rate to $\Delta\mu$; committor
distributions for 2.5 nm and 5 nm films; density profiles for the physically
confined films; volumetric nucleation rates obtained by seeding and cooling
ramps compared to literature values; and preliminary work using seeding
comparing bulk to film samples in TIP4P/ice(8.5).

##  S1 Methodology

All simulations were performed with LAMMPS,Thompson _et al._ (2022) with
particle positions propagated with the velocity verlet algorithm.Swope _et
al._ (1982) The timestep, $\delta t$, was set at $10\,\mathrm{fs}$ in mW
simulations and $2\,\mathrm{fs}$ in TIP4P/ice(8.5) simulations. Simulations of
thin water films employed a slab geometry under periodic boundary conditions,
with the temperature maintained with a Nosé-Hoover thermostat with a damping
constant $100\delta t$. Simulations of bulk water were performed under
periodic conditions in the NPT ensemble; temperature was similarly maintained
with a Nosé-Hoover thermostatNosé (2006), and pressure maintained at 0 bar
with a Parinello-Rahman barostatParrinello and Rahman (1981) with damping
constant $1000\delta t$.

###  S1.1 Preparation of thin mW films

For unsupported thin films, the thickness $W$ was defined as

$W=N/A\bar{\rho}$ (S1)

with $N$ the number of molecules, $A$ the cross sectional area of the film,
and $\bar{\rho}=33.3774\,\mathrm{nm^{-3}}$ the bulk density of mW molecules at
220 K. By running initial tests at 6000 and 8192 molecules in the 6 nm film,
we found that there were no significant finite size effects present for a
critical seed in a 6000 particle system, akin to previous investigations that
selected similar sized lateral dimensionsHaji-Akbari _et al._ (2014); Lü _et
al._ (2013). $W$ was therefore controlled by selecting the $x$ and $y$
dimensions ($l_{x}$ and $l_{y}$ respectively) of the film so that
$l_{x}=l_{y}=\sqrt{N/W\bar{\rho}}=13.40\,\mathrm{nm}\times\sqrt{1\,\mathrm{nm}/W}$.

To generate an initial configuration of mW water molecules for a given film
width, the packmol softwareMartínez _et al._ (2009) was used to randomly pack
6000 mW particles in a cuboid between $(0,\,0,\,-W/2)$ and
$(l_{x},\,l_{y},\,W/2)$. The $z$ dimension ($l_{z}$) of the simulation box was
then increased by a factor of 3 to generate a liquid-vapor interface. The
films were equilibrated at 220 K for $10000\>\mathrm{steps}$ at timesteps of
$0.1\>\mathrm{fs}$, $0.5\>\mathrm{fs}$, $1\>\mathrm{fs}$, $2\>\mathrm{fs}$ and
$5\>\mathrm{fs}$ sequentially, before further equilibration for
$50000\>\mathrm{steps}$ at a $10\>\mathrm{fs}$ timestep. After this, the
systems were judged to have equilibrated, as the density of the center of the
film and enthalpy had converged.

Independent initial configurations of the mW films were generated by
performing simulations of the equilibrated film at $220\,$K and taking the
positions of all atoms every $100\,$ps.

###  S1.2 Preparation of bulk mW systems

An initial random configuration of mW particles was generated using packmol,
with $l_{x}=l_{y}=l_{z}=(N/\bar{\rho})^{1/3}$. The system was then
equilibrated with the protocol described in Sec. S1.1, with the density and
enthalpy adjudged to have equilibrated. Independent configurations were then
sampled as described in Sec. S1.1

###  S1.3 Preparation and simulation of mW seeding systems

Independent configurations of mW as hexagonal ice were generated by first
utilizing the GenIce package to create an initial configuration of TIP4P/ice
molecules in the ice Ih crystal structure. The hydrogens were stripped out
using python code, and the oxygens relabeled as mW particles, to generate a
system of mW particles in the Ih crystal structure with $N=16000$. The system
was allowed to equilibrate at 220 K and 0 bar for 10 ns. Independent
configurations were then obtained by allowing the system to further evolve at
220 K and 0 bar, with the positions of all molecules taken every $100\,$ps.

To generate seeding systems, an initial configuration of mW particles of the
desired system (i.e., thin film or bulk liquid) is selected, and python code
utilizing the MDAnalysis software packageMichaud-Agrawal _et al._ (2011);
Richard J. Gowers _et al._ (2016) is used to cut a spherical cavity in the
liquid. A spherical cluster of Ih is created by cutting a similar sphere out
of an initial configuration of ice Ih. The spherical seed is placed in the
cavity cut in the liquid configuration to initialize a seeding system.

Seeding systems progressed by first allowing the liquid molecules to
equilibrate for 80 ps at 220 K with the molecules in the seed held fixed; for
a seed immersed in the bulk liquid, pressure was also maintained at 0 bar.
This was done to allow a more realistic liquid/crystal interface to develop,
with 80 ps selected to allow enough time for this surface to relax, but not so
long that the seeds close to the critical size grew substantially. Then, all
particles were allowed to evolve until either $n_{\mathrm{cl}}\gg
n_{\mathrm{cl}}^{\ddagger}$ (taken as $n_{\mathrm{cl}}\geq 400$) or
$n_{\mathrm{cl}}\approx 0$ (taken as $n_{\mathrm{cl}}\leq 20$). Information on
how the $n_{\mathrm{cl}}$ was identified and monitored is given in Sec. S2.

###  S1.4 Obtaining n${}_{\rm cl}^{\ddagger}$ and f+ from seeding simulations

The probability of growth, $p_{\rm ice}(n_{\rm cl}^{\ast};W)$ for a given
cluster size $n_{\rm cl}^{\ast}$ is found by considering all initial clusters
of size $n_{\rm cl}^{\ast}-9\leq n_{\rm cl}(t=0)\leq n_{\rm cl}^{\ast}+9$. The
probability of growth for this group of simulations was found as $p_{\rm
ice}(n_{\rm cl}^{\ast};W)=n(\mathrm{Grown})/n(\mathrm{Total})$. The bin size
of 19 was selected as clusters were regularly observed to change size by $\sim
9$ molecules between frames in analysis, and corresponded to a 5% change in
radius between the largest and smallest extremes of the range at criticality.
$p_{\rm ice}(n_{\rm cl}^{\ast};W)$ was identified for all integer $n_{\rm
cl}^{\ast}$, and combined to generate Fig. 3a. $p_{\rm ice}(n_{\rm cl};W)=0.5$
was then identified to find $n_{\rm cl}^{\ddagger}$. The ranges presented as
in Tab. 1. were found by splitting the entire data set (for each system)
randomly into three subsets and repeating the analysis with each subset. The
largest and smallest values obtained were quoted as the range.

Once $n_{\rm cl}^{\ddagger}$ was identified, the attachment rate as found by
identifying all simulations for a given film width of initial cluster size
$n_{\rm cl}^{{\ddagger}}-5\leq n_{\rm cl}(t=0)\leq n_{\rm cl}^{{\ddagger}}+5$,
and initiating simulations from these configurations. $\langle(n_{\rm
cl}(t)-n_{\rm cl}^{{\ddagger}})^{2}\rangle$ was found by assuming $n_{\rm
cl}(t=0)\approx n_{\rm cl}^{{\ddagger}}$ and taking the average from all
simulations. $f^{+}$ was then found by applying a linear fit, and taking half
the gradient.

###  S1.5 Preparation of TIP4P/ice systems

All simulations of TIP4P/ice were performed in a similar manner to their
corresponding mW simulation, with electrostatic interactions computed with a
particle-particle particle-mesh Ewald method,Hockney and Eastwood (2021) 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 Å.Kolafa and Perram (1992) For simulations of liquid water in
contact with its vapor, we set the electric displacement field (along $z$)
$D=0$, using the implementation given in Refs. 87; 88. The rigid geometry of
the molecules was maintained with the RATTLE algorithm Andersen (1983).

TIP4P/ice systems were generated using a similar methodology to mW. The only
film generated was made to match the film considered in Ref. 30, which had
$N=3072$ and $l_{x}=l_{y}=4.8096\,\mathrm{nm}$. Given the bulk density of
TIP4P/ice(8.5) at 230 K was measured in this study as
$\bar{\rho}=30.9366\,\mathrm{nm^{-3}}$, this gives $W=4.3\,\mathrm{nm}$.
Similar to mW, the film was generated by creating an initial random
configuration of TIP4P/ice molecules with packmol, before allowing the system
to equilibrate. The bulk system was also generated to match Ref. 30, with
$N=4096$. It was generated and equilibrated in a manner similar to Sec. S1.2,
but with the relevant fixes described above. In addition, to facilitate
equilibration at this low temperature, we used replica exchange with
temperature 230, 236, 242, 248 and 254 K; swaps between replicas were
attempted every 1 ps.Swendsen and Wang (1986); Earl and Deem (2005)

###  S1.6 Simulations with cooling ramps

Cooling ramp simulations were performed for
$W=\infty,\,6.0\,\mathrm{nm},\,3.5\,\mathrm{nm}$. The initial configurations
generated in Secs. S1.1, S1.2 were taken, and were cooled at a rate of
$0.2\,\mathrm{K\,ns^{-1}}$ using a Nosé-Hoover thermostat with damping
parameter 1000 fs. The systems were monitored every 100 ps, and the system was
adjudged to be frozen when $n_{\mathrm{cl}}\geq 200$. In the rare case that a
system recrossed 200 molecules, the later temperature at which
$n_{\mathrm{cl}}$ passed 200 molecules was taken as the freezing temperature.
Information on how the $n_{\mathrm{cl}}$ was identified and monitored is given
in Sec. S2.

###  S1.7 Preparation of systems with physical confining boundaries

mW films were confined between the (111) face of a Lennard-Jones crystal, with
parameters described in Tab. S1. This surface had been previously investigate
by Davies _et al._ Davies _et al._ (2022), and shown to nucleate ice only for
temperatures $T<201\,\mathrm{K}$, significantly below the temperature of 220 K
used for the seeding investigations performed in this study.

Parameter | Value
---|---
$\epsilon$ | 0.23 $\mathrm{kcal\,mol^{-1}}$
$\sigma$ | 2.5233 Å
LJ Cutoff | 7.53 Å
Lattice Constant | 4.1564 Å
Table S1: Parameters for the LJ crystal used to confine the thin film. This
crystal face is shown to only nucleate ice below 201 KDavies _et al._ (2022),
and so we consider it nucleation inactive for this study, with any nucleation
events occurring at 220 K expected to happen far from the interface

1920 Lennard-Jones atoms were packed into a close-packed layer, and 5 such
layers were combined to generate an FCC crystal with its close-packed face
exposed, perpendicular to the $z$ axis. This resulted in a crystal structure
with dimensions $5.5778\,\mathrm{nm}\times 24.4336\,\mathrm{nm}\times
1.1950\,\mathrm{nm}$. This crystal was placed in a simulation cell with its
uppermost layer at a position $z=-H/2$. An identical structure was placed with
its lowermost layer at $z=H/2$.

Between these two layers, packmol was used to generate a random configuration
of $N(H)$ mW molecules in a cuboid between $(0,\,0,\,-H/2)$ and
$(5.5778\,\mathrm{nm},\,8\,\mathrm{nm},\,H/2)$.
$N(H)=1570\times\frac{H}{[1\,\mathrm{nm}]}$ was selected so that the
physically confined film spanned the $x$ dimension of the LJ crystal, and
approximately $1/3$ of the $y$.

The system was then allowed to equilibrate, and independent initial
configurations sampled, in a manner similar to Sec. S1.1. Seeding systems were
the generated and investigated in a manner similar to Sec. S1.3. The critical
nucleus size was identified in a manner identical to Sec. S1.4.

The particles that make up the LJ crystal were not included in the dynamics,
and did not contribute towards the simulation beyond their interaction with
the mW particles.

##  S2 Order Parameters used for Identifying Cluster Sizes

Clusters of ice-like particles were identified by first classifying each
particle in the system as either liquid-like or crystal like, using the
Steinhardt parameters defined in the main paper:

$\bar{Q}_{6}^{(i)}=\frac{Q_{6}^{(i)}+\sideset{}{{}^{\prime}}{\sum}_{j}Q_{6}^{(j)}}{\nu^{(i)}+1},$
(S2)

with $Q_{6}^{(i)}$ defined with the 6th order spherical harmonics:

$Q_{6}^{(i)}=\frac{1}{\nu^{(i)}}\sqrt{\sum_{m=-6}^{6}\sideset{}{{}^{\prime}}{\sum}_{j,k}Y_{6m}^{\ast}(\hat{\mathbf{r}}_{ij})Y_{6m}(\hat{\mathbf{r}}_{ik})},$
(S3)

For all mW simulations, particle $j$ was considered adjacent (and so included
in the summation of Eqns. S2 and S3) if the distance between particles $i$ and
$j$ $|\textbf{r}_{ij}|\leq 3.5\,$Å. For this we used the implementation in
PlumedBonomi _et al._ (2019); Tribello _et al._ (2014), with the switching
function $\sigma(|\textbf{r}_{ij}|)$:

$\sigma(|\textbf{r}_{ij}|)=\begin{cases}1&\text{if }|\textbf{r}_{ij}|\leq
3.5\,\text{\AA}\\\
\exp\left(-\frac{(|\textbf{r}_{ij}|-3.5\,\text{\AA})^{2}}{2\times(0.05\,\text{\AA})}\right)&\text{if
}3.5\,\text{\AA}\leq|\textbf{r}_{ij}|<3.51\,\text{\AA}\\\ 0&\text{if
}3.51\,\text{\AA}\leq|\textbf{r}_{ij}|\\\ \end{cases}$ (S4)

This switching function leads to the selection of the $\nu^{(i)}$ nearest
atoms to particle $i$, where $\nu^{(i)}$ is defined as

$\nu^{(i)}=\sum_{j\neq i}\sigma(|\textbf{r}_{ij}|)$ (S5)

for all $j$ mW particles in the system.

Once $\bar{Q}_{6}^{(i)}$ is identified for each mW particle $i$ in the system,
particles were identified as being ice-like or liquid-like through the use of
a cutoff value $\bar{Q}_{\rm 6C}=0.4021$. $\bar{Q}_{\rm 6C}$ was found by
obtaining the probability distributions $P(\bar{Q}_{6}^{(i)})$ for both liquid
mW and hexagonal ice mW at 220 K and 0 bar,333Cubic ice was also considered,
and the corresponding distribution was found to be skewed to higher values of
$P(\bar{Q}_{6}^{(i)})$. Since all systems are being seeded purely with
hexagonal ice, we choose to base of derivation of $\bar{Q}_{\rm 6C}$ on the
liquid and Ih phases.  shown in Fig. S1. $\bar{Q}_{\rm 6C}$ was found as the
value for which an equal proportion of Ih and liquid particles are
misidentified.Sanz _et al._ (2013); Espinosa _et al._ (2016)

$\int_{0}^{\bar{Q}_{\rm 6C}}P(\bar{Q}_{\rm 6,\,Hex}^{(i)})d\bar{Q}_{\rm
6,\,Hex}^{(i)}=\int_{\bar{Q}_{\rm 6C}}^{1}P(\bar{Q}_{\rm
6,\,Liq}^{(i)})d\bar{Q}_{\rm 6,\,Liq}^{(i)}$ (S6)

Figure S1: $P(\bar{Q}_{\rm 6}^{(i)})$ for both liquid and hexagonal ice.
Inset: The left hand side (blue) and right hand side (orange) of Eqn. S6 shown
as a function of $P(\bar{Q}_{\rm 6C})$.

Once all particles have been classified as liquid-like or crystal-like, the
largest cluster of adjacent crystal-like particles is identified, using a
depth-first clustering search and identifying particles $i$ and $j$ as
adjacent if $|\textbf{r}_{ij}|\leq 3.5\,\text{\AA}$. The number of molecules
of this cluster is taken as $n_{\mathrm{cl}}$. Performing a similar analysis
for TIP4P/ice at 230 K (with $|\mathbf{r}_{ij}|$ the distance between oxygen
atoms on molecules $i$ and $j$ etc) gives $\bar{Q}_{\rm 6C}=0.3946$.

##  S3 Sensitivity of the nucleation rate from seeding to $\mathbf{\Delta\mu}$

Figure S2: $J(W)$ obtained with (a) $|\Delta\mu|/k_{\rm B}=126$ K, (b)
$|\Delta\mu|/k_{\rm B}=122$ K (same as Fig. 3d in the main text), and (c)
$|\Delta\mu|/k_{\rm B}=118$ K. The results are virtually indistinguishable
from each other.

##  S4 Committor distributions of 2.5 $\mathbf{nm}$ and 5 $\mathbf{nm}$ films

To test the suitability of $n_{\mathrm{cl}}$ as the “reaction coordinate”, we
present the committor distribution for $n_{\mathrm{cl}}$. The committor
distribution was found by first selecting clusters for which
$n_{\mathrm{cl}}(t=0)\approx n_{\mathrm{cl}}^{\ddagger}$ (taken in this case
as $n_{\mathrm{cl}}(t=0)=n_{\mathrm{cl}}^{\ddagger}\pm 3$). This resulted in
72 initial clusters for $W=2.5\,\mathrm{nm}$ and 53 for $W=5\,\mathrm{nm}$.
For each cluster, 20 different simulations were performed where the velocities
of each particle was randomized to create the expected velocity distribution
at 220 K, and allowed to evolve until the seed grew or shrank. The probability
$P(p_{\mathrm{B}})$ that a given cluster would grow to encompass the entire
system is presented below.

Figure S3: Committor distributions for (a) 5 nm and (b) 2.5 nm films.

##  S5 Density profiles for different values of $H$

Figure S4: $\langle\rho(z)\rangle$ for
$H=6.0\,\mathrm{nm},\,H=5.0\,\mathrm{nm},\,H=4.5\,\mathrm{nm},\,H=4.0\,\mathrm{nm}$.
$\bar{\rho}=33.3774\,\mathrm{nm^{-3}}$, the bulk density for mW at 220 K and 0
bar, is included as the dashed line.

##  S6 Volumetric nucleation rates obtained from cooling ramps

Volumetric nucleation rates for the cooling ramps investigation, found by
dividing the freezing rates of the bulk system by the $Nv$ where $N=6000$ is
the number of molecules in the system, and $v=3\times
10^{-29}\,\mathrm{m^{-3}}$ is the volume of a single mW molecule under these
conditions. These results are compared to previous literature values.Haji-
Akbari (2018); Russo _et al._ (2014); Moore and Molinero (2011); Sanchez-
Burgos _et al._ (2022)

Figure S5: $J(\infty)$ for this work and a range of previous studies.

##  S7 Preliminary Work on TIP4P/ice Films

We argue in the main paper that mW is suitable for this study, due to the
mechanisms of ice formation being similar to TIP4P/ice(8.5), and similar
convergence of structural order to bulk-like values within $\sim 1\,$nm from
the surface. In this section, we present preliminary work using the
TIP4P/ice(8.5) model for both bulk and film ($W=4.3\,\mathrm{nm}$). Seeding
was implemented in a manner similar to Sec. S1.3, but at 230 K and a 200 ps
equilibration time instead of 80 ps. Additionally, the same method was used to
identify ice-like molecules, with the oxygen in each molecule acting as the
point which $\textbf{r}_{ij}$ vectors point to and from. Due to the different
conditions and model, $\bar{Q}_{\rm 6C}=0.3946$.

Figure S6: $n_{\mathrm{cl}}(t)$ for 8 independent seeding simulations in a) a
bulk system, $\langle n_{\mathrm{cl}}(t=0)\rangle=232$ and b) a 4 nm film
$\langle n_{\mathrm{cl}}(t=0)\rangle=200$. In both cases, clusters of similar
initial size grow and shrink, and some clusters spend significant time at or
around their initial sizes. Both facts indicate these systems are at or near
criticality.

This investigation was extremely limited in scope due to the high
computational cost of TIP4P/ice(8.5). However, we show that both film and bulk
appear to be near critical for $n_{\mathrm{cl}}(t=0)\approx 200$. The fact
that bulk and film systems appear to report similar critical nuclei further
supports our mW-based study.