A vertex based approach to crystal facet modelling in phase field

This paper seeks a vertex based approach to faceted anisotropy in phase field modelling of crystal growth. We examine Wulff shapes and the connection they have with phase field formulations. On inspecting current approaches to modelling facets within phase field we observe that there are two distinct approaches: one imple- ments the faceting purely in the kinetic parameter thus avoiding the complications of taking gradients of discontinuous functions; while the second is based upon regularisation of the anisotropy function within the free energy function. Armed with our new insight into the operation of anisotropy within phase field we refine the second of these and advocate a vertex based approach to facet anisotropy modelling. Results include regular and irregular morphologies and hill-valley growth. We also present high undercooling effects on faceting, which can cause breaks in facets and, in the case of irregular shapes, distortion of the underlying Wulff shape.


Introduction
Since the establishment of a thermodynamic basis [1] and the pioneering efforts of [2], the application of phase field to crystal growth has become almost mandatory over the last decades. This is largely because the phase field methodology allows a computationally continuous transition between different phases of matter, avoiding the need for interface tracking, which is problematic for complex crystal morphologies, particularly where the topological connectivity changes, such as during melting or coalescence. The computational price to pay is the increased degrees of freedom needed at the transition, which in 3D is particularly computationally expensive. This led to models such as [3] which exhibit a high degree of phase field width independence. Nonetheless it is vital that the formulation used for the simulations is as efficient as possible and one of the aims of this paper is to improve efficiency in the anisotropy implementation.
The phase field simulation of complex crystal morphologies, uses the phase variable, ϕ, and necessarily involves the inclusion of an anisotropy function, A(∇ϕ), which reflects, at the mesoscopic scale, the underlying crystallography of the ordering at the atomic scale. It was soon realised that by varying the interface width as a function of the normal to the surface, that anisotropy can be introduced into the model [4], though this may now change with the introduction of models that keep the interface width constant [33]. With anisotropy the modelling of conceptually simple morphologies, such as dendrites, began in earnest, initially in 2D. However, the practical implementation of anisotropy in phase field involves reducing A to a somewhat less elegant non-linear function of gradients of the phase field, ϕ , which upon brute force variational computations can become an unrecognisable jumble of partial derivatives. This is particularly the case for complex morphologies, such as crystals faceting on multiple non-orthogonal growth planes in 3D. Consequently, the simulation of faceted morphologies is much less common in phase field than is that of morphologies with a continuous interface, such as dendrites, and mixed growth where continuous and faceted phases co-exist has hardly been addressed at all.
The interplay between co-existing faceted and continuous crystals is a key feature in explaining both the performance and potential failure of numerous advanced engineering alloys, particularly as many intermetallic compounds form as faceted phases. Examples of the deliberate inclusion of faceted phases in continuous matrix materials range from the common place, such as most Al − Si casting alloys to the exotic such as Re-rich hexagonal δ-phase in Ru containing high temperature superalloys [5]and the non-faceted/faceted Mg-Mg 2 Ni eutectic proposed as a potential hydrogen storage material [6]. Similarly, the formation of unwanted faceted phases in continuous matrix alloys can have a severe deleterious impact on performance, with for instance iron contamination causing widespread problems across a range of commercial Al and Al-Si alloys with the formation, among others, of needle Al 13 Fe 4 , Chinese script α-Al 8 Fe 2 Si and tetragonal plate β-Al 4 FeSi 2 [7].
The over-arching principle behind phase field modelling is the minimising of free energy, subject to diffusion parameters such as solute diffusion and mobility. Much of the modelling is within the formation of the free energy, which splits neatly into two parts: the surface energy; and the bulk energy. The latter, bulk modelling, is the more closely linked to thermodynamics, whereas surface modelling is completely divorced from bulk properties. This paper concerns surface modelling. Surface modelling and phase field appear simultaneously through the introduction of an interface width, being a measure of the size of the transition region. The most ubiquitous and extensively explored model is the two-dimensional 4fold anisotropy: where, A is the anisotropy and ν is the angle the normal to the solid--liquid surface makes with some fixed direction. Typically the horizontal is used as a reference direction, but multiple reference directions may be used for different particles. e.g. [8]. A phase field parameter, ϕ ∈ [0,1], is utilised, where ϕ = 1 2 represents the heart of the transition region and thus, formally, the surface. The normal direction is given by ∇ϕ which when normalised is n = ∇ϕ/|∇ϕ| and the link with the angle, ν, is given by where ϕ x ≡ ∂ϕ ∂x etc, for Cartesian coordinates, x,y. A major contribution to the understanding of anisotropy in phase field is presented in [9] which associates a vector field (termed the ξ-vector) with the anisotropy and the resulting equilibrium crystal shape -the Wulff shape: [9] also develops, to some extent, the method for implementing this in phase field. We have a further look at this term from the perspective of a differential form in this paper.
We then inspect facet models. A feature of all facet models is to create anisotropy with forbidden angles. To understand this requires knowledge of Wulff shapes and their relation to anisotropy. For the dendritic models the anisotropy (of which A = 1 +∊cos4ν is an example) is such that A ′ (ν) is always a continuous function of ν. However, as the magnitude of the anisotropy, ∊, is increased, the Wulff shape takes on a retrograde motion with magnitude A ′ (ν) and develops features known as "ears", which are in effect discontinuities. The onset of these discontinuities is when A + A νν = 0. Now, since the Wulff shape is defined as not including the ears, increasing ∊ without bound can give increasingly faceted crystals. The problem though is that the highly non-regularised (discontinuous) anisotropy function would be very unstable and problematic if kept in the free energy, and its adoption in the kinetic term became common, e.g. [10][11][12][13][14]. Some modellers stayed loyal to free energy implementation [15][16][17][18][19], through the approach of regularisation. The first two papers, produce faceting by increasing the strength of the anisotropy to produce a Wulff shape with pronounced cusps (ears), whereas the latter three papers work with anisotropy with formally flat sides, no cusps and rounded 'vertices'. The example in [17] is a square and that of [18] is focussed around hexagonal growth, but in [19] is a more general methodology that includes arbitrary shapes in 2 and 3 dimensions.
The approach of [19] is focussed on the facets of the shape in 3D. Once the normals to the facets are specified, then a surface energy can be formed that is a function of these normals and the direction of growth. This approach is to be compared with the approach advocated in this paper, where it is the vertices of the polyhedra or polygons that form the basis of the surface energy construction.
Facet models in phase field tend to begin with an "ideal" anisotropy function that typically contain unappealing discontinuities which, in turn, will likely cause numerical instabilitysee [17]. The anisotropy function with 4-fold symmetry in 2D, given in [17], is where, one notes the absolute value function. If we rescale this as and increase ∊ from zero, the associated Wulff shape moves progressively: from a circle; to a square with rounded corners; to a perfect square. Thus, one is led to conjecture that the anisotropy of a perfect square is formally Fig. 1. The anisotropy of a square with ∊ = 0 (isotropic) and as ∊ = 1, 2, 3 the anisotropy approaches that of being related to a perfect square Wulff shape (not shown).  2. The corresponding Wulff shape to Fig. 1 of a square with ∊ = 0 (isotropic) and as ∊ = 1, 2, 3 the anisotropy approaches that of being related to a perfect square Wulff shape (not shown).
A plot of Eq. 4 for four equally spaced values of ∊ is shown in Fig. 1 with its companion plot of Wulff shapes in Fig. 2. A closely related function that also gives 4-fold faceted symmetry is This has a similar Wulff shape, but what is noticed on closer inspection of this function is that there is a critical value, ∊ = 1/3, where the vertex becomes a true point, and the adjacent sides are not straight near the vertices. The anisotropy function Eq. 4 is thus preferential to Eq. 6. This raises the question: what properties of the anisotropic functions give clear faceted Wulff shapes with sharp (point) vertices, straight sides and nothing else? Developing phase-field models for sharp corners and facets is challenging because sharp corners are necessarily associated with negative stiffness and facets with a Wulff shape with missing orientations (A ′ (n) undefined). This leads to modifications of the model by rounding vertices, [17], and convexification, [20], and also by adding higher order terms, [21]. Higher order terms considerably increase the numerical difficulty and so we explore a convexification approach in this paper along with the rounding the vertices. We argue that if the equilibrium (Wulff) shape is identical for two candidate schemes, then the metastable results must also be identical, because the Wulff shape is the input to the dynamic phase field model. It follows therefore that a canonical formulation is the most natural starting point of the different formulations, leading to the most straightforward geometric interpretation. Ultimately, it must be modified for direct implementation, the end result being a formally convex Wulff shape with approximately sharp corners and near flat edges (in 2D). An ideal phase field formulation would perhaps include only measurable parameters, and so, in principle, one would acquire surface energy values from, for example, experiment or molecular dynamics simulation. The approach taken here is to inspect equilibrium morphology from observation and infer the surface energy from the vertices of the faceted object using the approach we detail later in this paper. In principle, when the equilibrium shape is reproduced in the simulation, metastable shapes will be explored at higher driving forces, possibly with the aid of anisotropic kinetic parameters, in order to make predictions of the non-equilibrium shapes found in experiment.
In this paper we address the above questions via a discussion of Wulff shapes in Section 2, facet modelling of simple shapes in Section 3 and more generally in Section 4. We give a short section on basic differential geometry notation, which we use in parts of the paper, in appendix A; details of the phase field model that we use for an inter-metallic alloy are provided in appendix B; and a discussion of the numerical method we adopt is in appendix C.

