Solid – liquid interfacial free energy of ice Ih , ice Ic , and ice 0 within a mono-atomic model of water via the capillary wave method

We apply the capillary wave method, based on measurements of fluctuations in a ribbon-like interfacial geometry, to determine the solid-liquid interfacial free energy for both polytypes of ice I and the recently proposed ice 0 within a mono-atomic model of water. We discuss various choices for the molecular order parameter, which distinguishes solid from liquid, and demonstrate the influence of this choice on the interfacial stiffness. We quantify the influence of discretisation error when sampling the interfacial profile and the limits on accuracy imposed by the assumption of quasi one-dimensional geometry. The interfacial free energies of the two ice I polytypes are indistinguishable to within achievable statistical error and the small ambiguity which arises from the choice of order parameter. In the case of ice 0, we find that the large surface unit cell for low index interfaces constrains the width of the interfacial ribbon such that the accuracy of results is reduced. Nevertheless, we establish that the interfacial free energy of ice 0 at its melting temperature is similar to that of ice I under the same conditions. The rationality of a core-shell model for the nucleation of ice I within ice 0 is questioned within the context of our results.


I. INTRODUCTION
The crystallisation of ice I from liquid water continues to be a benchmark problem for molecular simulation.][3][4] Homogeneous nucleation rates have been computed using the framework of classical nucleation theory [5][6][7] (CNT) or directly with path sampling approaches. 8,9][12][13] One particular subtlety presents a challenge to the accuracy of molecular simulations.As yet no simulation has demonstrated a propensity for ice I to crystallise entirely in its hexagonal polytype ice Ih (ABABAB stacking of molecular bilayers) at weak supercooling.Instead, simulations consistently produce stacking disordered mixtures of the cubic polytype ice Ic (ABCABC stacking) and ice Ih, as seen experimentally only at very low temperatures.The origin of this stacking disorder, and the preference for ordered hexagonal stacking in ice formed at weak supercooling, has been the subject of much discussion in the literature.For details, we refer the reader to a recent review by Malkin et al. 14 In a recent communication, one of us (DQ) presented a thermodynamic argument that low temperature critical nuclei of ice I are most stable with stacking disordered structure. 15his model is based on nearest neighbour interactions between ice bilayers.In the context of nucleation from the liquid, rather than transformation from other ice phases, experiments a) Current address: School of Mathematics and Physics, University of Lincoln, Lincolnshire LN6 7TS, United Kingdom.b) d.quigley@warwick.ac.uk suggest that the stacking sequence is entirely random (no correlations involving more than three bilayers) and hence the nearest-neighbour approximation is reasonable. 16At this level of approximation, any difference in solid-liquid interfacial free energy between cubic and hexagonal ice is neglected.
Historical experimental crystallisation data at low temperature can be interpreted as contradicting this assumption, based on CNT estimates for the interfacial free energy γ c between cubic ice and supercooled water. 17Calculations of γ c from this data require a figure for the cubic-hexagonal free energy difference ∆G ch .Calculations of this quantity are subject to considerable variation, ranging from a few J mol 115 in the popular coarse-grained mW model 18 to 627 ± 300 J mol 1 in first principles calculations with lattice vibrations treated anharmonically. 19For any positive estimate of ∆G ch , a value for γ c , lower than the interfacial free energy between hexagonal ice and water γ h , is needed to explain the preferential nucleation of cubic ice within CNT. 20 Despite the likelihood that these historical experiments yielded stacking-disordered rather than pure cubic ice, the notion that variations in stacking sequence influence interfacial free energy is worthy of investigation.
2][23][24] This can require simulations on large length and time scales but is increasingly tractable.The method is implemented entirely via postprocessing of simulation trajectories and does not require invasive modifications to simulation packages.Notably, such calculations have been largely confined to crystals with cubic symmetry; however Benet, MacDowell, and Sanz 25 have recently computed γ h for the popular TIP4P/2005 26 atomistic model of water via this approach.
In the present work, our primary objective is to compute the solid-liquid interfacial free energy for both cubic and hexagonal ice in the mono-atomic mW model of water 18 and to explore the factors which control and limit accuracy of the capillary wave method when applied to ice-water interfaces.A recent study 27 using the "mold integration technique" provides a useful benchmark for comparison.This was unable to resolve differences in interfacial free energy between facets of ice Ih in the mW model or between the orientationaveraged γ c and γ h when using the TIP4P/Ice 28 model.No comparison between the two polytypes was made for the mW model.
The mW model represents the water molecule with a Stillinger-Weber type potential, favouring tetrahedral local environments and parametrised to capture structural and thermodynamic properties of the ice-water system.The model has been used in a number of ice nucleation simulations, in all cases leading to crystallites of mixed cubic and hexagonal character. 3,8,29,30This model represents an interesting case, as the difference in bulk stability between cubic and hexagonal ice is very small. 15This leads to the expectation of stacking disordered nuclei at all temperatures where homogeneous nucleation is achievable, should the assumption that γ h ≈ γ c hold.
In a recent work, Russo, Romano, and Tanaka 31 have obtained an improved CNT fit to low-temperature nucleation simulations (also using the mW model) by including a low density metastable ice phase (ice 0) as the first nucleated species within a 3 parameter core-shell nucleation theory.The interfacial free energy γ z between ice 0 and water is wrapped into one of these fitting parameters and must be somewhat lower than both γ h and γ c to offset the lower bulk stability of ice 0. 32 As a secondary objective, we compute the interfacial free energy between ice 0 and water directly and determine how this influences the interpretation of the core-shell fit in Ref. 31.

