## Abstract

The melting of pure axisymmetric ice crystals has been described previously by us within the framework of so-called *geometric crystal growth*. Non-equilibrium ice crystal shapes evolving in the presence of hyperactive antifreeze proteins (hypAFPs) are experimentally observed to assume ellipsoidal geometries (‘lemon’ or ‘rice’ shapes). To analyse such shapes, we harness the underlying symmetry of hexagonal ice *I*_{h} and extend two-dimensional geometric models to three-dimensions to reproduce the experimental dissolution process. The geometrical model developed will be useful as a quantitative test of the mechanisms of interaction between hypAFPs and ice.

## 1. Introduction

Understanding how the basic tenets of surface thermodynamics control crystal shapes is a major focus of condensed matter science. Using a ubiquitous material such as ice as a test bed facilitates the advancement of these goals, while having immediate relevance to a wide range of applications. By independently controlling the temperature and pressure under which an ice crystal is grown, a range of equilibrium ice crystal shapes may be formed, from fully faceted surfaces to rough or rounded surfaces (Wettlaufer 2001; Cahoon *et al.* 2006).

It is rare to observe crystals in equilibrium; usually their shapes are probed in conditions of disequilibrium during which they are growing or melting, however slowly. Understanding the relationship between growth and equilibrium shapes (Elbaum & Wettlaufer 1993) is basic to unravelling the kinetic processes on surfaces and their relationship to lower-dimensional phase transitions, physisorption and chemisorption (Gilmer & Bennema 1972; Wettlaufer 2001). Equilibrium crystal shapes can be determined from the Wulff (1901) construction, which minimizes the surface energy for the fixed total volume of the sample. This geometrical construction requires as specified (or experimentally determined) the specific surface-free energy as a function of crystallographic orientation relative to the surface normal . For clarity, we note that for a liquid is the surface tension and hence is isotropic, independent of the surface normal and single valued. One then performs the transform where ** r** is a radial vector from the origin to the crystal surface. The equilibrium crystal shape emerges by taking the interior envelope of the set of planes lying perpendicular to radial rays that intersect the polar plot of . The overall size is characterized by

*λ*, but the shape is independent of it (Elbaum & Wettlaufer 1993). An equilibrium shape may be fully facetted (at low temperatures), it may be fully rough (at temperatures exceeding the roughening transitions for all of its facets) or both types of surface structure may coexist on a single crystal. Additionally, one might prepare an equilibrium crystal shape and use the Wulff construction to estimate the specific surface-free energy . However, owing to the long relaxation times of typical materials, it is experimentally challenging to observe equilibrium crystal shapes. Notable successes include crystals of helium in contact with its superfluid parent phase, due to the negligible latent heat of fusion and large thermal conductivity (Lipson 1987), micrometre-sized crystals of metals (Heyraud & Metois 1980), owing to their small sizes and the ability to maintain ultra-high vacuum conditions over long periods of time, and ice crystals in contact with water in a well-designed precision high-pressure apparatus (Cahoon

*et al.*2006; Maruyama 2011).

Knowledge of the kinetic processes of molecular attachment, detachment, and surface diffusion and incorporation are fundamental to predicting the evolution of a crystal interface on length scales much larger than the lattice parameter. Of significant basic interest and broad applicability is to gain an understanding of the relationship between the macroscopic shape of a pure ice crystal and the microscopic behaviour of water molecules. Moreover, the growth habits of ice can be influenced by impurities or additives such as antifreeze proteins (AFPs) that have specific affinity to certain crystallographic planes (Yeh & Feeney 1996; Wathen *et al.* 2003; Strom *et al.* 2004; Strom *et al.* 2005; Pertaya *et al.* 2007). Thus, quantitative studies of the non-equilibrium shapes of ice in the presence of AFPs may shed light on the kinetic processes underlying their action.