Anisotropy in phase field
Part of the material of this section is also discussed in different notation in [9] which applies the Cahn-Hoffmnan ξ vector [22] to phase field. The ξ vector is defined as the gradient of the radial coordinate times the anisotropy ξ = ∇(RA) but we prefer to use the more flexible and primitive concept of a differential 1-form The relationship between the two is where G ij are the components of the metric of flat space. The most natural coordinates to use are Cartesian when the metric matrix G ij = δ ij is diagonal giving However, for polar coordinates the components of the differential 1form, W and the ξ vector do not equate unless presented in an orthonormal basis.
We should point out that the more familiar differential operator, d (italisized) is treated and thought of as a scalar in normal vector calculus, and means d(RA) ≡ ∇(RA)⋅dr = ξ⋅dr (11) where r is the vector position,[x,y]. Hence, though all intimately related, the scalar differential, d(RA), is formally a different object to both the vector ∇(RA) and the one-form, W = d(RA). Please see the first of our appendices in Section A, which explains in more detail the geometric significance of a 1-form field and other related notation used throughout this paper. The following subsections illustrate and explicate the concept of a Wulff shape and its relation to the anisotropy -Section 2.1, and then how this mathematical object relates to a phase field formulation -Section 2.2.

Wulff shapes
To understand the connection of Wulff shapes with the anisotropy, A, it turns out that we require the function where ν is the angle of the vector X = R[cosν, sinν] to the horizontal so that R = |X| is the length of the vector. The angle ν is not to be confused with the angle, θ, at the origin where the vector X is located -see Fig. 3 Given, f, one can generate a Wulff shape by where the one-form (gradient) is interpreted as lying in the 2D domain centred at the origin for a particular R, say, R = 1. The simplest example is the isotropic case, f = R so that with, X = Rcosν, Y = Rcosν which is a circle. More general anisotropy such as generates the Wulff shape via where e R = dR = cosνi +sinνj and e ν = R dν = − sinνi +cosνj are unit vectors in the radial and azimuthal direction respectively (see Fig. 3).
We are also at liberty to use Cartesian coordinates, (X, Y) in our evaluation of the Wulff shape which gives This means we can equate ∂f ∂X The Wulff shape is given by the parametric curve using the Cartesian components of df: This also makes clear that the coordinates of the Wulff shape, (x, y) are distinct from the components of the normal n = (X/R, Y/R).