II. METHODS
2][23][24] The interfacial free energy is dependent on the crystallographic surface exposed to the melt and can be represented as a function on the unit sphere, Here the functions S i represent a basis set of symmetryadapted spherical harmonics appropriate to the point group of the solid crystal.Provided the free energy is only weakly anisotropic, the symmetry adapted nature of these functions ensures that the series can be truncated after the first few terms.In the case of crystals with cubic symmetry, the appropriate basis functions are the Kubic harmonics.In the general case, generation of appropriate functions in a suitable form is non-trivial.Implementation of an automated method for doing so, with only the crystallographic unit cell as input, 34 is available via the "gencs" utility distributed with the Alloy Theoretic Automated Toolkit (ATAT). 35This enables application of the capillary wave method to a much wider range of systems.Appropriate functions S i for cubic and hexagonal ice I, and for ice 0, are given in the supplementary material as Equations ( S3)-(S5).
As the interfacial free energy cannot be determined directly from simulations, the coefficients i are obtained via fitting to multiple measurements of the interfacial stiffness γ along directions in the tangent plane perpendicular to the surface normal (see Figure 1).The stiffness along a particular direction û in the tangent plane is related to the interfacial free energy via 33 γ (θ, where H γ (θ, ϕ) is the Hessian matrix (in polar coordinates) expressing the second derivatives of γ (θ, ϕ) evaluated at angles θ and ϕ.
The stiffness itself is computed via simulation of a "ribbon-like" interfacial geometry along û, under conditions of equilibrium coexistence.With reference to the example geometry in Figure 1, the interface is divided into N z columns of width δz, centred on discrete values z n , containing an element of solid-liquid interface of area δA = l y δz.In each column the quantity h (z n ) is computed, indicating the x position of the interface relative to its mean position.
Following earlier authors, 24 we denote the simulation setup via the Miller indices of the (hkl) crystallographic plane exposed to the liquid and the direction [h k l ] along which the simulation box is shortest (the y direction in Figure 1).For example, a simulation geometry for the (110) surface of cubic ice in which the short direction is taken as [001] will be denoted as (110) [001].
When comparing solid-liquid interfacial free energies between two crystal structures, it is important to calculate FIG. 1. Simulation geometry.Cartesian vectors i, j, and k denote the crystal coordinate system such that the polar angles θ and ϕ define a surface normal.Molecular dynamics simulations are performed in a local coordinate system with the x direction chosen along this normal.In this case, the direction û = e θ in the tangent plane is chosen as the z direction.
interfacial height profiles with the same choice of local order parameter ξ, used to locate the x position of the interface and hence h (z n ) in each column.For the cubic ice Ic structure, the commonly used q 6 order parameter performs well but recognises only three solid neighbours in ice Ih leading to an erroneous large difference between the interfacial free energies of these two polytypes.
A better choice, which is appropriate for ice I interfaces only, is the CHILL+ algorithm of Nguyen and Molinero. 36his is based on the ubiquitous Steinhardt order parameter. 37 per-molecule vector q 3 is constructed with 7 components, each of which is a sum of the l = 3 spherical harmonics calculated for the unit vectors connecting a molecule to its four nearest neighbours.The correlation of q 3 between each molecule and its neighbours is used to classify molecules as solid, liquid, or interfacial.We refer readers to Ref. 36 for details and note that a very similar order parameter was used to construct a reaction coordinate for ice nucleation in earlier work. 8In our work, we assign ξ = 10 to molecules identified as ice-like by CHILL+, ξ = 0 to molecules identified as liquid, and ξ = 5 to molecules (if any) identified as interfacial.We then plot (as a function of x) the quantity ξ as the average value of ξ over all molecules within each of the N x elements of width δx along the column perpendicular to the area element δA.A tanh function is then fitted to these data and the x position of the interface is taken as the inflexion point.
To assess any variability of results arising from the choice of local order parameters, we also compute the interfacial position using the more complex order parameter of Russo, Romano, and Tanaka. 31This order parameter is also able to distinguish ice 0 from liquid.Here the vector order parameter q 12 is used, accumulated over the 16 nearest neighbours of each molecule.The vector q 12 for each molecule is then replaced with the average over these 16 neighbours.Molecules are deemed to be connected if the dot product of their respective normalised q 12 vectors exceeds a threshold value of 0.75.Here we set ξ equal to the number of such connections on each molecule and locate the interface by fitting ξ as a function of x to a tanh function as above.Example data for both methods of locating the interface are shown in Figure 2.An example interfacial height profile is shown in Figure 3.   10] geometry.The computed interface shape is plotted as red points, each of which indicates the interfacial position at each of 50 discretised points in the z direction.Water molecules are coloured according to their number of solid-like connections, in this case identified via the q 12 order parameter.Particles with more than 12 connections are coloured grey, with the remainder coloured turquoise.Red beads mark the position of the interface (also identified by q 12 ) in each interfacial column.System dimensions are those reported in Table I.
It is apparent that the two order parameters identify interfacial positions which can differ by an amount comparable to the intermolecular spacing in the crystal.Should this difference prove to be anything but a constant offset, then the shape of h(z n ) (and hence the calculation of the interfacial free energy) will be order-parameter dependent.We investigate this explicitly in Section IV.
For small wavenumbers q, the discrete Fourier transform h (q n ) of these profiles can be connected to the interfacial stiffness via where A = l y l z is the interfacial area.By plotting measurements of this quantity against q n on a double logarithmic scale, one obtains a linear plot from which the interfacial stiffness can be extracted from the intercept.The requirement of small wavenumber restricts linear behaviour to small values of q n δz such that sin (q n δz) ≈ qδz. 22In our work, we restrict the linear fit to modes for which q n δz < 0.5.We also note that Equation ( 3) is dependent on the interfacial ribbon having a width sufficiently small that no energy is present in modes along the short or y direction in Figure 1.For facets with large surface unit cells, that requirement may be difficult to reconcile with that of surface periodicity, leading to a spurious dependence of the measured stiffness on the ribbon width l y .We calculate γ for two mutually orthogonal choices of the surface direction û at each of the several surface planes, averaging the resulting stiffness when these are equivalent by symmetry.The coefficients (γ 0 , 1 , 2 . . . ) are then chosen to minimise the total difference between the measured stiffness values and those obtained by substituting Equation ( 1) into (2).For further details of the method we refer the reader to Ref. 21.
All molecular dynamics simulations were performed using the LAMMPS 38 package.Interfacial configurations were sampled from 20 ns trajectories following 5 ns careful equilibration in the isothermal-isobaric ensemble 39,40 under atmospheric pressure.h(q n ) was sampled at intervals of 10 ps during subsequent simulations in the NPH ensemble, in which the enthalpy is conserved to within fluctuations of the barostat kinetic energy.This has the advantage of constraining the longterm average solid fraction via a feedback mechanism between latent heat and temperature.The three cell dimensions l x , l y , and l z were allowed to vary independently subject to constant 90 • cell angles.This ensures that the interface is unstrained at all points of the feedback cycle.All simulations used a time step of 1 fs.Simulations of ice I interfaces were conducted at the mW melting temperature using 30 000-50 000 molecules with an interfacial ribbon of length 180-200 Å. Ice 0 simulations required substantially larger ribbon lengths to capture sufficient wavevectors for accurate fitting.In this case simulations used 80 000-180 000 molecules at the ice 0 melting temperature of 240 K, with ribbon lengths of 570-620 Å.