The driving force for both growth and dissolution is defined as the chemical potential difference between the parent phase and the solid, Δ*μ*. When the drive is sufficiently small that gradients in the rate-limiting diffusion field possess a characteristic length that is large relative to the size of the crystal, then interfacial kinetic processes control the shape of the moving interface. Such a circumstance is common in a wide range of systems and settings and is modelled within the framework of *geometric crystal growth* (Wettlaufer *et al.* 1994; Wettlaufer 2001; Tsemekhman & Wettlaufer 2003; Cahoon *et al.* 2006). One can consider geometric crystal growth as a kinetically grounded kinematic approach. This is because the local interfacial velocity is specified based on the local crystallographic orientation , and on the kinetics understood to be associated with it, for a given Δ*μ*. Such approaches are successful in two-dimensional growth having explained quantitatively the evolution of axisymmetric shapes (Wettlaufer 2001; Cahoon *et al.* 2006). Here, motivated by experiments involving AFPs, which substantially influence the interfacial kinetics and hence the overall crystal shape, we extend geometric approaches to three dimensions and focus specifically on the shapes of ice crystals during dissolution. Note that, while strictly speaking the shapes we model geometrically are dissolution shapes, here we use dissolution and the more intuitive term melting interchangeably. The distinction is irrelevant for either the actual shape or the method by which we model it.

The geometric approach explains the transient evolution of equilibrium forms containing both rough and faceted surfaces (Wettlaufer *et al.* 1994; Maruyama *et al.* 2000; Tsemekhman & Wettlaufer 2003). Because the rough orientations lack an interfacial nucleation barrier whereas the faceted orientations grow by an activated process, the general growth process is one whereby the rough orientations grow out of existence and leave the shape dominated by slow-growing facets. This is termed ‘global kinetic faceting’. In treating the melting of an initially faceted shape, growth–melt asymmetries seen in ice crystals were captured (Cahoon *et al.* 2006). There, an apparent rotation of the principal facets by 30^{°} was explained by the underlying kinetic anisotropy of due to step propagation. Pertaya *et al.* (2007) observed similar experimental phenomena in AFP solutions, attributing it to the inhibition of crystal growth by AFPs bound to ice crystal surfaces. While the prism facets in pure ice in contact with water disappear above −16^{°}C and at pressures below 160 MPa (Cahoon *et al.* 2006), hyperactive antifreeze proteins (hypAFPs) in solution created highly faceted hexagonal ice crystals viewed along the *c*-axis and ellipsoidal or ‘lemon-shaped’ prismatic ice crystal faces at normal pressures and at temperatures near 0^{°}C (figure 1). Modelling these complex three-dimensional shapes requires an extension of two-dimensional geometric models. Here, such an extension has been developed and allows us to simulate the observed hypAFP-modified shape evolution.

### (a) Geometric crystal growth

A small, spatially uniform and steady driving force characterizes near equilibrium growth or melting (Elbaum & Wettlaufer 1993), both of which are described by geometric crystal growth (Elbaum & Wettlaufer 1993; Wettlaufer *et al.* 1994; Maruyama *et al.* 2000; Wettlaufer 2001; Tsemekhman & Wettlaufer 2003; Cahoon *et al.* 2006). By ‘small’ it is meant that while the drive is large enough to change the shape of a crystal, if the drive is removed the shape relaxes to the equilibrium value on an experimentally observable time scale. As in the Wulff construction, where the equilibrium shape is determined by the orientation dependence of the specific surface-free energy , in geometric growth or melting the shape is determined by the orientation dependence of the local interfacial velocity , itself dependent upon the growth drive. For two-dimensional growth, one describes the surface itself by a time-dependent vector ** C**[

*x*(

*u*,

*t*),

*y*(

*u*,

*t*)], in Cartesian coordinates

*x*(

*u*,

*t*) and

*y*(

*u*,

*t*) depending on an arc-length parameter,

*u*, related to the arc-length

*s*by the metric in the usual manner;

*ds*=|∂

**/∂**

*C**u*|

*du*. Thus, for an inward pointing unit normal and an anisotropic velocity function , the surface evolves according to 1.1where we suppress the explicit representation of the Cartesian two-vector because the arc-length convention implies reference to a point on an evolving two-dimensional surface. Note that the time derivative is taken at fixed arc length. An equivalent means for dealing with interfacial kinetics as a function of is to characterize that direction in terms of the angle

