A new phase-field model for strongly anisotropic systems

We present a new phase-field model for strongly anisotropic crystal and epitaxial growth using regularized, anisotropic Cahn–Hilliard-type equations. Such problems arise during the growth and coarsening of thin films. When the anisotropic surface energy is sufficiently strong, sharp corners form and unregularized anisotropic Cahn–Hilliard equations become ill-posed. Our models contain a high-order Willmore regularization, where the square of the mean curvature is added to the energy, to remove the ill-posedness. The regularized equations are sixth order in space. A key feature of our approach is the development of a new formulation in which the interface thickness is independent of crystallographic orientation. Using the method of matched asymptotic expansions, we show the convergence of our phase-field model to the general sharp-interface model. We present two- and three-dimensional numerical results using an adaptive, nonlinear multigrid finite-difference method. We find excellent agreement between the dynamics of the new phase-field model and the sharp-interface model. The computed equilibrium shapes using the new model also match a recently developed analytical sharp-interface theory that describes the rounding of the sharp corners by the Willmore regularization.


Introduction
The formation of faceted pyramids on nanoscale crystal surfaces is an important phenomenon that has been attracting wide attention owing to its role in the self-organization of quantum dots. In crystalline films, these structures arise through competition among surface and bulk forces that result in instability. Instability may originate from a lattice misfit between film and substrate, from strong surface anisotropies or from kinetic surface fluxes. Facets on crystalline thin films, however, arise for a thermodynamic reason. In particular, such facets arise because the surface free energy is non-convex with respect to the surface normal (Herring 1951). Quantitative modelling of self-organization processes thus requires a detailed description of surface energies and mass transport mechanisms along the crystal surface.
Thermal faceting (spinodal decomposition) of thermodynamically unstable crystal surfaces caused by strongly anisotropic (non-convex) surface free energy densities and driven by surface diffusion has been considered in Stewart & Goldenfeld (1992), Liu & Metiu (1993) and Savina et al. (2003). However, all of these theoretical and numerical treatments are restricted to long-wave approximations based on small variations in surface orientation. This introduces a clear limitation in their quantitative predictive power as many experimentally observed facet angles in thin crystalline films are not small.
Only recently have numerical methods been proposed to deal with the full geometric evolution. The evolution law follows from a non-dimensional surface free energy of the form (Gurtin & Jabbour 2002) E½G Z ð G gðnÞ C b 2 H 2 dG; ð1:1Þ with g(n) the classical surface free energy density with n the surface normal, H the mean curvature and ffiffiffi b p a small scale over which corners are smeared out. The energy can be viewed as a geometric Ginzburg-Landau-type energy with a nonconvex functional and a gradient term. This is most apparent in the one-dimensional setting, where the energy can be written as E½GZ Ð G gðqÞC ðb=2Þ j v s q j 2 dG with q the tangent angle (angle between the tangent vector s and the x -axis) and s is the arc length (Herring 1951;DiCarlo et al. 1992). Thus, it is also possible to interpret the gradient (curvature) term in equation (1.1) as the next higher order term in a more general surface free energy densitygZgðn; H ; .ÞZ gðnÞC ðb=2ÞH 2 C/ : The model we investigate here is the H K1 gradient flow of the energy (1.1). This leads to the (non-dimensional) surface diffusion equation where V is the surface normal velocity; D G is the surface Laplacian; n is the mobility; H g is the weighted mean curvature; and K is the Gaussian curvature. The weighted mean curvature is defined through H G ZV$x, where xZDg(n) is the Cahn-Hoffman (1974) vector. When the surface energy is sufficiently anisotropic, equation (1.2) is ill-posed for bZ0. When bO0, Siegel et al. (2004) and Hausser & Voigt (2005a) solved this highly nonlinear sixth-order equation numerically. However, both approaches are restricted to curves as severe numerical problems occur for evolving parametric surface meshes. A level-set approach circumventing some of these difficulties is discussed in Burger et al. (2007). Here, we are interested in a phase-field approximation of the problem. The phase-field method is capable of describing the evolution of complex, topology-changing surfaces and also provides a general framework to add further physical effects such as elasticity or compositional differences. In the phase-field approach, sharp interfaces are given a finite thickness and the energy (1.1) is modified appropriately. As such, this is also known as a diffuse-interface method.
Attempts towards a phase-field approximation have been made in Eggleston & Voorhees (2002), Wise et al. (2005Wise et al. ( , 2007,  and Wheeler (2006). However, none of these models approximate the geometric evolution described above in equation (1.2) correctly as the diffuse-interface width tends to zero. In an alternative phase-field approach, Eggleston et al. (2001) and Eggleston & Voorhees (2002) dealt with the ill-posedness of the anisotropic problem by convexifying the energy. This method does not reproduce equation (1.2) in the sharp-interface limit, and it also modifies the physics by suppressing the nucleation of new facets if the film is started with an orientation within the nonconvex region (see Gurtin & Jabbour 2002;Fried & Gurtin 2004). Wise et al. (2005Wise et al. ( , 2007 and Wheeler (2006) used a regularization of the ill-posed equation based on adding the square of the Laplacian of the phase-field variable to the energy (instead of the mean curvature squared as shown above). This leads to a rounding of corners and edges and allows for spinodal decomposition of unstable orientations into stable ones. However, the asymptotic limit of the associated phase-field/diffuse-interface system is not equation (1.2) and the dynamics and equilibrium shapes of the two models are different.
A correct phase-field approximation of equation (1.2) requires an approximation of the Willmore energy. Attempts to construct phase-field approximations for the Willmore energy date back to De Giorgi (1991), who conjectured an appropriate approximation, which has since been analysed and simplified using different methods by Loreti & March (2000), Du et al. (2004Du et al. ( , 2005, Röger & Schätzle (2006) and Wang (2007). An alternative approach has been used by Biben & Misbah (2003), Biben et al. (2005) and Jamet & Misbah (2008). In , the De Giorgi approach to approximating for the Willmore energy, which we follow in this paper, is proposed as a regularization for the strong anisotropic surface energy in a phase-field model. Adaptive numerical simulations using this approach were performed by Wise et al. (2007). However, the specific forms of the Willmore term and the phase-field/diffuse-interface approximation of the anisotropic surface energy term used in these works were incompatible. Thus, this approach also does not give equation (1.2) as the asymptotic limit. The incompatibility arises because of an inappropriate combination of anisotropic and isotropic terms in the energy. In particular, the orientation-dependent interface thickness introduced by the surface energy anisotropy (Kobayashi 1993) is not compatible with the Willmore approximation that assumes a uniform interface thickness. This motivated us to develop a new way of approximating anisotropic free energies in phase-field models. A key feature of our approach is a diffuse-interface thickness that is independent of the orientation. This makes it straightforward to combine isotropic and anisotropic terms. The resulting system is highly nonlinear and is sixth order in space.
Using matched asymptotic expansions, it can be shown that our new model converges to the classical sharp-interface model in the limit of vanishing interface thickness (see the electronic supplementary material). We present simulations in two and three dimensions, which show excellent agreement between the dynamics and equilibrium shapes simulated by the sharp-interface model and by our new phase-field model. At equilibrium, the computed shapes also match the results of an analytical sharp-interface theory by Spencer (2004), which describes the rounding of sharp corners by the Willmore regularization.
The paper is organized as follows. In §2, we discuss equilibrium shapes of crystals with strong anisotropies. In §3, we derive the phase-field model for the Willmore regularization of strongly anisotropic surface diffusion. In §4, we briefly describe the numerical scheme used to solve the highly nonlinear sixth-order phase-field/diffuse-interface equations. In §5, we present simulations of the evolution of facetted curves and surfaces. We compare the equilibrium shapes with the asymptotic results in Spencer (2004). We also demonstrate mound formation and the coarsening of surface structures. Finally we draw conclusions in §6. The formal matched asymptotic expansion showing the convergence to the sharp-interface model is given in the electronic supplementary material.

Equilibrium shapes
The equilibrium (Wulff ) shape of a crystal is defined as the shape that minimizes the surface free energy E[G] under the constraint of fixed enclosed volume. Depending on the details of the free energy density g, the equilibrium shape may contain corners and edges. In such cases, it is energetically favourable to exclude high-energy orientations and allow for missing orientations. In two dimensions, the surface energy g(n) can be written as gZg(q), where q is the tangent angle. Corners, edges and missing orientations occur if the stiffness becomes negative for some orientations. The stiffness is defined asgZ gC g qq . In three dimensions, with gZg(q, f) depending on two angles, no such analytical criterion is known. Instead, the polar plot of 1/g has to be used and it has been shown in Sekerka (2005) that the onset of missing orientations occurs at a convex-to-concave transition in the 1/g plot. Using the Cahn-Hoffman (1974) x-vector, missing orientations occur when 'ears' and 'flaps' form in the x-plot (Sekerka 2005). Consider the fourfold symmetric model-type anisotropy where a is the anisotropy strength and dZ2, 3 is the space dimension. In two dimensions, this is equivalent to gðnÞ Z 1 C a cos 4q: ð2:2Þ In this paper, we use these anisotropy functions exclusively although our approach of course works for any such function. The anisotropy function (2.1) is convex for values a!1/15 and non-convex for aR1/15. In figure 1, the corresponding Wulff shapes that are plotted (in three dimensions the x-plot is shown), which in three dimensions resemble a double-sided pyramid, are shown for convex (aZ0.06) and non-convex (aZ0.2 and 0.3) surface energy functions g from equation (2.1). Note that the actual Wulff shapes do not contain the unphysical ears and flaps (that appear in the three-dimensional x-plot) and instead contain corners and missing orientations (e.g. Sekerka 2005).
The H K1 gradient flow of the corresponding surface energy E½GZ Ð G gðnÞ dG defines the following model for anisotropic surface diffusion: where D G is the surface Laplacian (Laplace-Beltrami operator) and n is the mobility. In two dimensions, D G Zv ss , where s is the arc length, and vE vG Z ðg C g 00 ÞH ; ð2:4Þ where g 0 Zdg/dq. When the anisotropy is sufficiently strong such that there are missing orientations in the Wulff shape, the evolution by surface diffusion is inherently unstable. In fact, the evolution equations are actually ill-posed because the equations are backward parabolic for these orientations owing to the non-convexity of the anisotropic surface energy. One way to overcome this ill-posedness is to regularize the equation by adding a curvature-dependent term to the interface energy. This was already proposed on physical grounds in Herring (1951), and later mathematically introduced in DiCarlo et al. (1992) and Gurtin & Jabbour (2002). Such a curvature-dependent term introduces a new length scale on which sharp corners and edges are rounded. In , it is argued that even higher order terms in the energy are necessary to prevent a surface from forming corners and edges; however, here we will only consider energies as in equation (1.1). Minimizing the surface energy thus becomes a compromise between a large curvature at the corners and edges, which decreases orientations with large surface energy but increases the regularization term, and small curvature at the corner, which decreases the regularization term but increases orientations with large surface energy. The amount of corner rounding is therefore determined by these two competing energy terms. The plausibility of such a regularization is clear, but its effect on the equilibrium shape was only recently analysed for curves and is still open for surfaces. For curves, Spencer (2004) gave an asymptotic analysis that shows the convergence to the sharp-corner (Wulff ) results as the regularization parameter b/0 and provides an analytical formula for the equilibrium shape near a rounded corner for bO0 (figure 2). We will use this asymptotic solution, referred to as a regularized Wulff shape, later to compare with equilibrium corner profiles for the new phase-field/diffuse-interface model.

Phase-field/diffuse-interface approximation
In the phase-field/diffuse-interface approach, narrow transition layers replace sharp interfaces and an order parameter 4 that denotes the phases in a multiphase system is introduced. The order parameter is constant (0 or 1) in each phase and the interface (transition layer) between phases corresponds to the region where 4 varies from 0 to 1. The classical formulation of the phasefield/diffuse-interface equations for isotropic systems is based on the energy (Cahn & Hilliard 1958) where e is a small parameter that is a measure of the interface transition layer thickness and fZf(4) is a double-well potential. A phase-field approximation for surface diffusion based on this energy was introduced by Cahn et al. (1996). The evolution law is where M(4) is the mobility which is localized near the interface, e.g. M(4)Z 44(1K4). Equation (3.3) is a degenerate Cahn-Hilliard equation that converges as e/0 to motion by surface diffusion VZD G H. An improved version has been given in , where m is replaced by Gm, with GZG(4) localized near the interface. A detailed asymptotic analysis for this and other approaches, not discussed here, is given in Gugenberger et al. (2008). In the isotropic case, numerical simulations indicate that using Gm rather than m does not significantly influence the evolution (A. Rätz & A. Voigt 2008, personal communication).
Assuming that the double-well potential is f(4)Z4 2 (1K4) 2 /4, then in all cases the profile across the interfacial layer is where d denotes the signed distance to the centre of the interfacial layer. Following the classical approach to incorporate anisotropy (Kobayashi 1993), the anisotropic energy is where the anisotropic gradient m is given by m Z g 2 ðnÞ V 4 C gðnÞ jV 4 j PV n gðnÞ; ð3:8Þ and the unit normal n and the projection matrix P are given by and P Z I Kn5n; where I is the identity matrix. In equation (3.8), V n represents the gradient with respect to the components of the normal vector. In , the convergence of this model to motion by anisotropic surface diffusion VZD G H g is shown by the method of matched asymptotic expansions as e/0. In this approach (see also Wheeler et al. 1996), the interface thickness is seen to be a function of orientation Starting from these equations various attempts have been made to extend the approach outlined above to strongly anisotropic (non-convex) surface free energy densities g. As in the sharp-interface case, equations (3.6)-(3.8) become ill-posed (e.g. Wise et al. 2007). Wise et al. (2005Wise et al. ( , 2007 and Wheeler (2006)  to the energy. In , an approach is proposed which uses instead a De Giorgi-type phase-field approximation of the Willmore energy as a regularization. In this case, the additional term in the energy is b 2 ð U 1 e 3 ðf 0 ð4ÞKe 2 D4Þ 2 dU: ð3:11Þ This was implemented and investigated in Wise et al. (2007). Such phase-field/ diffuse-interface approximations of the Willmore energy have been shown to converge to the Willmore functional in the limit e/0 (see Loreti & March 2000;Du et al. 2004Du et al. , 2005Röger & Schätzle 2006;Wang 2007). However, it can be seen that the Willmore term leads to a constant interface thickness and is thus not compatible with the phase-field/diffuse-interface approximation of anisotropy as outlined above. In an alternative approach to approximating the Willmore energy, used in Biben & Misbah (2003), Biben et al. (2005) and Jamet & Misbah (2008), the curvature may be directly calculated by HZV$n, with n given as above, and the integrand in equation (3.11) may be replaced with H 2 jVfj. While this approach is feasible to use here, it leads to additional nonlinearity compared with the De Giorgi approach. We now reformulate the anisotropic energy (3.5) such that the interface thickness is independent of orientation. For this reason, we consider the anisotropic surface energy as Note that this form of the energy differs from the form in equation (3.5) in that g multiplies both the f and jV4j 2 terms. This is a natural formulation if we interpret the term ð1=eÞð f ð4ÞC ðe 2 =2Þ jV 4 j 2 Þ as approximating the surface delta function as in the original Cahn-Hilliard system. The corresponding evolution law becomes To further simplify, we may use the asymptotic result that near interfaces f ð4Þ w ðe 2 =2Þ jV 4 j 2 thereby giving m wgðnÞV 4 C jV 4 j PV n gðnÞ; ð3:15Þ which is different from equation (3.8) by a factor of g. As is shown in the electronic supplementary material, the resulting interfacial layer is independent of orientation and is given by the isotropic formula (3.4). This allows the anisotropic surface energy (3.12) to be combined consistently with the phasefield approximation of the Willmore regularization. Our regularized free energy thus reads This is now a sixth-order system. We should point out that this evolution does not strictly guarantee that dE/dt%0 because of the asymptotic approximation (3.15) of m. If we did not use this approximation, the energy would be nonincreasing. However, the term f (4)/jV4j can be problematic when e is not sufficiently small, since in principle jV4j could vanish in regions where f does not. So using the approximation, we have to check the behaviour of the energy numerically. In all cases, we find that the energy is non-increasing. In the electronic supplementary material, we use the method of matched asymptotic expansions (following Pego 1989;McFadden et al. 1993;Leo et al. 1998;Loreti & March 2000), which shows that the system (3.17)-(3.20) converges to anisotropic surface diffusion as e/0.

Numerical methods
To numerically solve the sixth-order system derived in §3, a straightforward modification of the algorithm developed by Wise et al. (2007) suffices. We use the mobility and free energy functions given by where t is a small positive parameter (tZ10 K6 ) introduced to keep the mobility from vanishing. We solve a coupled system of three second-order equations for 4, m and u using second-order accurate-centred finite differences in space. A Crank-Nicholson scheme is used to discretize in time. The anisotropic gradient is discretized in conservation form in space and is lagged in the multigrid V-cycle.
To solve the resulting nonlinear equations, a nonlinear multigrid method is used. We further use adaptive block-structured Cartesian mesh refinement (see Wise et al. (2007) for details).

Numerical results
We begin by comparing the results we have obtained from our new phasefield/diffuse-interface model (3.17)-(3.20) with the sharp-interface model (2.3) at different stages of the evolution in two dimensions. In order to study the sharpinterface model described in §2 dynamically, we use a non-stiff front tracking technique (e.g. Hou et al. 1994;Leo et al. 2000;Siegel et al. 2004;Li et al. 2005). The idea is to use a high-order (spectrally accurate) parametrization of the interface together with a special tangential velocity to enforce dynamic equal spacing in arc length and a non-stiff time integration algorithm for the interface tangent angle. We consider the evolution of a circle towards equilibrium under surface diffusion with a strong fourfold anisotropy given in equation (2.1), with aZ0.1, Willmore regularization bZ0.01 and interface thickness parameter eZ0.03. The phase-field description of the initial circle is 4ðx; y; t Z 0Þ Z 1 2 1:0Ktanh The phase-field computational domain is UZ(0, 6.4)!(0, 6.4). The simulation is performed with a root-level grid size hZ0.1 (64!64 grid points); there are two levels of refinement and the finest resolution is hZ0.025. There are six or seven points across the interface during the simulation. The time step is sZ0.001.
In the sharp-interface model, the parameters a, b and n are matched to the new phase-field model using the asymptotic formulae in the electronic supplementary material. In the sharp-interface simulation, there are 256 nodes on the interface and the time step size is sZ0.125!10 K8 .
In figure 3, we present snapshots of the evolving 4Z0.5 level set (solid curve) at different times together with the sharp-interface results (dashed curve). The corresponding adaptive meshes and 4Z0.5 level sets are shown in the figure 3(ii). The initial circular interface transforms under the evolution to a fourfold symmetric shape with slightly rounded corners at tZ100, which is almost the equilibrium shape. There is an excellent agreement between the new phase-field model and the sharp-interface model throughout the evolution. In figure 4, we compare the corresponding energies for the two models, where the scaling of the sharp-interface energy is consistent with the asymptotic formulae in the electronic supplementary material. Both energies are non-increasing in time and there is excellent agreement between the new phase-field (solid curve) and the sharp-interface (dashed curve) results.
Next, we compare the new phase-field model results near equilibrium with the sharp-interface Wulff shape and the regularized Wulff shape near the corner (Spencer 2004). To emphasize the effects of anisotropy, we consider here aZ0.3. The initial shape is a circle and the domain and all other parameters are the same as in figure 3. We also compare the results of the new phase-field model (figure 5a) with previous phase-field models resulting from the classical anisotropy model equation (3.5) combined with either the Willmore regularization (3.11) (figure 5b) or the linear regularization (3.10) (figure 5c). Observe that the linear regularization not only affects the corner angle but also results in sides that are more curved than the Wulff shape and the Willmore regularization. The 4Z0.5 level curve is plotted at time tZ200 where the corresponding energy profile in figure 6 suggests that the shapes are nearly equilibrated. The Wulff and regularized shapes are shown as the black and grey dashed curves, respectively. Near the corners, in contrast to the Wulff shape, both the phase-field models with Willmore regularization yield interfaces with smoothed corners, which match the regularized Wulff shape. This is the first time we are aware of that such a comparison has been performed in the phase-field/diffuse-interface context. For the new model, away from the corners, there is very good agreement between the numerical results and Wulff shape. But there is a mismatch between the numerical results from previous phase-field models and Wulff shape away from the corners. It seems that the angle between the two sides is larger in these models than that in the Wulff shape. This deviation is real and is not due to the possibility that the shape is still evolving since the energy (shown in figure 6) indicates that the shape is very close to equilibrium. To quantify the deviation, in figure 7, we examine the orientation angle q (tangent angle) as a function of arc length s near the corner. Observe that the new and previous phase-field models with Willmore regularization match the regularized Wulff shape (Spencer 2004) at the corner, where qZ0, but away from the corner only the new model matches the asymptotic solution and hence the Wulff angle. The previous phase-field model with the linear regularization gives a different result both near and away from the corner. Note that the (Wulff ) angle between the two corner sides in the shapes from figure 5 is equal to pK[max(q)Kmin(q)]. The deviation between the previous models and the Wulff angle is apparent.  In figure 6, the energy is shown for the new (solid curve) and previous phasefield models with Willmore (dashed curve) and linear (dot-dashed curve) regularizations. The total energy for all the models is decreasing. At tZ40, the decrease in energy becomes very small and after some time the energy is almost constant. Thus, we consider tZ200 to be very close to equilibrium. Also note that the energies for the previous models are larger than that of the new model, with the linear regularization being quite a bit larger. It is interesting to look at the energy at early times (figure 6b) where the energies for the models with Willmore regularization are shown. As indicated by arrows, there is a short period of time where the rate of decrease in energy is fairly small. Right before this time, there is a rapid drop in the energy. This phenomenon occurs in all the phase-field models so we focus only on the new model here. As we show, this rapid decrease is associated with the morphological evolution of the shape. In figure 8, snapshots from the evolution of the initially circular shape are shown at early times. By time tZ1.5, a three-hill/two-valley structure is formed. This is seen clearly in the close-up and orientation plots shown in figure 8b,c. Note that the hills and valleys correspond to the points where qZ0; in our orientation v s q!0 at qZ0 corresponds to a hill. This hill-valley structure emerges because the strong thermodynamic driving force to remove the highenergy orientations on the circular interface creates instability, manifest by the hill-valley structure, which is mediated by the (Willmore) regularization. The drop in the energy seen in figure 6b before the arrow corresponds to the formation of the hill-valley structure. Next, the hill-valley structure further evolves to form the single hill structure seen at tZ3.5; the energy decreases slowly during this process. To characterize this in more detail, we plot in figure 9 the x coordinates of the hills (circles) and valleys (pluses) at different times. A single hill forms first by tZ0.4. Then, by tZ0.8, the three-hill/two-valley structure emerges. Then, shortly after tZ1.5, the hill and valley located around xz2.8 merge leaving a two-hill/one-valley structure remaining. The hill and valley located near xz3.6 then merge just after tZ1.9 leaving the single hill located at xz3.2. We now investigate the effect of the Willmore regularization parameter b on the equilibrium interface morphologies using the new phase-field model. We consider three different values, bZ0.004, 0.01 and 0.018. The other parameters are aZ0.1 and eZ0.03. The mesh and time step are the same as in the previous simulations. Using the energy as a guide (not shown), we use the time tZ100 shapes as approximations of equilibrium. The near-equilibrium 4Z0.5 level curves at tZ100 are shown in figure 10a. The evolution from an initially circular shape to these shapes is similar to that shown in figure 8, but as b decreases sharper and smaller corners form in the equilibrium morphologies. This is further quantified by the plot of the interface orientation angle in figure 10b. Away from the corner, the corresponding three interfaces merge to a shape that matches the sharp-interface Wulff shape. In figure 11, we consider the evolution of a two-dimensional model of a thin film geometry using the new phase-field model. In particular, the evolution of an open, initially flat interface with periodic boundary conditions is shown. The initial geometry is given as 4ðx; yÞ Z 1 2 1Ktanh y K 1:6 2 ffiffi ffi 2 p e :

ð5:2Þ
The domain is UZ(0, 6.4)!(0, 3.2). The fourfold anisotropy parameter is aZ0.5, the regularization parameter is bZ0.01 and the interface thickness parameter is eZ0.03. The mesh parameters are the same as in previous simulations. In figure 11, the 4Z0.5 level curves are shown at early times during the evolution. Initially, the flat interface decomposes into a hill-valley structure as the interface tries to rid itself of high-energy orientations. This corresponds to a rapid decrease in energy as seen in figure 12 and this process is similar to that described earlier for closed shapes. As the system evolves, the hill-valley structure begins to coarsen as hills and valleys are removed from the system. The arrows in figure 11 indicate the sequence in which the hills and valleys are removed. As seen from figure 12b, these coarsening events are associated with rapid drops in the energy. To investigate the coarsening process in more detail, we consider the interface orientation angle q. In figure 13, we observe that the hill-valley structures with the smallest [max(q)Kmin(q)] disappear first, i.e. those with wider corner angles. In particular, the hills marked by 1 at tZ0.5 disappear earlier than the other hills. Next, the hills marked with 2 at time tZ2.0 disappear. At longer times, the coarsening process is shown in figure 14 where the x coordinates of the hills (circles) and valleys (pluses) are shown at different times. As seen in the figure, the coarsening process may occur symmetrically by either two hills merging with a valley to produce a hill or by two valleys merging with a hill to produce a valley. As the interface coarsens, both the hills and the valleys tend to match the regularized Wulff shape (Spencer 2004). This symmetry (e.g. the shape of a hill is the inversion of that of a valley) is seen in figure 15a where the orientation angle for the new phase-field model is shown at a late time tZ350 (solid curve) together with the asymptotic form of the sharp-interface Willmore regularization (dashed curve) from Spencer (2004). Note that three of the hill-valley structures have smaller [max(q)Kmin(q)] values than the others. These are the next to coarsen, although the process takes a very long time. The corresponding interface morphology is shown in figure 15b. We further note that symmetry in both the coarsening process and the shapes of the hills and valleys may be broken when additional effects are considered, such as deposition or elastic forces (e.g. Hausser & Voigt 2005b;Zhou et al. 2008).
We next turn to three dimensions. In figure 16, the evolution of an initially spherical shape towards the Wulff shape is shown, using the new phase-field model with aZ0.3, bZ0.01 and eZ0.03. Observe that the Wulff shape resembles a double-sided pyramid. The computational domain is UZ(0, 3.2)!(0, 3.2)! (0, 3.2). The simulation is performed with a root-level grid of 32!32!32 (with mesh size hZ0.1) and there are two levels of refinement. The resolution on the finest mesh is hZ0.025 and the time step size sZ0.001. As in two dimensions, there are approximately six or seven mesh points across the interface during the simulation. In figure 16a, the 4Z0.5 surface is plotted. In figure 16b, the yZ1.6 slice is shown together with the corresponding slices of the bounding boxes of the three-dimensional adaptive mesh. Each box contains a mesh that is one-half of the size of that in the box that contains it. In figure 17, we compare the yZ1.6 slice (solid curve) at tZ20 with the corresponding slice of the three-dimensional Wulff shape (dashed curve; which also corresponds to the two-dimensional Wulff shape since the slice is through the centre of the three-dimensional shape), and also with the asymptotic Willmore regularized solution (dot-dashed curve; Spencer 2004). As in two dimensions, away from the corner, there is excellent agreement between the Wulff shape and the result from the new phase-field model. Near the corner, there is excellent agreement with the asymptotic solution. In figure 18, we present a final three-dimensional example corresponding to an open, periodic three-dimensional interface, which serves as a model for a thin film. The new phase-field model is used with aZ0.3, bZ0.01 and eZ0.024. The computational domain is UZ(0, 3.2)!(0, 3.2)!(0, 3.2). The simulation is performed with a root-level grid of 16!16!16 (with mesh size hZ0.2). There are three levels of refinement. The finest resolution is hZ0.025 and the time step size sZ0.001. There are approximately six mesh points across the interface during the simulation. In figure 18a, the 4Z0.5 surface is plotted at time tZ40; in figure 18b, the corresponding yZ1.6 slices are shown together with the slices of the bounding boxes of the three-dimensional adaptive mesh. As in two dimensions, the slightly perturbed initial interface rapidly decomposes into a hill-valley structure that coarsens leaving a central smoothed pyramid surrounded by a periodic hill-valley structure. A comparison (not shown) indicates that the yZ1.6 slices of the hills and valleys match the sharp-interface theory of Spencer (2004).

Conclusion
We have presented a new phase-field model for strongly anisotropic crystal growth using regularized, anisotropic Cahn-Hilliard-type equations. A new formulation was used for the anisotropic surface energy and the regularization was performed using a phase-field/diffuse-interface approximation of the Willmore energy (square of the mean curvature). A key feature of our approach was that the interface thickness is independent of orientation. This removed an incompatibility between the anisotropy and the Willmore regularization proposed and used in earlier phase-field models. An asymptotic analysis presented in the electronic supplementary material demonstrates the convergence of the new phase-field model to the sharp-interface model. We presented two-and three-dimensional numerical results using an adaptive, nonlinear multigrid finite-difference method. We found excellent agreement between the new phase-field and the sharp-interface dynamics. The computed equilibrium shapes using the Cahn-Hilliard approach also match a recently developed analytical sharp-interface theory that describes the rounding of the sharp corners by the Willmore regularization.
In the future, we will extend the model described here to simulate heteroepitaxial thin film growth. We will incorporate deposition, misfit strain, a substrate and wetting effects. Through continued improvements in the efficiency of the numerical algorithm, we plan to study the long-time coarsening of thin film systems. In addition, one of the important problems we seek to address is the control of the formation and coarsening of nanoscale structures to achieve desired spatial orderings.