III. GEOMETRY AND CONVERGENCE
Table I reports the simulation geometries used in reconstructing the interfacial free energy of all three ice structures.In this section, we carefully consider the convergence of the measured interfacial stiffness for example surfaces.We examine how the spatial resolution at which the interfacial ribbon is sampled and the simulation run length govern accuracy of stiffness measurements.We also establish how the width of the ribbon (y direction in Figure 1) influences stiffness measurements.

A. Spatial resolution
We first consider convergence with respect to the spatial sampling resolution δz along the interfacial ribbon.It is important to note two distinct effects of this parameter.Smaller δz leads to a more accurate representation of the interface geometry but also affords a larger maximum wavevector to be included in fitting Equation (3). Figure 4 plots the stiffness measured for a selection of interfacial geometries as a function of N z when including all wavevectors q for which qδz is less than 0.5.Jumps in stiffness are observed at values of N z for which new wavevectors commensurate with the periodicity of the simulation box become available.Convergence to within statistical uncertainty is achieved at N z = 50 for both ice I polytypes, and at N z = 100 for ice 0, which for both cases corresponds to a δz close to interatomic spacing.We note that in the case of ice 0, sampling with N z = 20 leads to a stiffness less than half the converged value.
Figure 5 plots the measured stiffness versus N z for the ice Ih (basal) [11 20] geometry at a fixed number of wavevectors.Comparison with the central panel of Figure 4 indicates that enhanced resolution in the representation of h(z n ), rather than inclusion of additional wavevectors in the analysis, has the greater effect on convergence.
Figure 5 also shows the effect of spatial resolution in the direction normal to the interface, demonstrating that a sufficiently high N x must be used to accurately capture the FIG. 5. Interfacial stiffness as a function of N z at fixed N x = 80 (left) and as a function of N x at fixed N z = 50 (right).In both cases, the simulation geometry is ice Ih (basal) [11 20].Simulation dimensions are those reported in Table I.Red points indicate calculations using height profiles computed via the CHILL+ algorithm and black via the q 12 solidity criterion.In all cases, γ is computed from only the 8 shortest wavelengths which satisfy qδz < 0.5.
interfacial height h(z n ) and hence the stiffness.Errors introduced by insufficient resolution along z and x are similar in magnitude but of opposite sign.Sampling at very close to the ultimate limit of spatial resolution (the interatomic spacing) is required to achieve converged results.We adopt the strategy of sampling at as close to this limit as possible for the remainder of this paper.