*θ*that the surface tangent makes with the positive

*x*-axis, leading to

*V*(Δ

*μ*,

*θ*). For a given velocity function, the evolution of an initial shape is given by the solution to equation (1.1) using the method of characteristics (Wettlaufer

*et al.*1994).

Independent of whether one is describing the evolution of a two- or three-dimensional crystal, as long as the conditions of geometric growth are met, the global evolution of an initial crystal is governed by the physics embodied in a local velocity function *V* . However, a two-dimensional surface bounding a three-dimensional crystal must be parametrized by two arc-lengths modifying equation (1.1) as follows:
1.2where the arc-length parameters *u* and *v* correspond to the longitudinal and latitudinal positions on the evolving crystal surface. Thus, by parity of reasoning with the two-dimensional case, the arc-lengths *u* and *v* carry the same information as do the angles *θ* and *φ* that define the tangent plane perpendicular to the surface normal . The tangent plane is defined by the unit surface tangent , whose angle *θ* is measured from the positive *x*-axis, and the unit binormal , perpendicular to , which has a projection onto the equatorial plane taking an angle *φ* relative to the *x*-axis (figure 2). Whence, the local normal velocity has the following dependences *V* (Δ*μ*,*θ*,*φ*), it is intuitive that a unique correspondence between *u*, *v* and *θ*, *φ* exists globally for continuous convex shapes. When a discontinuity forms at a corner, the transformation is piecewise valid on the continuous regions separating corners. A surface that becomes faceted, and hence has zero curvature, must be treated separately such as described in Tsemekhman & Wettlaufer (2003), the results of which we extend presently, by rewriting equation (1.2) as
1.3described in more detail in the appendix. Owing to the fact that, for a given driving force, *V* (Δ*μ*,*θ*,*φ*) depends solely on the orientation of the normal, the solution ** C**(

*θ*,

*φ*,

*t*) of equation (1.3) is given along straight characteristic lines along which the normal is preserved 1.4All orientations on the initial surface

*C*_{0}(

*θ*,

*φ*) move along characteristic lines and at each time

*t*the locus of all points

**(**

*C**θ*,

*φ*,

*t*) associated with them defines the crystal shape. When two characteristics meet, a discontinuity in the surface slope can occur which can be described as a shock (Maruyama

*et al.*2000; Tsemekhman & Wettlaufer 2003). The trajectory of these shocks has been analysed in the two-dimensional case (Tsemekhman & Wettlaufer 2003). Simply reversing time allows for the treatment of melting.

Now we construct a simple velocity function *V* (Δ*μ*,*θ*,*φ*) geared towards our interests in the anisotropic melting shape of an ice crystal *viz.*,
1.5where *V*_{f} is the velocity normal to faceted surfaces, such as the basal plane, and *V*_{r} is that along rough surfaces, such as all prismatic planes which here become rough over time. The function *ξ*(*θ*) captures the transition from faceted to rough growth along the surface and is periodic with a period of *n*, and thus continuously varies between zero and one: 0≤*ξ*(*θ*)≤1. For simplicity, we chose , where *m* is an even integer. Motion of a molecularly rough interface is not hindered by activated processes, such as is the case at a facet, because for growth or melting there is no nucleation barrier; the surface can be thought of heuristically as a liquid–vapour interface. Thus, *V*_{r} is a relatively simple function of the growth drive Δ*μ*, which under experimental conditions is either proportional to the undercooling at fixed pressure or the overpressure at fixed temperature (Maruyama *et al.* 2000; Cahoon *et al.* 2006; Maruyama 2011). Using the intuition gleaned from both growth and melting studies of pure materials (Wettlaufer *et al.* 1994; Maruyama *et al.* 2000; Cahoon *et al.* 2006), we extend the linear response regime to the two-dimensional surface as follows:
1.6in which the linear response is shown in the first term with the growth drive scaled by a constant *c*_{r}. Moving along a surface from rough orientations towards vicinal orientations, the normal motion of an interface is influenced both by detachment of molecules along with the migration of surface molecules away from the facets (Wettlaufer *et al.* 1994; Cahoon *et al.* 2006). This behaviour is captured by the two angle-dependent terms, in which *p* and *q* are even integers, and *m*≥*p*,*q*.