Anisotropy and Wulff shapes in phase field
On implementation into phase field the anisotropy enters the energy functional via the relations X = Rcosν, Y = Rsinν and the function For example, (20) gives f defined in terms of X, Y as . The surface energy functional in phase field is given as Applying standard variational techniques on F gives ∂f ∂∇ϕ (22) The question arises: What is the relation between the Wulff shape and the phase field formulation? To answer we compare the Cartesian components, in 2D, of Eq. 18, with the components of ∂f ∂∇ϕ , with the identification From now on we shall use the concise notation We finish this subsection with an evaluation of the expression: where we have in mind The most straight forward evaluation of this expression is to use Cartesian components and the chain rule to write We can write this without components as ∇⋅Dg = ∇∇ϕ : DDg. where It is convenient from this point on to use the symbol f for both AR and f ∇ϕ where the context is unambiguous. We can write Eq. 29 in the following way to make the geometric connection with Wulff shape clear by using the notation of [9], ξ ≡ Df, so that which duplicates a result found in [9] and where clearly these results apply in three dimensions (and higher). We make particular use of the well known expression (see [9] for example), Eq. 27, in the facet model we advocate in this paper. There is much discussion in [23], of the geometry of Wulff shapes and a number of models, including many in this paper, that are realised using level set methods. [20] show that the Wulff shape can be defined by for all n and use this to prove results in phase field including the behaviour of a facet model at a vertex. It is one of the purposes of this paper to verify whether equilibrium results for facet models do indeed coincide with the faceted Wulff shapes.

Facet modelling
This section states and looks at the consequences of a definition of anisotropy relating to surface orientation, n and crystal vertices, p i , given in Eq. 33, which is the central equation of this paper. But in order to make a connection with a particular model in the literature due to [17] we work out in detail two equivalent applications of this model to this simpler square case in Section 3.1. Other approaches can be found, for example, see that advocated in [25,26] which use continuous polynomial functions of angle with powers typically in excess of 50: where p i specifies the vertices and the dot product is associated with an angle, and Hi() is the Heaviside function. Eq. 32 has many advantages, the directions of minimum energies and their values can be controlled, as well as the metastability and instability of the missing orientations.
We will now establish a connection between this model and the one we develop in this section, which we state here to be: for n vertices. Now the max function of a sequence of values, also known as the infinity norm, is well approximated by for large w and q i > 0. So we see that the expression Eq. 33 would be similar to Eq. 32, where a 0 = 0, a i = 1 if a power of 1/w would be included. Although expression Eq. 32 has the advantage over Eq. 33 of being continuous in all its derivatives, we find that Eq. 33 is ultimately more simple in application.
Another modelling technique might be termed large ∊ modelling whereby the Wulff shape has large ears leaving the desired crystal shape inside the envelope. In all implementations the resulting Wulff shape with ears is regularised. Examples of this modelling technique are found in [27][28][29][30]. The application in [28] is to that of snow crystal and a variety of crystal structures are produced by varying 4 model parameters. However, the model suffers from the disadvantage of not being clearly derived from a variational of the free energy and contains, for example: an elongation parameter, Γ that appears equivalent to a post processing of the crystal by a transformation in the z direction; also present in the model is an anisotropy function multiplying the bulk driving force. Despite these caveats the range of snow crystal morphologies produced is potentially of great interest to metallic solidification (with theoretical modification if necessary).
An important contribution to the literature is found in [31]. Here a level set method is used to produce 2D and 3D shapes including a dodecahedron. The method is to specify the geometry using angles using both large ∊ and absolute value functions. A discussion and classification of different approaches to faceted equilibrium crystal models is found in [20], but in doing so does not single out any one approach as advantageous or more mathematically or physically justified.