B. Thickness of interfacial ribbon
The standard analysis which leads to Equation (3) assumes a pseudo-one-dimensional geometry in which the width of the interfacial ribbon l y is assumed to be sufficiently thin that no energy is present in fluctuations along y.For the lowest index surfaces of crystals with cubic or orthorhombic symmetry, very thin interfacial ribbons can easily be constructed in the geometry of Figure 1 for which l y is commensurate with the periodicity of the surface unit cell.For higher index surfaces, or for crystals (such as ice Ih) with non-orthogonal symmetry, undesirable (for purposes of analysis) skewed simulation cells are required to achieve the minimum possible ribbon thickness.In practice, one normally accepts a thicker ribbon to achieve orthogonal supercells and to satisfy the minimum image convention for efficiency.The effect of using these larger ribbon thicknesses has been explored previously in the context of aluminium 22 but never (to our knowledge) for ice.
To establish a benchmark property of a flat solid-liquid interface, we simulate the ice Ih (basal)[prism] geometry using a reduced cross section of l y = 13.3Å, l z = 15.4Å, over a sufficiently short period of time that capillary waves are not observed.The root mean square deviation (RMSD) of the interfacial height profile in this simulation is measured as 0.69(8) Å, which we take as the reference value for an atomically "flat" interface, i.e., the minimum possible deviation from a constant interfacial height resulting from thermal fluctuations of the atomic position about their lattice sites only.
In Figure 6 we plot the convergence of the interfacial stiffness as a function of N z for three thicknesses of the interfacial ribbon in the ice 1h (basal)[prism] geometry.In each case, we FIG. 6. interfacial stiffness as a function of N z for three thicknesses of the interfacial ribbon in the ice Ih (basal) [11 20] geometry.Other dimensions are as reported in Table I.The RMS deviation from a perfectly uniform profile across the interfacial ribbon (y direction in Figure 1) is also indicated.The RMSD for a "flat" surface is 0.69 (8), indicating that for l y = 13.3Å no capillary waves are present in the y direction.The interface position is located via the q 12 solidity criterion.
further divide each discrete interfacial element into sections of width δy and measure the RMSD from the mean height profile over these sections.The resulting quantity can be compared to the value of 0.69(8) Å to establish if the geometry is sufficiently narrow to suppress capillary waves across the ribbon.An appreciably different stiffness is observed for the narrowest ribbon, which is the only geometry in which the above criterion is satisfied.For thicker ribbons, the presence of energy in modes across the ribbon results in an overestimate of the pseudo-one-dimensional stiffness.
As the stiffness of ice I solid-liquid interfaces is expected to be approximately isotropic, we adopt the strategy of keeping the interfacial ribbon at this thickness or smaller, while remaining sufficiently wide to avoid the self-interaction of molecules with their periodic images.For ice 0, this is problematic.The tetragonal unit cell of ice 0 results in larger minimum ribbon thicknesses even for low index interfaces due to periodicity constraints.Furthermore, as the interfacial stiffness of this phase is expected to be lower, one expects the criteria for the suppression of capillary waves across the ribbon to impose a narrower geometry than that for ice I.We must therefore acknowledge that for the ice 0 system sizes reported in Table I, the resulting stiffness and free energy estimates represent an upper bound.