In the experiments reported here and elsewhere (Maruyama *et al.* 2000; Cahoon *et al.* 2006; Maruyama 2011), it is seen that ice crystal shapes are twofold symmetric parallel to the basal plane and sixfold symmetric about the *c*-axis, and hence we take the simplest form of equation (1.6) that describes the symmetry of the crystal; *p*=*q*=2, *n*_{θ}=2 and *n*_{φ}=6 (figure 3). Under the conditions of geometric growth or melting, wherein the driving force is weak, *V*_{f}≪*V*_{r}, and hence, for example, the advance of the basal plane is much slower than that of the prismatic planes in experiments.

Given an initial shape, each point on the crystal advances according to the orientation of the surface. Depending on the initial shape, trajectories of two points may collide in a shock event (Maruyama *et al.* 2000; Tsemekhman & Wettlaufer 2003), according to the rubric of the method of characteristics. When a shock occurs, the participating characteristic orientations are eliminated from the crystal surface (Tsemekhman & Wettlaufer 2003). The right-hand side of equation (1.4) displays the wide range of possibilities, depending on the initial shape and the detailed form of *V* (Δ*μ*,*θ*,*φ*), for shock formation.

## 2. Results

Experiments on ice crystal growth in the presence of hypAFPs revealed highly faceted crystals in the basal plane (Bar *et al.* in press). AFPs inhibit ice crystals from growing within a certain temperature range, thereby decreasing the freezing point below the equilibrium melting temperature (Yeh & Feeney 1996). When crystals grown in solutions containing AFPs are continuously cooled and the magnitude of the freezing point is exceeded, sudden and rapid growth occurs due to the magnitude of supercooling imposed. In contrast, melting can be controlled using relatively small temperature gradients, and shape evolution is readily observed under a microscope (Bar *et al.* 2008). We used the three-dimensional model described earlier to simulate the melting process and to demonstrate how ice crystals melt in hypAFP solutions. The parameters were chosen to fit the simulation results to the experimental results of melting of crystals in *Tenebrio molitor* AFP (*Tm*AFP) solutions. We note that in such a manner comparison with other AFP solutions will offer future insight into system specificity versus generality of this approach.

Figure 4 shows model results of a dissolving ice crystal viewed parallel to the *c*-axis. As time progresses, the initially rounded shape evolves into a hexagonal structure and the faceted basal plane diminishes in overall area. This is a transformation commonly observed in nanolitre osmometer experiments with *Tm*AFP solutions (Bar *et al.* in press), distinct from the pure case in which the basal facet is maintained (Maruyama 2011).

Figure 5 displays the evolution perpendicular to the *c*-axis. The oblate spheroidal, or ‘drum-like’, initial shape forms ridges during melting. The relatively slow ablation normal to the basal plane creates a dynamically faceted oblate spheroid of the opposite aspect ratio that we refer to as a truncated lemon. Such an evolution is seen in the experimental observations shown in figure 1. Clearly, the morphology of these crystal surfaces is not arbitrary, but arises from a general set of underlying rules. Finally, as is seen particularly well in figure 5*c*,*d*, regions of the crystal become faceted during ablation and hence their local curvature must decrease. The evolution of the local curvature on certain regions of a crystal surface in which only one of the principal orientations is important can be analysed in ostensibly the same manner as in two dimensions (Wettlaufer *et al.* 1994; Maruyama *et al.* 2000; Tsemekhman & Wettlaufer 2003; Cahoon *et al.* 2006). For example, this might be treated by fixing the angle *φ* associated with the binormal and solving a local equation for the curvature *κ*_{θ} as a function of the angle *θ* associated with the tangent for a one-dimensional boundary, *viz*.
2.1In this equation, ∂/∂*τ* is taken at constant *θ*, the so-called ‘velocity stiffness’ is determined at fixed *φ*, which is, *mutatis mutandis*, the same as in Wettlaufer *et al.* (1994). Thus, for fixed *φ*, the evolution of the curvature as a function of *τ* from the initial value *κ*_{θI} is
2.2We note however that any local axisymmetry of the two-dimensional surface that allows this heuristic interpretation of the observed curvature decrease will in general be limited to a finite time interval. After such time it is possible that the curvature evolution in the other direction will intervene, as well as issues associated with the formation of shocks associated with regions on the surface with negative surface stiffness as discussed above and in the earlier studies (Wettlaufer *et al.* 1994; Maruyama *et al.* 2000; Tsemekhman & Wettlaufer 2003; Cahoon *et al.* 2006).