Square crystal growth
The purpose of this section is to examine probably the simplest case of faceting, the square crystal, from two different but equivalent perspectives: the angle coordinate and the Cartesian coordinate approaches. Traditional anisotropy is formulated in angle coordinates and thus it would seem natural to use this approach for faceting too. However, we will see that working with Cartesian coordinates is much more natural and offers simpler formulation. Because the square Wulff shape is relatively simple it allows special treatment and leads to a relatively easy form for implementation in phase field. This also allows comparison with the work of [17] and also serves to illustrate why applying either coordinate approaches to arbitrary faceted crystals in 3D will likely lead to difficulties upon implementation.

Using Cartesian coordinates
We define the surface energy term g = 1 2 f 2 using the anisotropy (strictly speaking f = RA, where A is the anisotropy, but there seems to be no name for the function, f) and we define f for a square: where X, Y are Cartesian coordinates. We are interested in evaluating terms like and where so that which can be written For the repeated derivative in the same direction we obtain ∂ ∂X which is not defined at ∊ = 0. Hence, we choose an approximation to this by assuming ∊ small. See Fig. 4 for how the function is approximated. We also exploit the symmetry of the shape to work in the upper right quadrant only to give the final result for ∇∇g as The final expression for the gradient term is As a footnote to the above it is interesting to create a model without the Dirac function for comparison. This also simulates a pseudo square but with less straight facet sides and more rounded corners. This pseudo square model, does not appear to have a corresponding free energy function. That said, it is a useful first step to ignore the Dirac function on implementation for any desired crystal shape.
Having considered the use of Cartesian coordinates we now look at the equivalent representation based upon an angle coordinate. We now look at an alternative way of incorporating this model that we can relate directly to the formulation found in [17]. We write and then introduce some regularization to write where the effect of a chosen constant, λ, is to round off the corners of the square: λ = 0 is a circle; λ = ∞ is a square; intermediate values give a square with rounded corners. Using Eq. 48 we require (A 2 ) νν . Considering the positive (X > 0, Y > 0) quadrant and its neighbouring quadrants and boundaries: we find that Now focus just on the discontinuity at ν = π/2 and write and then approximate using a small angle δν≪1 so that for angles near π/2 which is less than unity. This then gives Finally we require away from the boundary and for π/2 − δν < ν < π/2 +δν we obtain Using the expression Eq. 55 we can work backwards to give the effective, A that Eq. 55 implies: This is illustrated in Fig. 5 which shows the original, A as a solid red line and the approximation using a dotted blue line. We also show the function A +A νν for this same approximating function, and clearly as δν→0 we would develop infinite spikes. Clearly the smaller the chosen value of δν or ∊ the better the Dirac delta is approximated and an upper limit is the angle between the centre of the facet side and the nearest vertex; in this case, ϕ/4. We find δν ≈ 0.01 is stable, but have not explored this issue beyond choosing a reasonably small value that does not cause undue stability and performance compromise. Both A = |cosν| +|sinν| and Eq. 55 can be used as demonstrated by [17]. However, it is much easier using the Cartesian computation, Eq. 47 plus a necessary provision for the discontinuity as we shall see in Section 4.
This section has discussed two alternative approaches for working with a square Wulff shape. The Cartesian approach seems slightly simpler, but both approaches are likely to be difficult to apply for general shapes in 3D. The difficulty being that the location of the discontinuities along faces and edges can be a substantial geometric undertaking for a general collection of vertices.