C. Simulation length
Each measurement of stiffness should be made over a sufficient number of fluctuations to accurately represent the ensemble of interfacial height profiles in the current geometry.In Figure 7 we plot the convergence of the interfacial stiffness as a function of time for representative geometries of ice Ic, ice Ih, and ice 0. In the case of ice I phases, the interfacial stiffness measurement is converged to within statistical error after approximately 80 ns of simulation.In the case of ice 0, where the use of larger l z leads to additional spatial averaging, convergence is achieved at approximately 25 ns.For all interfacial geometries, we carefully ensure that the stiffness measurement is converged with respect to the simulation length to a similar accuracy.FIG. 7. Calculated interfacial stiffness as a function of simulation time.The left panel plots the stiffness for the ice Ih (basal) [11 20] (black) and the ice Ic (111)[1 10] (red) geometries which expose identical crystal layers to the liquid.The right panel shows much more rapid convergence for the ice 0 (001)[010] system.The interface position is located via the q 12 solidity criterion in all cases.

IV. RESULTS
We simulate ice-water interfaces for a number of surfaces corresponding to large interplanar spacings D hkl , these being expected to have the lowest interfacial free energy and hence be most commonly expressed in the crystal morphology. 41

A. Ice Ih
Plots from which the stiffness of ice Ih interfaces were extracted are shown in Figure 8 for the case where the interface was located using the q 12 solidity criterion.All plots should share a common gradient of −2 on the logarithmic scales shown, which is applied as a constraint in the fitting procedure.
When fitting these data to second derivatives of Equation ( 1), we find that the contribution from the term in § 2 is negligible.Any improvement upon including either S 3 and/or S 4 is marginal.The best fit is achieved when including the isotropic term plus S 1 and S 4 only, constraining 2 and 3 to zero.This suggests that the weak deviations from an isotropic interfacial free energy are almost entirely captured via the S 1 term.Details of the optimised parameters and the fitted versus measured stiffnesses are reported in Tables S1, S4, and S7 of the supplementary material.FIG. 8. Logarithmic plots of the average Fourier-transformed interface height squared weighted by the area A divided by thermal energy k B T as a function of q for various simulations of interfaces between ice Ih and water.Solid lines are a linear fit to the data points satisfying qδz < 0.5.In this case, the q 12 solidity criterion is used to locate the interface.The final estimates of interfacial free energy for ice Ih at coexistence are given in Table II, where the results of Espinosa, Vega, and Sanz 27 computed via the mold integration technique are shown for comparison.

B. Ice Ic
Plots from which the interfacial stiffnesses of ice Ic were extracted are shown in Figure 9.
FIG. 9. Logarithmic plots of the average Fourier-transformed interface height squared weighted by the area A divided by thermal energy k B T as a function of q for various simulations of interfaces between ice Ic and water.Solid lines are a linear fit to the data points satisfying qδz < 0.5.In this case, the q 12 solidity criterion is used to locate the interface.When fitting these data to second derivatives of Equation ( 1), we find a systematic improvement in the fit when including terms up to and including S 3 .Including the S 4 term leads to no discernible change to within the statistical uncertainty of the resulting interfacial free energies.Calculated interfacial free energies for the four ice Ic interfaces studied are presented in Table III.
Details of the optimised parameters and the fitted versus measured stiffnesses are reported in Tables S2, S5, and S8 of the supplementary material.
As with ice Ih, interfacial free energies calculated via the use of the q 12 solidity criterion are systematically higher than those obtained via CHILL+.There is no significant difference in the overall magnitude of the free energy between the two polytypes.The (111) plane of ice Ic and the basal plane of ice Ih expose identical surface structures to liquid up to the second sub-surface ice double-layer.Within our calculations, the free energies of these two surfaces are identical to within statistical error for the CHILL+ order parameter but differ by a small but measurable amount when using q 12 .

C. Ice 0
Plots from which the interfacial stiffnesses of ice 0 were extracted are shown in Figure 10.Note that results are calculated using only the q 12 method of identifying the solid-liquid interface, as CHILL+ does not usefully recognise ice 0 as a solid.
Fitting the measurement stiffnesses to second derivatives of Equation ( 1) retained terms up to and including S 5 .Details of the optimised parameters and the fitted versus measured stiffnesses are reported in Tables S3, S6, and S9 of the supplementary material.Fitted interfacial free energies for ice 0 are presented in Table IV.

