## Abstract

We present a linear stability analysis of a planar ice interface during unidirectional solidification of a hard-sphere colloidal suspension. We find that the interface can become unstable due to constitutional supercooling, yielding a new mechanism for pattern formation in colloidal systems. The interfacial stability is shown to depend strongly on the size and concentration of the particles. Increasing the particle radius tends to stabilize the interface, while increasing the concentration has a destabilizing effect. Additional effects that may influence the stability and morphology of such a system are described.

## 1. Introduction

Directional solidification of colloidal suspensions is of technological importance in a growing number of applications, including the fabrication of microaligned porous materials (Zhang *et al*. 2005), the cryopreservation of medical tissue and foodstuffs (Mazur 1970; Mashl *et al*. 1996) and the remediation of contaminated land (Gay & Azouni 2002; Northcott *et al*. 2005). An important property of colloidal systems is that they do not freeze uniformly. Instead the ice segregates, aligning the expelled particles into different configurations depending on the freezing conditions (figure 1). One of the major challenges for the future development of these technologies is to predict conditions under which the various ice patterns are expected to occur.

Previous approaches to modelling freezing particulate systems have generally focused on two aspects: the interaction of a particle with an advancing ice front (Uhlmann *et al*. 1964; Rempel & Worster 1999) and the freezing of rigid or elastic porous media (O'Neill & Miller 1985; Rempel *et al*. 2004). For example, Jackson & Chalmers (1958) considered a system in which ice is separated from a particle by a thin layer of water. This layer, now called a premelted film, is caused by long-range intermolecular interactions (e.g. van der Waals or electrostatic) between the ice and the particle surface (Dash *et al*. 2006). A temperature gradient causes some of the water to freeze, reducing the thickness of the film. Repulsive interactions between the ice and particle then push the particle towards warmer temperatures. Explicit calculations of the intermolecular forces between ice and various particle surfaces have led to quantitative predictions of particle and multiparticle dynamics (Wettlaufer & Worster 2006).

Contemporaneously with these studies have been developments of homogenized continuum theories focussing on the formation of discrete layers of ice such as those illustrated in figure 1*c* (O'Neill & Miller 1985; Fowler & Krantz 1994; Watanabe *et al*. 2000; Rempel *et al*. 2004). However, less has been done to account for the more complex ice patterns shown in figure 1.

Here, we take a different and complementary approach based on our modelling of freezing colloidal suspensions (Peppin *et al*. 2006). We previously obtained equations describing the solidification of a suspension of hard spheres and found that under certain conditions the suspension ahead of the ice–suspension interface can be below its freezing temperature (‘constitutionally supercooled’), leading potentially to an instability in the shape of the interface. The aim of the present work is to investigate the stability of the system and obtain conditions under which an initially planar ice–suspension interface will begin to deform. We find that increasing the particle radius tends to stabilize the interface, while increasing the concentration has a destabilizing effect. We obtain a regime diagram giving the critical wavelength at the onset of instability as a function of concentration.

## 2. Formulation

Figure 2 illustrates an experimental setup commonly used to study the solidification of colloidal suspensions (Mashl *et al*. 1996; Mutou *et al*. 1998; Watanabe 2002). A suspension of bulk volume fraction *ϕ*_{∞} is enclosed between two long parallel plates. The plates are laid across two heat exchangers, such that the freezing temperature of the suspension lies between *T*_{L} and *T*_{H}. The heat exchangers are kept fixed in place and the suspension is pulled towards the cold side at constant speed *V*. As the suspension freezes, a boundary layer of rejected particles builds up against the solidification front. In principle, a steady state can be reached in which the volume fraction *ϕ*_{s} of particles in the ice is equal to the volume fraction of the bulk suspension.

In this section, we first write down the conservation equations and boundary conditions describing the system and solve for the steady-state configuration. We then perturb the steady profile via normal modes, linearize the governing equations in the perturbed quantities and solve the resulting ordinary differential equation to determine conditions under which the interface is stable (perturbations decay in time) or unstable (perturbations grow in time).

### (a) Governing equations

We assume constant densities and thermal properties of ice, water and particles. Conservation of mass ahead of the interface can then be written, in a frame of reference moving with the pulling speed *V*,(2.1)where *h* is the height of the deforming interface above the steady planar position and *D* is the mutual diffusion coefficient, given by the expression (Russel *et al*. 1989)(2.2)where is the Stokes–Einstein diffusivity and(2.3)Here *k*_{B} is Boltzmann's constant, *T* is the absolute temperature, *R* is the particle radius, *μ* is the fluid dynamic viscosity and *Z* is the compressibility factor. For *Z*, we use the following hard-sphere equation of state (Peppin *et al*. 2006)(2.4)where *a*_{1}=4−1/*ϕ*_{p}, , *a*_{3}=18−10/*ϕ*_{p}, , and *ϕ*_{p}=0.64 is the volume fraction at random close packing. Equations (2.3) and (2.4) are based on sedimentation experiments using polystyrene beads and molecular dynamics simulations of hard spheres (Russel *et al*. 1989). Figure 3 shows *D* as a function of *ϕ* for several different particle radii.

Neglecting diffusion in the ice, conservation of mass at the interface is given by(2.5)where *ϕ*_{I} and *ϕ*_{S} are the volume fraction of particles on the unfrozen and frozen sides of the interface, respectively, is the speed of the deformed interface normal to itself, *x* is position along the interface and is the normal to the interface at *x*. The interfacial concentrations are related by the equation(2.6)where *k*_{s} is the segregation coefficient. In general, *k*_{s} depends on the interfacial concentration and on the solidification velocity *V*. For example, Rempel & Worster (1999) have shown that at small freezing velocities a spherical particle will be rejected by a growing ice phase (*k*_{s}≈0), while at larger velocities the particle will be engulfed (*k*_{s}≈1). Here, we consider relatively slow freezing so that the segregation coefficient is near zero and approximately constant. Far from the interface the boundary condition is(2.7)

One advantage of the experimental configuration in figure 2 is that the temperature field can be assumed constant in a frame of reference fixed with respect to the heat exchangers and given by the linear expression(2.8)where *T*_{i} is the temperature at the steady planar interface and *G*_{T} is the constant temperature gradient imposed on the system. Equation (2.8) is the ‘frozen-temperature approximation’ (*e.g*. Davis 2001) and is strictly valid when the thermal properties of ice and suspension are the same, and the growth velocity is slow enough that the effects of latent heat can be neglected at the freezing interface. Though these conditions are not well met in ice/water systems, we employ the frozen-temperature approximation here for mathematical convenience to expose the principle mechanisms for instability straightforwardly. This limitation will change only the magnitude of the effective temperature gradient in the suspension and not the nature of the instability in the linear regime.

When local equilibrium prevails, the temperature at the planar interface is equal to the thermodynamic freezing temperature *T*_{f}(*ϕ*), which for a hard-sphere suspension is given by (Peppin *et al*. 2006)(2.9)Here, *ϕ*_{i}=*ϕ*_{∞}/*k*_{s} is the volume fraction at the steady planar interface, *T*_{m} is the freezing temperature of the pure fluid and , where is the volume of a particle, *ρ* is the density (assumed constant) and *L*_{f} is the latent heat of fusion of ice. The temperature at the deforming interface is(2.10)where *ϕ*_{I} is the volume fraction at the deforming interface, is twice the mean curvature of the interface and *γ* is the interfacial surface energy. The surface energy term on the right-hand side of equation (2.10) represents the Gibbs–Thomson effect, i.e. the change in temperature at the interface due to the curvature there. For equation (2.10) to apply macroscopically, the radius of curvature of the interface must be much larger than the radius of a particle: |1/*K*|≫*R*. Here, we assume that the surface energy does not depend on the particle concentration or the local orientation of the interface, and treat *γ* as an averaged constant. We return to this issue later in the paper.

By neglecting the slight temperature dependence of *D*_{0}, we can introduce the dimensionless variables , , , and where *δ*_{L}=*D*_{0}/*V*, *δ*_{t}=*D*_{0}/*V*^{2}, *G*_{T} is the temperature gradient and *G*_{ϕ} is the concentration gradient at the planar interface obtained from equation (2.5) as(2.11)Equations (2.1)–(2.8) become, upon dropping the carets,(2.12)(2.13)(2.14)(2.15)(2.16)where *M* is the morphological number,(2.17)*Γ* is a surface energy parameter,(2.18)and is the gradient in the freezing temperature at the planar interface,(2.19)

### (b) Basic state

In the steady planar case, *h*=0 and equations (2.12) and (2.13) become(2.20)(2.21)where *C*_{i}=*C*_{∞}/*k*_{s}. Equation (2.20) can be integrated from *z*=0 to *z*=∞ using equation (2.21) as a boundary condition to yield(2.22)Equation (2.22) can be solved using the boundary condition *C*=*C*_{i} at *z*=0 to yield an implicit expression for *C*(*z*). In the basic state, *T*(*z*) is given by equation (2.16). Figure 4 shows steady-state temperature and concentration profiles for three different interfacial concentrations. In figure 4*a*,*d*, the concentration gradient, and hence the freezing temperature gradient, at the interface is small and the temperature ahead of the interface is always warmer than the freezing temperature. In figure 4*c*,*f*, the temperature ahead of the interface is below the freezing temperature and the suspension is constitutionally supercooled and potentially unstable.

### (c) Linear stability

The variables describing the steady state are now perturbed (except for *T* which is assumed steady because of the disparity between thermal and particle diffusivities)(2.23)where and are the solutions to the steady planar case obtained in §2*b*. The perturbations *C*′ and *h*′ have the form of normal modes (Chandrasekhar 1961; Davis 2001)(2.24)where *σ* is the growth rate of the disturbance and *α* is the wavenumber of the normal modes along the interface. The growth rate *σ* determines whether or not a particular disturbance will grow in time. Assuming exchange of stabilities, the condition Re(*σ*)=0 defines a locus of parameter values separating stable (Re(*σ*)<0) from unstable (Re(*σ*)>0) states. In the absence of exchange of stabilities, there exists the possibility that the imaginary part of *σ* is non-zero, in which case the interface will oscillate in time.

Inserting equations (2.23) and (2.24) into equations (2.12)–(2.16) and linearizing in the perturbed quantities leads to the following nonlinear ordinary differential equation for the disturbance amplitude:(2.25)where *D*_{C} is the derivative of the diffusion coefficient with respect to concentration. The boundary conditions on (2.25) are(2.26)and(2.27)where(2.28)

Equations (2.25)–(2.28) represent an eigenvalue problem that can be solved numerically using the shooting method (Press *et al*. 1992).

## 3. Results and discussion

Before solving the full equations (2.25)–(2.27), it is instructive to consider the dilute limit so that *D*→1. In this case, equations (2.25)–(2.27) can be solved analytically to yield the characteristic equation(3.1)Equation (3.1) is similar in form to the equation derived in the case of dilute alloys (Mullins & Sekerka 1964; Davis 2001), with the important distinction that here the coefficients depend on the particle radius. We first concentrate on the case in which *σ* is zero. Equation (3.1) can then be solved for the morphological number *M* to yield(3.2)Figure 5 shows plots of *M* versus the wavenumber *α* for *k*_{s}=1×10^{−6} and various particle radii. Below the curve the interface is stable, while above the curve it is unstable. The minimum in the curve gives the critical morphological number at which instability first sets in. The physical origin of this instability is analogous to the alloy case (Tiller *et al*. 1953; Mullins & Sekerka 1964). A perturbation of the interface causes a portion of the ice phase to be further into the suspension. In the absence of supercooling, such a perturbation would melt back. If the suspension is supercooled, however, the perturbation will tend to grow into the supercooled region and the interface is unstable.

Figure 6*a* shows the solution of equations (2.25)–(2.27) for increasing particle concentrations. As the concentration increases, the interface becomes more unstable. Also shown is the dilute result (dashed line) and for *ϕ*_{i}<0.01 there is good agreement. At higher concentrations, the full numerical solution must be used. Figure 6*b* is a regime diagram giving the critical wavelength at the onset of instability as a function of concentration for several particle radii. This diagram estimates dendrite spacings as a function of parameters, such as concentration and particle size, yielding valuable insight into the ice morphologies expected in freezing hard-sphere suspensions.

As noted by Mullins & Sekerka (1964) and others (Kurz & Fisher 1998; Davis 2001), the surface energy term in equation (3.1) tends to stabilize the interface. In the present case of a hard-sphere suspension *Γ* is proportional to *R*^{5}, so that increasing the particle radius tends to quickly stabilize the interface. Increasing *R*, however, reduces *D*_{0} and hence the particles are less able to diffuse relative to the moving interface. It is likely that increases in *R* will eventually lead to loss of local equilibrium at the interface. That is, as the particle radius increases at a given freezing velocity, non-equilibrium effects such as viscous flow at the interface become more important, leading eventually to particle engulfment (Uhlmann *et al*. 1964; Chernov *et al*. 1976; Rempel & Worster 1999). In the future, it may be possible to model this behaviour by allowing the segregation coefficient *k*_{s} to depend on *R* and *V*. It has been shown by Coriell & Sekerka (1983) that if *k*_{s}=*k*_{s}(*V*), the interface can oscillate in time (Im(*σ*)≠0), thereby trapping impurities (particles) during the rapid part of the cycle and rejecting them during the slow part. This type of oscillatory motion of the interface has been observed during the solidification of corn starch suspensions (Mashl *et al*. 1996) and silica microspheres (Watanabe 2002).

Before concluding, we briefly discuss some limitations of the theory and potential directions for future research. First of all, we note that salt ions were likely present in the systems illustrated in figure 1. Such impurities can themselves cause a morphological instability of the ice interface (Mullins & Sekerka 1964). The effect of this instability on concentrated colloidal suspensions remains to be determined. We note that, due to the presence of a Debye layer, a likely influence of solutes is to induce electromigration of the particles and hence change the effective diffusion coefficient. Our purpose here was to show that the colloidal suspension can become unstable *even in the absence* of electrolytes. This basic result can then serve as a benchmark for more complex systems. Experiments on monodisperse silica or polystyrene microspheres dispersed in an ultrapure liquid would provide the most rigorous test of the theory.

Another factor complicating the application of the theory to ice systems is the highly anisotropic crystal growth of ice, with much slower growth perpendicular than parallel to the basal plane (Wettlaufer 2001). This effect could, in principle, be modelled via a surface energy that depends on orientation (Koo *et al*. 1991; Davis 2001). We believe that anisotropy is a crucial aspect of a quantitative explanation of the features in figure 1. Other directions for future research include accounting for the different thermal properties of ice and suspension, the temperature dependence of the diffusion coefficient, diffusion and thermal diffusion (regelation) in the ice phase, a velocity-dependent segregation coefficient, nonlinear effects in the stability analysis and hydrodynamic effects in the suspension.

## 4. Conclusion

The freezing interface in a hard-sphere suspension can become unstable due to constitutional supercooling. This is an important new mechanism for ice segregation and pattern formation in freezing colloidal suspensions. Increasing the particle radius appears to have a stabilizing effect, while increasing the concentration is destabilizing. We have obtained a regime diagram indicating the critical wavelength at the onset of instability as a function of concentration and particle radius. Finally, we have laid bare a number of outstanding issues and opportunities for future work, including *inter alia*, nonlinear effects, the influence of electrolytes and crystallographic anisotropy.

## Acknowledgments

This research was supported by a Gates Cambridge Scholarship, a grant from Trinity College and funding from the Leverhulme Trust. J.S.W. acknowledges support from the US National Science Foundation (0PP0440841), the Department of Energy (DE-FG02-05ER15741) and both Trinity College and DAMTP at the University of Cambridge where this work evolved.

## Footnotes

- Received August 18, 2006.
- Accepted October 26, 2006.

- © 2006 The Royal Society