The vertex based facet model
In this section we look afresh at the anisotropy function for faceting, taking into account that the model of [17] is an example of the correct approach applied to a square geometry. The role of ∊ in the model Eq. 4 is to control the curvature at the vertices and thereby provide more stability upon implementation.
The Wulff shape associated with anisotropy, A, in 2D, given by where we use the notation, W rather than ξ because we are thinking of the whole curve with e R = [cosν, sinν], e ν = [ − sinν, cosν]. For faceted morphologies the curve is not continuous. For example, the anisotropy A = |cosν| +|sinν| has a Wulff shape consisting of vertices of a square with no joining edges. Another function that might be used to generate a square is the function A = 1 + ∊|sin2ν|. For values of 0 < ∊ < 1/3 a faceted shape emerges but it does not have the property, for any ∊, that the Wulff shape only consists of vertices. So, in some sense, the formulation of [17] where A = |cosν| +|sinν| is more true to the Wulff shape than, for example, A = 1 + ∊|sin2ν|. More generally, for any desired Wulff shape: what is the property of the anisotropy that gives precisely the vertices of the desired Wulff shape? For this we require all angles from the centre of one facet to the centre of the next facet to map to the one vertex common to both under the mapping Eq. 59. The following geometric figure for producing a regular hexagon, Fig. 6, illustrates the property we require. The solid blue circle arc in Fig. 6 is given by giving Hence, for − π/6 < ν < π/6, we have That is, all angles in this range map to one common point, [1,0]. Intuitively, as the angle, ν swings round between 0 and π/6 the vector, A ν e ν grows in such a way as to keep W at one point. We have the added bonus in this model that the angle formulation is readily and directly written in terms of Cartesian components of ∇ϕ (we shall see shortly why this is In 2D this implies that and in 3D ∇⋅D The above result is not complete because of the discontinuities in the first derivatives at X = 0 and Y = 0. The final result has been given previously in Section 3.1, but we also cover the more general case in the next sections. What we have achieved in this section is a description of a direct and simple route from an arbitrary faceted shape defined by its vertices to the surface energy, which we now illustrate. Figs. 7-12 show some hexagonal Wulff shapes, beginning with a regular hexagon in Fig. 7, and    an irregular hexagon in Fig. 8. In both cases the vertices lie on a common circle (centred at the origin). This has the consequence that the anisotropy (shown in dashed blue) meets at the centre points of each facet. In the phase field simulation of faceted Wulff shapes it is often advantageous to regularise these shapes to include rounded vertices as in Fig. 9. Another irregularity that can occur is when the vertices are not equal distance from a common point such as in Fig. 10. Here the anisotropy still forms a series of circle arcs that meet on a facet, but in these cases they do not meet at the middle of the faceted side. The extreme event where the anisotropy fails to meet on the facet is not allowed (see Fig. 11), because the two surface energies no longer have a common value on the facet. On the other hand if we place one of the vertices too close to the origin it will simply be ignored as in Fig. 12.
Despite the description of surface energy being simple, 3D phase  Fig. 12. When a vertex is placed so as to produce a concave shape, it is simply ignored. For example, a vertex ( − 0.2, 0.2) has no effect on the anisotropy, and, by implication, on the surface energy. There is no concave Wulff shape.  field implementation for arbitrary shapes is potentially difficult. We now discuss a straightforward way of circumnavigating the potential difficulties.

Arbitrary vertices and implementation
The generalisation to arbitrary irregular convex polygons as in Fig. 8 is readily achieved (although the physical basis for such irregularity may be questioned). But we ultimately seek a model construction relating to crystal growth in 3D, which in its generality will also relate to 2D crystals.
Our starting point is Eq. 33, which gives: where the application of this expression for arbitrary (and acceptable) vertices is non-trivial because of the discontinuities in the derivatives for certain values of X, i.e. normal to a facet. See Fig. 13 for the anisotropy of a cube. The versatility of expression Eq. 33 is illustrated in Fig. 14 for a dodecahedron and in Fig. 15 for a skewed rectangular prism. Of particular note in Fig. 15, is that the anisotropy of a skewed cube is not the skewed anisotropy of a cube. The simplest way to see that this is true is to use the knowledge that all anisotropies of 3D faceted Wulff shapes consist of interlocking spheresnot ellipses, for example.
We require ∇∇g, where g = 1 2 f 2 . In the neighbourhood of vertex p i our function and so (in this neighbourhood) and when X points in the direction of p i , or even in its vicinity, i.e. closer than neighbouring vertices.
We now concentrate on two adjacent points, p i and p j , joined by a line and separated by a plane surface, S, such that the position X on the plane obeys: That is n ij ≡ (p i − p j )/|p i − p j | is orthogonal to the plane. We state as a postulate that: 1. The solution to Eq. 73 must include a vector X pointing between p i and p j .
2. The solution to Eq. 73 must include a vector X normal to, and located on, a facet.
These statements have the effect of restricting the choice of vertices to not only form a convex shape but also arranges the facets to be normal to some radial vector from the origin.
Consider a function where the whole domain is divided into two regions A and B separated by a surface, S and furthermore on the joining region. Consequently This situation arises in the description of faceted anisotropy. For regular shapes, one can use lines and surfaces of symmetry to address the discontinuities that exist on these surfaces. But the complications accumulate in 3D even for regular shapes leaving this approach only useful for simple shapes in 2D like squares and hexagons. In fact, even in the simplest case of a square the problem is not trivial due to the presence of discontinuous derivatives as [17] demonstrates.
We seek an easy-to-implement method where the input are the vertices, p i ,i ∈ [1, n], of the desired Wulff shape only. The essence of the approach is: once f is defined as a maximum of a series of n values (equal to the number of vertices) for any given point in the tangent space, X = ∇ϕ, one avoids seeking an analytical expression for its derivative and, instead, uses a finite difference approximation. Moreover, one can avoid finding the second derivative altogether by using a finite volume approach in a finite difference code, as we will describe. An alternative, which also avoids a second derivative, but is not explored here, is the Finite Element Method, The finite difference approximation is not directly associated with the underlying mesh because we are seeking a derivative on the tangent space, i.e. on the space of X. The use of a small number, h on this space, gives an exact result for df when X points towards or near a vertex, p i , because f is linear in X,Y. However, when X points midway between two vertices (or more), the finite value of h has the effect of smoothing the discontinuity. One is of course free to choose smaller and smaller values of h→0, but ultimately this will lead to numerical instability. The penalty for too large an h is that the facets will become curved.
We first explore the simplest implementation, but then follow this with a more stable finite volume approach.

Implementation in Phase field
We give full details of our phase field model in Appendix B, where we adopt the Allen-Cahn equation for the phase field and Cahn-Hilliard for solute diffusion, following minimisation of the free energy. The free energy consists of bulk and surface energy, the latter of which controls the Wulff shape via the anisotropy, A(n).
For faceted shapes the anisotropy, A(n), leads to ill-posed equations, [21], and regularisation, for example, by adding a curvature term to the energy is demonstrated in [21,32]. Building on these approaches [33] establish the connection between sharp interface theory and phase field, and demonstrate the coarsening dynamics of the faceting instability during growth. The main advantage of these approaches is that it allows the phase field model to incorporate metastable anisotropy, i.e. that leads to Wulff shapes with 'ears' (also referred to as 'wings' in [21]). Arguably, a simpler approach is given in [17], which retains the Allen-Cahn phase field equation for the simulation of a square Wulff shape, and is the approach we build on in this paper for general shapes in two and three dimensions. A feature of this latter approach is to avoid a Wulff shape with ears by recognising that the flat faceted faces (edges in 2D) can be approximated by faces with near zero curvature and using (rounded) 'vertices'. The resulting Wulff shape being formally convex to allow stability in the computation. We can associate the model with a tensor as defined in [24] (using only the gradient contributions): which, comparing with the expression (5.10) in [21], reveals it to be the simpler model. For ease of reference we rewrite our surface energy function where, as in Eq. 4, ∊ controls the rounding of corners to aid numerical stability as noted by [17]. A finite difference approximation to a second derivative in the X direction is given by Now, in the continuous part of the domain, away from the forbidden directions, this result is exact for quadratic functions of X, as we have here. For X lying on a discontinuity, the Dirac function is essentially approximated using the chosen value of h, and clearly as h→0 the true value is approached. Similarly The above certainly has the virtue of simplicity and is easily imple-mented in both 2D and 3D. That said, despite ∇∇g being approximated, it is still an approximation to a singularity. So we detail an alternative approach.

Regularization and implementation
This section details a finite volume approach to finite difference using our knowledge of the Wulff shape and its relation to phase field. The approximation to the singularity is avoided but the problem of h→0 is not avoided. This is our preferred numerical method which we detail for 2D only, but also a finite element method will have the same advantages.
1. We wish to implement, using finite differences, the gradient energy where for any X = ∇ϕ, f satisfies for the n vertices, where n is the normal to the boundary of volume, V, denoted ∂V. In 2D on a Cartesian mesh labelled i, j this becomes four integral contributions for each side of the square surrounding the point (see Fig. 16), one of which is  On the right-hand side of Eq. 88 is a function of X = ϕ x and Y = ϕ y which needs to be a evaluated at the point (i + 1 2 ,j + 1 2 ). On a regular Cartesian grid, with Δx = Δy, we use: The point (i + 1 2 , j − 1 2 ) is dealt with in a similar way. At the point 3. Now a crucial key to the method is the evaluation of the partial derivatives ∂f ∂X , ∂f ∂Y , which are not evaluated analytically. We use the stencil Please note the distinction between h and Δx. Δx is the grid distance between mesh points (assumed equal in x and y) but h, related to a typical maximum value of |∇ϕ|, can be chosen arbitrarily taking into account machine accuracy and how closely one wants the numerical scheme to follow the approximation to the max function. A smaller h gives lower curvature on the facets. For example, for a grid size of Δx = 0.78 and maximum |∇ϕ| ≈ 0.1/Δx, we explored values of h ∈ [0.001, 0.02], being roughly in the range of 1/5 to 1/100 the max |∇ϕ|. The larger h giving noticeably curved sides, and the former smaller value reducing the stability so that the simulation is an order of magnitude slower to run. It is beyond the scope of this paper to quantify numerical stability, but h does seem to result in curvature, κ̃h/max(|∇ϕ|).

The final scheme is then
assuming a square cell of size Δx 2 . Fig. 16 illustrates the scheme with the part boundary, C Right shown as a dashed line.

Recommended scheme: approximating the max function
We now explore, use for our results and recommend, the use of a power law approximation to the max function. That is, the approximation of the max function using higher order polynomials where we assume that |p i | = 1 for all the vertices. Use of this continuous function avoids the difficulty, present in the previous scheme, of having to have a large, h, on the tangent space. This is because, in the previous scheme, the parameter, h, served two purposes: one, a convenient way of implementing arbitrary vertices, and; two, a way of approximating the max function. In Eq. 93 the approximation to the max function is already made, and so the use of a finite difference, h, only serves for the convenience of dealing with arbitrary vertices. Crucially, h > 0 can be made arbitrarily small (and is not restricted by the grid size). A note on implementing Eq. 93: for w = 512, for example, we find it far more efficient to use ((x 8 ) 8 ) 8 rather than x 512 or exp(wlnx) in the Linux 17.0 Intel FORTRAN we employ. It may be argued that the final scheme we recommend is very similar to the power law schemes of [25,26]. The difference, though, is that our proposed scheme incorporates an approximation to the max function, and does not just resemble it. The main ingredient common to all facet modelling is to generate discontinuities in the first derivative in the direction of the tangent to each facet.
We have given a method for generating surface anisotropy given the vertices of an equilibrium shape, but what of smooth shapes that cannot be specified by vertices? To some extent this question is beyond the scope of this paper, but since any surface can be approximated by a series of smaller and smaller facets, any convex shape can be approximated with this method. For example, in Eq. 33, set all vertices, p i , i = 1..N, to lie equidistant on the unit circle. Then is follows that, as N→∞ the maximum value is attained when p i = n and so A = 1 and the isotropic model is reproduced.

Phase field simulation results
Motivated by curiosity to see the effect of the discontinuity in the anisotropy model, we first inspect the pseudo model given in Eq. 47, where all the complications of the discontinuity in the anisotropy are ignored. In Fig. 17, we show one of the effects of neglecting this term. Here, the curvature, κ, is output at a point in the middle of a facet of the desired square Wulff shape, with and without the correct treatment.
Comparing the results, we see that the facet is more quickly achieved and remains flat with the discontinuity included via the approximation of the Dirac delta function. In short, the correct treatment of the discontinuities aids in establishing the Wulff shape quickly and Fig. 17. Two plots of the evolution of curvature, κ, on a facet side, with and without the correct treatment of the discontinuity. R is a measure of a increasing size of the crystal. The first derivative of the anisotropy is discontinuous and so the second derivative, needed in phase field, formally produces a singular (Dirac) function. The plot that lies predominantly at κ = 0 (blue) contains the correct treatment of the discontinuity. Neglecting this term can lead to concave (κ < 0) sides and surfaces as can be seen as the (red) line crosses κ = 0. maintaining its shape.
Unlike rapid solidification of dendritic growth, near equilibrium growth of a regular faceted crystal tends towards the Wulff shape. In practice "near equilibrium", in our context, is when the bulk driving force is sufficient for growth. Simply setting this term to zero results in melting and so, by trial and error we find a temperature just low enough to give growth.
However, its evolution towards that end state is interesting in itself and can be compared with non-phase field models such as [34,35], being discrete models for growth and dissolution respectively. As is normal practice, the phase field simulation seeds the growth in a metastable energy state as a circle or sphere.
We begin by considering grid-induced effects on a growth of a regular hexagon in Figs. 18 and 19 which is coloured according to values of the solute field, c, with the internal region being at the solid concentration. We employ the approximation, Eq. 34 with w = 512 to Eq. 33 in Eq. 79 and the implementation methods of Section 4.1. We use ∊ = 0 for Fig. 18 and ∊ = 0.1 for Fig. 19. Part of the quarter domain grid used is shown to illustrate that the hexagonal crystal has been rotated by π/3 to aid the eye and demonstrate that the larger value of ∊ effectively eliminates grid anisotropy effects. Not shown in the plots, and tangential to the main argument of this paper, is the observation that partitioning is that the maximum value of c is found in the centre of each facet and approaches the equilibrium value for c L as given by the common tangent construction to the free energy curves used. From this point on all simulation results concentrate on the morphology. Fig. 30 is a superposition of a series of snapshots of crystal growth starting from eight seeds of varying sizes. The results are similar to those reported in [33] and, more uniformly, in [32]. The most notable aspect of these simulations is the formation of stable concave sharp valleys from a convex Wulff shape.
It is also instructive to consider shapes without regular Wulff shapes to see the evolution of the faceted morphology towards or away from the Wulff shape. An example of this is the Wulff shape in Fig. 8 using a near equilibrium phase field simulation at a series of times in Fig. 20. However, when the undercooling is increased by 100 K or more departures Superimposed is the Wulff shape as a dashed line, with dotted lines from the origin to aid the eye. To get the best fit the centre of the Wulff shape is located slightly above the centre of the phase field initial condition.   from the Wulff shape are evident-see Fig. 21. We also notice in Fig. 21 that there is a delay before the final shape is formed.
By adopting the WBM model, [2], with inter-metallic Al 13 Fe 4 we found, even with a high undercooling ΔT/T M ≈ 1/6, the resulting morphology returned a scaled version of the equilibrium shape, e.g. the inner regular hexagon shown in Fig. 22. Changing the model to that of KKS (see [36] and appendix B) we found that we were able to generate greater growth rate and consequently were able to achieve the facet breaking reported in [17,18]. Extending the run time results in Fig. 23, which has established a dendritic structure with the facet orientations no longer in the allowed (equilibrium) direction of a regular hexagon.
In summary to get a true Wulff shape in a phase field simulation we require near equilibrium and a high order polynomial approximation to the max function, Eq. 33.
In 3D, when changing the anisotropy to a rectangular prism facet model with ratio 1:1:2 the initial spherical seed evolves towards the 1:1:2 state quite slowly, arriving at   Simulation of (equilibrium) hexagonal shape using two different models: WBM and KKS. Both model results are shown at the same time step indicating a higher growth velocity for the outer shape. Using WBM method the growth velocity is low and the equilibrium shape, the hexagon is reproduced. For the KKS method a greater growth velocity results in facet breaking. Fig. 24. Phase field simulation of a rectangular prism. the facets and is intended as a visual aid. Our final, 3D simulation, which illustrates the flexibility of our proposed method is of the Archimedean truncated octahedron shown in Fig. 29. Here again, as an aid to the eye we have coloured the faces according to their orientation. The faces in the figure contain both regular hexagons and squares, and although it is indeed the equilibrium Wulff shape desired, it is not a forgone conclusion that such a result should emerge intact without distortion by the phase field method. Fig. 29 thus confirms the method works, despite the attendant regularisation, smoothing of the interface and rounding of the vertices and edges.
The main aim of this section has been to illustrate and confirm the generic nature of our proposed method to accommodate any Wulff shape, in 2D or 3D, within phase field. We have also illustrated the importance of the regularisation introduced in [17] to combat grid anisotropy by 'rounding' the vertices. We have also illustrated the applicability of the method to irregular shapes with less symmetry and to, to a limited extent, the effects of non-equilibrium dynamics on the resulting departures from the underlying equilibrium shape.

Summary
We have presented a very general method for the phase-field simulation of arbitrary faceted crystals in both 2-and 3D, illustrated for a range of both regular and irregular morphologies. The method lends itself both to a relatively straightforward geometrical interpretation of the simulated morphology, via the Wulff shape, and is of comparable ease to implement for arbitrary vertices as, for example, [25,26], but offers a more canonical approach.
The paper is concerned with a new approach to facet modelling and, as a result, claim a phase field method which is both flexible for arbitrary    faceted shapes and simple to implement. As such, we believe this approach offers a radically new way of simulating realistic morphologies for complex, faceted crystals in phase-field. Moreover, due to the simple way in which the crystal morphology is specified in the model, i.e. by the specification in real-space of the co-ordinates of the vertices, experimentally determined crystal shapes may be easily parameterised for simulation.
The essence of the approach is contained in Eq. 33, which has the flexibility that the input to the model are the vertices of the Wulff shape only. The main results of the paper are: the facet anisotropy Eq. 33 and its approximation, Eq. 34; and the implementation methods given in Section 4. Also, to mitigate against grid anisotropy and help with numerical stability, we advocate the application of a term originating in [17] to give rounded corners in Eq. 4, and generalised in Eq. 79 as illustrated in Fig. 19.

Data availability
The authors confirm that the data supporting the findings of this study are available within the article [and/or] its supplementary materials.

Declaration of Competing Interest
The authors declare that they have no known competing financial interests or personal relationships that could have appeared to influence the work reported in this paper.
6. Now consider the differential, df df : By this definition, when f is a coordinate, say, f = x i , we have Consequently, multiplying both sides by ∂f ∂x i and summing produces which implies, formally, that Thus definition Eq. A7 is consistent with the standard definition of a differential. 7. When v is a basis vector, v = ∂ ∂x j we have the result that where δ i j takes the value 1 when i = j, otherwise zero. 8. Now it also follows from the definition in Eq. A7 that the operation of df on vectors is linear: for any function λ, which means that df is a first rank tensor with basis in R 2 dx 1 , dx 2 (A13) 9. We can define any dual vector, α, in this basis by and its operation on vectors by hold. This is illustrated in Fig. 31 where, assuming the contour are for heights in unit steps then df(v) ≈ − 1.7

Appendix B. Phase field model
Throughout this paper we employ the following non-dimensional model for the phase field: where on a domain, Ω, Here, δ is an input parameter controlling the interface width, λ is the dimensionless undercooling, where T M is the melting temperature, and Δf is to be explained in Section B.2, and effectively relates latent heat of fusion, L, by L = Δf/λ. f B is the bulk free energy associated with an intermetallic model for Al 13 Fe 4 in this case as this alloy exhibits faceted crystalline structures.

B.1. Bulk thermodynamics
We begin with free energy of the liquid and solid  (B7) The coefficients vary depending on the temperature, T, and we use T = 1100K and the arrays of coefficients: and Gibbs functions The above coefficients come from CALPHAD data base, [38].

B.2. Quadratic common tangent approximation
The bulk free energy is formally contructed but we recognise the difficulty of working in phase field with a function that is only defined in the range c ∈ [0.2350, 0.27246]. Moreover, the above definitions add significantly to computation. So we adopt an approximation to the above by constructing two quadratic functions, f L ,f S that agree at the common tangent points in the following manner: We define Δf from the common tangent construction, being the difference in values of f B at the two common tangent points, c L , c S , i.e.
Use of a quadratic approximation eases the application of the KKS model of [36], which we employ in producing the results depicted in Figs. 22 and 23. In essence the bulk free energy is formed in a similar way to Eq. B15 but with the following difference f B (ϕ, c) = g(ϕ)f L (c L (c, ϕ)) + (1 − g(ϕ))f S (c S (c, ϕ)) (B17) supplemented by the condition For quadratic bulk free energy curves Eq. B18 is linear and so, after substitution, Eq. B17 forms a non-linear equation for, say, c S (c) (or c L (c)) given any value of ϕ and c, where it should be noted that, in the KKS model, c S and c L are fields. Further discussion of this model is beyond the scope of this paper.

Appendix C. Numerical method
The numerical method is described in detail in [39], but with one significant change, which increases the flexibility of the code. The non-linear, adaptive mesh, implicit, multigrid solver describe in [39] uses a pointwise Jacobian to solve the system. The construction of the Jacobian in [39] approximates the gradient term by ignoring the anisotropy. This treatment is modified in the present code as follows. Consider the system of non-linear algebraic equations giving the updated values at the same point, u * i . A further approximation made in [39] is to assume there are no cross terms, e.g.
where, for example, for two variables, ϕ, c we have In the present code we expand the Jacobian to cover, in 2D, 4 cells, and in 3D, 8 cells. For two field variables, ϕ and c, this results in the 16 × 16 Jacobian matrix, J in 3D. Now because of the complexity of an arbitrary anisotropy function and complicated data base functions, we find it too cumbersome to construct a new Jacobian with every change of model. We therefore adopted an approximation to the Jacobian using finite differences. The reader should note that this construction is not related to the mesh size. The approximate Jacobian may be denoted Here, F is a vector of length, 16, and δu is a variation of each of the two variables in all 8 cells. So, the resulting Jacobian is 16 × 16 as stated.