A. Ice Ih
We first discuss our results for the interfacial free energy of ice Ih in comparison to the mold integration calculations of Espinosa, Vega, and Sanz 27 for the mW model.Both our methods of identifying the solid-liquid interface consistently lead to a small but measurable anisotropy in estimates of γ h favouring the basal plane, a feature not evident in the mold integration results.While the statistical uncertainty in our estimates is smaller, there is a systematic difference between our two sets of data, the q 12 order parameter leading to slightly higher interfacial free energies at all surfaces.We note that the mold integration method does not require an order parameter for identifying molecules which lie within a solid-like environment, which from these results would appear to be a source of small but significant systematic error in capillary wave calculations.
The fitted isotropic component γ 0 = 34.5(2)mJ m 2 for the CHILL+ solid identification and 36.0(3)mJ m 2 when FIG. 10.Logarithmic plots of the average Fourier-transformed interface height squared weighted by the area A divided by thermal energy k B T as a function of q for various simulations of interfaces between ice 0 and water.Solid lines are a linear fit to the data points satisfying qδz < 0.5.The q 12 solidity criterion is used to locate the interface.using q 12 .These bracket the value of 35.5(3) mJ m 2 obtained by Espinosa et al. 7 for the mW model when extrapolating results of the seeding method to coexistence.Despite these small systematic differences from earlier work, our results confirm that the ice Ih-water interfacial free energy for the mW model lies toward the upper end of the experimental estimates and is somewhat higher than values obtained with detailed atomistic models. 7,25

B. Ice Ic
We now compare values of interfacial free energy between the two ice I polytypes.When using the CHILL+ method of identifying the solid-liquid interface, the (111) face of ice Ic and the basal plane of ice Ih have identical interfacial free energies to within the statistical uncertainty of our measurements and fits.It might be expected that these two surfaces are both perpendicular to the stacking direction and hence present identical molecular arrangements to the liquid within the first two bilayers.Hudait et al. 30 have recently demonstrated by other methods that the surface tension of these two surfaces is identical.In our work, the q 12 method of identifying the solidliquid interface resolves a slight difference in interfacial free energy between these two surfaces; however, this difference is smaller than the discrepancy between the two choices of order parameter and cannot be considered significant.We also note that the isotropic component of the interfacial free energy γ 0 is identical to within statistical uncertainty between the two polytypes when using q 12 (Tables S4 and S5 of the supplementary material), but there is a small systematic difference in this quantity when using CHILL+.Again this is smaller than the difference in interfacial free energy measured by the two choices of order parameter.
We therefore conclude, to within the error inherent in the choice of order parameter and finite sampling of the interfacial fluctuations, that the interfacial free energies of cubic and hexagonal ice I should be considered identical within the mW model.We have not considered the interfacial free energy of surfaces formed by stacking disordered ice but would expect this to be indistinguishable from the pure polytypes given the above results.

C. Temperature dependence of γ for ice I
The capillary wave method is limited to measurements of interfacial free energy at the melting temperature.Here the area of an infinite plane interface is well defined by the geometry of the simulation cell.As recently discussed by Cheng, Tribello, and Ceriotti, 42 away from coexistence the partitioning of nucleus free energy into volume and surface area terms is essentially arbitrary.If one chooses (as in the seeding method) to calculate the volume contribution as the number of molecules n in the solid nucleus multiplied by the bulk chemical potential difference between solid and liquid, then the resulting γ will depend on the interfacial area and hence how one chooses to define the (usually spherical) dividing surface.To be consistent with the usual assumptions of CNT, this surface should be chosen such that the enclosed volume is V = ρ s n, where ρ s is the number density of molecules in the bulk solid.Cheng, Tribello, and Ceriotti 42 further demonstrate that imposition of this criterion is equivalent to the Tolman curvature correction to γ.
We note here that the seeding method calculations of Espinosa et al. 7 are necessarily consistent with CNT by design and use only data at the critical nucleus size.The resulting fitted interfacial free energies are Tolman-corrected for the curvature of the critical nucleus, the size of which is temperature dependent.In Figure 11 we compare these data to interfacial free energies extrapolated from our measurements via the Turnbull correlation, 44,45 where ∆H(T ) is the temperature dependent enthalpy difference between solid and liquid and C is a constant determined by our measurement of γ 0 at coexistence.We use measurements of ∆H (eV/molecule) and ρ s (molecules per cubic angstrom) readily obtained via straightforward simulation.Here we use the isotropic component of γ h (in mJ m 2 ) computed using the q 12 order parameter for comparison to ice 0 in Sec.V D. In these units, this leads to a Turnbull coefficient of C = 38.5.In all cases, the fitted γ 0 (T ) obtained via the seeding method calculations of Espinosa et al. 7 lies below that estimated from the Turnbull correlation, consistent with a positive Tolman length and a lowering of the interfacial free energy as expected for ice nuclei. 46Furthermore, the magnitude of the difference increases with decreasing temperature, i.e., in inverse proportion to the critical nucleus radius of curvature.We conclude that the Turnbull correlation is useful in extrapolating interfacial free energies away from coexistence, provided that the resulting data are interpreted as lacking the Tolman correction.This is consistent with the work of Limmer and Chandler 43 who compute planar interfacial free energies away from FIG. 11.Isotropic component of the solid-liquid interfacial free energy for ice Ih and ice 0 as a function of the temperature below freezing of ice I. Temperature dependence is extrapolated from calculations at coexistence (solid dots) via the Turnbull correlation.The values fitted to simulations of seeded ice Ih critical nuclei by Espinosa et al. 7 are shown for comparison, as are calculations of ice I planar interfacial free energies by Limmer and Chandler. 43The dashed line indicates the temperature at which Russo, Romano, and Tanaka 31 have inferred nucleation via a core-shell model involving ice 0.
coexistence and favourably compare these to predictions from the Turnbull correlation.
It should be noted that improved fits to explicitly computed nucleation barriers have been reported when allowing the solid-liquid chemical potential difference to deviate from the bulk value, both for mW ice 3 and simple lattice models, 47 i.e., a value of γ 0 consistent with the usual assumptions of CNT may not be optimal.