## 3. Conclusions

We have extended two-dimensional *geometric models* for crystal growth and melting to three dimensions and have examined the implications for the problem of ice crystals melting (dissolving) in solutions of hypAFPs. It is found that, consistent with observations of ice crystals in solutions of *Tm*AFP, initially rounded crystals melted to form stable faceted ‘lemon shapes’. The detailed morphology of such shapes is expected to depend on environmental conditions, such as the thermal hysteresis promoted by the AFP, the concentration of AFP, and the driving force Δ*μ*, all of which affect the underlying velocity function *V* (Δ*μ*,*θ*,*φ*). Indeed, having demonstrated the ability of geometric theories to capture detailed morphological changes in this complex system, the next stage of research aims at systematically testing the environmental influences on *V* (Δ*μ*,*θ*,*φ*). In this manner, refinement of both model structure and protein–ice interactions is possible. Importantly, such an approach provides quantitative insight concerning how macroscopic shapes are controlled by the mechanisms of dissolution and melting as they are modified by the presence of AFPs.

## Acknowledgements

This work has been supported by the National Science Foundation (NSF) under grant no. CHE-0848081, the Condensed Matter and Surface Science programme at Ohio University (CMSS), the Israel Science Foundation (1279/10 and 1281/10), the Marie Curie International Reintegration Grant (256364) and by the European Research Council (281595). J.J.L. thanks the National Natural Science Foundation of China and the Science (31160188) and Technology Foundation of Ministry of Education of China (211030) for their research support. J.S.W. thanks the Helmholtz Gemeinschaft Alliance ‘Planetary Evolution and Life’, the Wenner-Gren Foundation, the John Simon Guggenheim Foundation, and the Swedish Research Council for generous support of this research.

## Appendix A. Derivation of the three-dimensional model

The approach here is a straightforward extension of Tsemekhman & Wettlaufer (2003), and thus we provide here only the major waypoints sufficient to allow the reader to follow the logic. The evolution equation (1.2) is
A1where the arc-length parameters *u* and *v* correspond to the longitudinal and latitudinal positions corresponding to angles *θ* and *φ* that define the tangent plane perpendicular to the surface normal . The plane is collinear with the unit surface tangent , whose angle *θ* is measured from the positive *x*-axis, and the unit binormal , perpendicular to , which has a projection onto the equatorial plane taking an angle *φ* relative to the *x*-axis (figure 2). In terms of *u* and *v*, we write the usual arc lengths as *ds*_{u}=|∂** C**/∂

*u*|

_{v}

*du*and

*ds*

_{v}=|∂

**/∂**

*C**v*|

_{u}

*dv*, and with the use of A2a reparametrization of the surface is generalized from Tsemekhman & Wettlaufer (2003) as A3and A4Here

*κ*

_{θ}(

*θ*,

*t*) is a compact form of the curvature as a function of

*θ*and

*t*for fixed

*φ*and the same nomenclature follows, the obvious changes being made, for

*κ*

_{φ}(

*φ*,

*t*). Combining equations (A2)–(A4) yields A5The appropriate generalization of equation (1.4) from Tsemekhman & Wettlaufer (2003) gives the pair A6Combining equations (A5) and (A6) yields A7which is equation (1.3), the solution to which is equation (1.4).

- Received December 9, 2011.
- Accepted May 25, 2012.

- This journal is © 2012 The Royal Society