D. Core-shell model for nucleation of ice 0
We now turn our attention to the results for ice 0 interfaces.Our results indicate an anisotropic interfacial free energy which would lead to non-spherical crystallites slightly compressed along the c-axis.Espinosa et al. 48have recently calculated a value for γ z = 35.5 mJ m 2 at the ice 0 melting temperature using mold integration.Details of the exact surface used are not given; however, this figure is remarkably consistent with our calculation at the (100) surface.
As calculated in Ref. 32, the value of γ z need only be 10% lower than γ h in order for the barrier to ice 0 nucleation to be smaller than for ice I. Russo, Romano, and Tanaka 31 interpret this preferential nucleation via a modified CNT in which nuclei have an internal ice I core structure coated in a shell, thickness ∆R, of interfacial ice 0. For small widths of the outer shell, this reduces to the following form of free energy ∆F of a growing nucleus relative to the metastable liquid, where n is the number of particles within the growing ice I nucleus and β is the inverse temperature.The quantity γ z appears in both b and c, which were treated as fitting parameters in Ref. 31.By repeating the fit of Equation ( 5 Within this interpretation, one can reasonably neglect the Tolman correction when comparing γ h to γ z as the radii of curvature for the ice 0 shell and the competing pure ice I nucleus are similar.Figure 11 includes the temperature dependence of the isotropic component γ 0 of the interfacial free energy between ice 0 and liquid water, extrapolated via the Turnbull correlation as for ice I.For ice 0, the Turnbull coefficient using the same units as above is C = 312, very much larger than for ice I.This is due to the much smaller difference in enthalpy between solid and liquid at coexistence.Consequently, the gradient of γ z (T ) is much steeper than for ice I.At 215.2 K it is clearly energetically favourable to form a liquid interface to ice 0 rather than to ice I.
Using the data from Figure 11 and other thermodynamic parameters as used in Refs.15 and 32, the CNT barrier heights to the nucleation of ice 0 and ice I are very similar (within a few percent) at 215.2 K, suggesting ice 0 nucleation is at worst an energetically competitive crystallisation route in the mW model.Explicit calculations of γ away from coexistence would be required to draw definitive conclusions which do not rely on the empirical Turnbull correlation.
It is however informative to revisit the above fit to a coreshell nucleation model using our extrapolated value for γ z at 215.2 K.The thickness of the ice 0 shell is related to the fitting parameter c by the relationship where ρ s is the number density of molecules in the ice I core.
Constraining the interfacial free energy between ice 0 and water to that in Figure 11 requires that the shell thickness ∆R is 3.3 Å.Only slightly larger values would result from Tolmancorrected values of γ z .We must acknowledge that our ice 0 results should be interpreted as an upper-limit on γ z due to the use of large interfacial ribbon widths, more accurate results may increase the calculated shell thickness further.Based on the trend for ice I in Figure 6, we suggest that this could not increase the shell thickness by a factor of two, required to accommodate the smallest unit cell dimension of ice 0. We suggest that the improved fit to nucleation data obtained in Ref. 31 is a result of the extra degree of freedom in the fitting procedure rather than an indication of core-shell behaviour.
In addition to their mold integration calculations at the mW ice 0 melting temperature, Espinosa et al. 48have also used simulations seeded with spherical nuclei to calculate the nucleation rate of pure ice 0 crystallites at a higher temperature of 226.5 K.We note that the use of spherical crystallites in that work may under-represent the lower energy (001) interface of ice 0. At this temperature, the isotropic estimate of γ z was calculated as 29.9 mJ m 2 , leading to a nucleation rate 200 orders of magnitude lower than that for ice Ih under the same conditions.This value of γ z is indeed higher than that predicted by our calculations, as would be expected for spherical nuclei, but (as discussed above) implicitly contains a Tolman correction for the critical nucleus size which would be expected to act against such a difference.We cannot conclusively state that the use of non-spherical nuclei in seeding calculations would lead to a higher nucleation rate for ice 0.
Taken in isolation, our results do not rule out the nucleation of pure ice 0 crystallites at very low temperatures in the mW model; however, we rely on use of the Turnbull correlation which has not been validated by other methods for ice 0. In any case, the lower stability of ice 0 in more detailed models suggests that the nucleation of this phase is not realistic. 32

VI. CONCLUSIONS
We have computed the solid-liquid interfacial free energy for both cubic and hexagonal polytypes of ice I via the capillary wave method, producing results which are in broad agreement with earlier work.The free energy is only weakly anisotropic.We have explicitly explored the convergence of interfacial stiffness measurements with respect to simulation geometry, discretisation, and simulation length, and compared results for ice I polytypes when using two different methods of assigning molecules to the solid and liquid phases.We find a small systematic difference between results calculated by these two order parameters, which is only resolvable due to careful minimisation of the statistical error.
We have also computed the interfacial free energy of ice 0, γ z , by the same method, with the caveat that interfacial ribbonwidths are constrained to be rather large by the surface unit cells of low index terminations.This may lead to an overestimate of γ z , implying that our calculations represent an upper limit on this quantity.One option for removing this limitation in future work would be to adopt a two-dimensional analysis of capillary waves. 49,50Here the amplitude of the observed waves is logarithmic rather than linear in the dimensions of the interfacial cross section, requiring larger simulations than those reported here before the level of statistical accuracy becomes comparable to the systematic error introduced by the finite thickness.
A limitation of the capillary wave method restricts our knowledge of these interfacial free energies to the melting temperature of each phase.However, by exploiting the empirical Turnbull correlation, we can estimate γ h and γ z at temperatures where the nucleation of ice via a core-shell mechanism involving ice 0 has been suggested.We show that the Turnbull correlation is consistent with previous measurements of the ice I interfacial free energy provided the resulting quantities are interpreted for flat interfaces, i.e., lacking a Tolman correction.
Comparison of γ h and γ z suggests that the preferential nucleation of ice 0 may be possible at low temperature.However our calculated values combined with existing thermodynamic data suggest that a core-shell nucleation process is not a viable description, as the resulting shell thickness is comparable to a single molecular diameter and significantly smaller than the ice 0 unit cell dimensions.We stress that this conclusion may well be particular to the mW model.Further calculations using models with atomic resolution may be required to draw unambiguous conclusions.

SUPPLEMENTARY MATERIAL
See supplementary material for full details of the symmetry adapted basis functions used to expand γ(θ, ϕ), the fitted coefficients, measured versus fitted interfacial stiffness, and free energy plots.

FIG. 2 .
FIG.2.Plot of the averaged local order parameter ξ for particles within volume elements running along columns perpendicular to an ice 1c (111)[1 10] interfacial ribbon geometry.Red crosses indicate ξ computed via the CHILL+ algorithm and black circles via the q 12 connections.Solid lines are tanh fits to the data, identifying the location of the interface as the point of inflexion at x = 205.0Å (black line) and x = 202.6Å (red line).

FIG. 3 .
FIG.3.Interface location in the ice Ic (111)[1 10] geometry.The computed interface shape is plotted as red points, each of which indicates the interfacial position at each of 50 discretised points in the z direction.Water molecules are coloured according to their number of solid-like connections, in this case identified via the q 12 order parameter.Particles with more than 12 connections are coloured grey, with the remainder coloured turquoise.Red beads mark the position of the interface (also identified by q 12 ) in each interfacial column.System dimensions are those reported in TableI.
) to the measured F(n) data reported in Ref. 31 at 215.2 K, we obtain a = −0.134,b = 0.120, and c = 7.404.

TABLE I .
Interfacial geometries used in calculating the stiffness of ice-water interfaces leading to reconstruction of the anisotropic interfacial free energy.Interfacial stiffness calculated via Equation (3) as a function of N z for three interfacial geometries.System dimensions are those reported in TableI.Red points indicate calculations using height profiles computed via the CHILL+ algorithm and black via the q 12 solidity criterion.

TABLE II .
Interfacial free energy at each of the three ice Ih surfaces studied.Results are reported for both the CHILL+ and q 12 methods of locating the solid-liquid interface.Results obtained by Espinosa, Vega, and Sanz 27 using the mold integration technique are shown for comparison.

TABLE III .
Calculated solid-liquid interfacial free energies of all ice Ic surfaces considered in this work.Results are shown for both methods of locating the solid-liquid interface when including m parameters (m − 1 symmetry adapted functions) in the fit to γc (θ, ϕ).

TABLE IV .
Calculated solid-liquid interfacial free energy for all planes of ice-0 considered in this work.