The hodograph equation for slow and fast anisotropic interface propagation

Using the model of fast phase transitions and previously reported equation of the Gibbs–Thomson-type, we develop an equation for the anisotropic interface motion of the Herring–Gibbs–Thomson-type. The derived equation takes the form of a hodograph equation and in its particular case describes motion by mean interface curvature, the relationship ‘velocity—Gibbs free energy’, Klein–Gordon and Born–Infeld equations related to the anisotropic propagation of various interfaces. Comparison of the present model predictions with the molecular-dynamics simulation data on nickel crystal growth (obtained by Jeffrey J. Hoyt et al. and published in Acta Mater. 47 (1999) 3181) confirms the validity of the derived hodograph equation as applicable to the slow and fast modes of interface propagation. This article is part of the theme issue ‘Transport phenomena in complex systems (part 1)’.


Introduction
Anisotropy of interfaces plays a crucial role in the formation of equilibrium shapes [1], changing of growth direction of crystals to the preferable one at a critical  governing parameter [2], selection of stable mode of dendritic growth [3] and facetting with the formation of microscopic defects in the bulk of crystals [4]. Figure 1 shows a region around a dendritic tip, the growth of which has been selected and faceted due to the existence of anisotropy of solid-liquid interface. From the near equilibrium properties, these anisotropic interfaces lead to crystallographically oriented properties such as free energy, magnetization and surface tension [5]. From the kinetic properties, the interface mobility and coefficients of atomic attachment to the interface strongly depend on the crystallographic faces orientation [6,7].
In the last century there has been great success in the description of equilibrium shapes and their motion. Herring [8] generalized the formalism of Wulff on thermodynamics of crystal surfaces. Cahn & Hoffman [9,10] introduced a useful construction, the so-called ξ -vector, to characterize equilibrium shapes and missing orientations of sharp interfaces. Caginalp     Kobayashi [11][12][13] were the first to introduce the anisotropy in the diffuse interface formalism with the surface tension anisotropy by allowing the gradient energy factor to be dependent on the orientation of the phase interface. Wheeler & McFadden [14][15][16] provided description of anisotropic interfaces using the Cahn-Hoffman ξ -vector within the framework of the phase field model. Avoiding unstable interfacial orientations, the growth of highly anisotropic interfaces with facetting and edges was analysed by Eggleston et al. [17] and Debierre et al. [18]. Finally, an overview by Sekerka [19] presented approaches in the description of anisotropic crystal growth.
In the present work, we extend the description of anisotropic interfaces to the highly rapid regimes of their motion with the appearance of local non-equilibrium effects [20]. With this aim, a kinetic phase field model [21,22], which would be reduced to the single equation of motion called 'the hodograph equation of interface', is used [23,24]. The hodograph equation predicts non-stationary regimes as well as steady-state regimes of interface motion and it has been applied, for instance, to quantitative estimations of non-stationarity periods of dendrite growth [25]. Using the formalism of the Cahn-Hoffman ξ -vector and the presently developed anisotropic phase field royalsocietypublishing.org/journal/rsta Phil. Trans. R. Soc. A 379:  [14] who obtained the sharp interface limit corresponding to the case for which the diffuse interface width is small compared to a characteristic macroscopic length scale.
The article is organized as follows. As an advancement of the isotropic phase field model ( §2) we introduce the kinetic phase field model with crystalline anisotropy in its hyperbolic formulation ( §3). The obtained hodograph equation in the form of Herring-Gibbs-Thomson-type equation is compared with the previous equations of equilibrium, accelerationvelocity-dependent equation as well as Born-Infeld and Klein-Gordon equations describing the propagation of isotropic and anisotropic interfaces ( §4). A quantitative comparison of linear and nonlinear equations for interface motion which follow from the obtained kinetic equation for anisotropic crystalline interface is given ( §5). A summary of our conclusions is presented ( §6). Finally, electronic supplementary material, [26] and appendix A add the material for the derivation of the present hodograph equation applicable to the slow and fast modes of interface propagation.

Isotropic phase field equation
The isotropic phase field equation and its solution were obtained in the work [23] and analysed in comparison with MD-data in [27]. In this section, we recollect the results of these works necessary to reproduce anisotropic advancement of the hyperbolic phase field model.
The free energy functional for the entire system of the volume v 0 in the isotropic case can be written in the form which yields the dimensional isotropic form of the dynamic equation as Here, φ is the phase field variable defining the phase state as φ = 0, in the liquid phase, 1, in the solid phase, τ φ is the relaxation time of the gradient flow for the phase field ∂φ/∂t, M φ is the mobility, C is the solute concentration and T is the temperature. The Gibbs free energy change on transformation G takes into account direction of transformation and the variety of transformations obtained for enthalpy of fusion, H f (T) [28], dilute mixtures [29], and functions obtained from thermodynamics databases [30]. Note that G i (T, C)(i = l, s) is the Gibbs free energy of the phase, and the indexes l and s are related to the liquid and solid phases, respectively. The interpolation function p(φ) and the double-well function g(φ) are defined by [31] According to definition (2.3), the interpolation function p(φ) varies monotonically from p(0) = 0 to with the diffuse interface stationary width and under the boundary conditions φ = 1 as x ≤ δ and φ = 0 as x ≥ δ. Then, the surface energy γ of the isotropic interface is given by which is proportional to the product of the interface width δ and the energy per unit volume associated with the barrier height W φ . In the dynamics, with the boundary conditions φ → 1 as x − Vt → −∞ and φ → 0 as x − Vt → +∞, with the constant velocity V limited by V φ as a maximum speed of phase field propagation, the velocity-corrected effective interface thickness, 11) and the mobility M φ related to the isotropic interface mobility μ as First, the particular solution (2.9)-(2.12) with the hyperbolic tangent function follows from the general set of analytical solutions of Allen-Cahn-type equations [33] which is given by equation (2.2). Second, the interface velocity, V, cannot exceed the maximum speed of disturbance propagation in the phase field, because the phase field itself dictates the interface shape and its velocity, i.e. V < V φ in the solutions to equations (2.9)-(2.12). Third, with regard to the effective interface thickness (2.11), one has to note two important issues: (i) with increasing interface velocity, should become smaller than the constant interface width δ that has been chosen as a reference for the interface thickness in equilibrium state, equation (2.7); (ii) within the limit V → V φ , one gets → 0, therefore, the phase field variation will be steeper with the tendency to build up a sharp interface as the velocity increases.
In the dimensionless isotropic form, equation (2.2) can be re-written taking into account equations (2.8) and (2.12) as To obtain the Herring-Gibbs-Thomson-type equation following from the hyperbolic phase field equation (2.13), we shall obtain the sharp interface limit corresponding to the case that the diffuse interface width is small compared to a characteristic macroscopic length scale , which can be taken to represent a typical radius of interfacial curvature in the non-planar case [14]. More specifically, we consider the distinguished limit while maintaining a finite value for the ratio (2.15) in order that the surface energy remains finite in the sharp interface limit. Using these scales, it is convenient to choose the length and coordinates in units of , the energy density G and energy barrier height W φ in unit of latent heat H f , the time t and relaxation time τ φ in units of a diffusive time scale τ consistent with the accompanying energy diffusion processes, the surface energy γ in units of H f and the interface kinetics given by μ in units of /( H f τ ). Then, the isotropic form of equation (2.2) follows from equation (2.13) as where plays the role of a small parameter in the subsequent asymptotic treatment and it is defined by This ratio means that the high value of barrier W φ (T, C) between phases provides a smallness of and the asymptotic limit → 0.
3. Hyperbolic phase field model with anisotropy (a) Free energy and scaled phase field equation Consider a binary system consisting of solvent and solute under isothermal condition, with the temperature T being constant in the overall system. The system is undergoing phase transition, solidification/melting, from the undercooled/overheated state. Taking the existence of the anisotropic diffuse solid-liquid interface into account [34,35], the free energy functional in units of H f for entire system of the volume v 0 is described as where length, time, Gibbs potential, G(T, C, φ, ∂φ/∂t) and homogeneous functions of surface energy, γ (n), are measured in units of , τ , H f and H f , respectively, and γ is the surface energy of the interface depending on the normal vector n pointing from solid to liquid, The homogeneous extension of γ (n) is given by [36] |∇φ|γ − ∇φ |∇φ| = γ (−∇φ), (3.3) such that the free energy functional (3.1) becomes where the Gibbs potential G(T, C, φ, ∂φ/∂t) is described by the expression [37,38] which has the local equilibrium contribution , (3.6) and the local non-equilibrium contribution In addition to the interpolation function p(φ) and the double-well function g(φ) given by equation (2.5), the contributions (3.6) and (3.7) include: the phenomenological coefficients α φ (T, C) proportional to the relaxation time τ φ of the gradient flow ∂φ/∂t, the barrier W φ (T, C) between phases as well as the Gibbs free energies G l (T, C) and G s (T, C).
A stable evolution of the entire system is given by the Lyapunov condition of non-positive change of the total Gibbs free energy. For the functional (3.1), this condition gives the inequality from which one finds the following phase field equation (appendix A) Note that M φ and ξ are measured in units of ( H f τ ) −1 and H f , respectively. Here, the ξ -vector of Hoffman & Cahn [9] is described by where the second equality in equation (3.10) defines the ξ -vector as a homogeneous function.

(b) The hodograph equation as a form of the Herring-Gibbs-Thomson-type equation
The anisotropic form of the phase field equation can be written if we assume that the surface energy γ and the interface mobility μ are homogeneous functions of the first degree in order to include their anisotropy, i.e. γ (−∇φ) = −γ (∇φ) and μ(−∇φ) = −μ(∇φ), as for ξ -vector equation (3.10). Indeed, taking into account these conditions and using the same treatments from §2 for equation (2.2), one can find that the anisotropic equation (3.9) takes the following form: where equation (2.17) holds for the definition of small parameter . If the inertial term is neglected, τ φ → 0 in equation (3.11), the difference in Gibbs free energy is taken as [28] G where μ(n) is the orientation dependent interface mobility, V n (t) and A n (t) = ∂V n (t)/∂t stand, respectively, for the velocity and acceleration normal to the interface, 14) and the divergency of the ξ -vector is is the term dependent on the curvatures κ i , surface energy γ (n), orientation-dependent interface mobility μ(n) and their first order derivatives with respect to θ i . In the particular case of a convex shape having negative principal curvatures κ i (κ i < 0 for i = 1, 2), the term D n (t) plays a role of deceleration. Equation   This equation also follows from the zero variational derivative from free energy (3.1) with ∂φ/∂t = 0 [4]. Such a condition defines the minimal thermodynamic potential with respect to the transfer of particles from one phase to the other. Therefore, equation (4.1) describes the shape of anisotropic particle in equilibrium, presenting a balance between a volumetric tendency to exchange and a surface tendency to save a shape of coexisting phases.

(b) Acceleration-velocity-dependent Gibbs-Thomson equation
In the case of the isotropic phase field model, the surface energy, the interface mobility and the maximum speed of disturbances are independent of the orientation, i.e. γ (n) = γ , μ(n) = μ and V φ (n) = V φ have constant averaged values. Equation (3.16) is then simplifying for i = 1, 2 due to the zero derivatives and Taking into account equations (2.7), (2.8) and (2.12) established in the case of the isotropic model, the following relations among the model parameters are obtained where we also used equation (4.3). Using solution (2.6), and taking into account equations (2.7), (2.8), (2.12) and (4.4), the time-dependent thickness of the interface (t) is now given by [23] (t) = 2 3 δ 1 − W 2 n (t), W n < 1. (4.5) The inequality in equation (4.5) has the same meaning as for solutions (2.9)-(2.12): because the phase field itself dictates the interface shape and its velocity, the interface velocity cannot exceed the maximum speed of disturbance propagation in the phase field, W n < 1. By the above considerations in the isotropic scheme, we recover from equation (3.16) our generalized acceleration-velocity isotropic Gibbs-Thomson equation found in [23] and analysed for (100)-Ni-crystal orientation in [27] with the same parameters, with the negative curvature K = κ 1 + κ 2 , the phase field assumes the value φ = 1 and φ = 0 in the solid and liquid phases, respectively. At small and moderate velocity, W n 1 and with the absence of the driving force, G = 0, equation (4.6) arrives at In two dimensions, this equation was used by Gurtin & Podio-Guidugli [39] to explain the evolution of interface possessing effective inertia consistently with experimental data on the oscillation of quantum crystals [40].

(c) Velocity-dependent Herring equation
In the case of absence of inertial effects, τ φ → 0 (V φ → ∞, A n = 0) and constant normal velocity V n = const, equation ( which describes the interface motion of anisotropic particle due to imposed Gibbs free energy change on transformation, G, and both principal curvatures, κ 1 and κ 2 . A couple of important particular cases can be outlined, which follow from equation (4.9). First, if the driving force is absent G = 0 (i.e. with the absence of supersaturation or supercooling), equation (4.9) predicts the interface motion by both principal curvatures (4.10) The latter effect can be seen as mean curvature flow described by the Allen-Cahn equation [36,42,43] if both principal curvatures are equal to the mean curvature of the surface, i.e. κ = κ i (i = 1, 2) and θ = θ i (i = 1, 2). Therefore, equation (4.10) reduces to Second, in the case of planar interface, κ 1 = κ 2 = 0, a simplest equation of motion, can be found from equation (4.9). This equation follows from the classic theory of irreversible processes and from traditional phase field model [11,34].

(d) Born-Infeld equation
One of the specific cases of nonlinear wave propagation can be found from equation (3.16) under the absence of driving force, G = 0, and with the infinite orientation-dependent interface mobility, μ(n) → ∞

17), is described by
with K = κ 1 + κ 2 standing for the negative curvature. Equation (4.14) can be directly obtained from equation (4.6) and it is used in nonlinear electrodynamics [44,45]. Such an analogy of nonlinear waves and the fast phase interface propagation can be clarified in forthcoming works.

(e) Klein-Gordon equation
Accepting τ φ → ∞ and M φ /τ φ ≡ const in equation (3.9), one gets the undamped Klein-Gordon equation the Universe [46]. In this sense, equation (4.15) can be considered as the Klein-Gordon equation for the anisotropic field propagation of the inflation stage of the matter. Using the definition of the ξ -vector of Hoffman & Cahn (3.10), equation (4.15) can be transformed to the isotropic version of the undamped Klein-Gordon equation [23]. A simplest form of the damped Klein-Gordon equation can also be found from equation (3.9) with τ φ ≡ const but within the limit γ (−∇φ) → 0: This equation describes the oscillatory motion of matter within the action of the potential function g(φ). Such a type of front propagation is used for qualitative analysis of the inflation stages of various fields [46].

Comparison with data of atomistic simulations
In this section, we shall transform from the dimensionless variables and functions (as described in §2, 3 and 4) to the dimension equations and expressions as is given by table 1.
At small driving force of solidification or melting, interface velocity has a linear dependence on undercooling or overheating, respectively [6] that is also predicted by equation (4.6) in comparison with MD-data [23,24,27]. At large driving force, MD-simulation predicts interface velocity with non-linearity of two types: velocity-undercooling relationship with saturation [47] and with maximum by the crystal growth velocity [48,49]. These nonlinear dependencies in growth kinetics of crystals are well described by the solution of equation (4.6) and by the hyperbolic phase field equation (2.2) [22]. The kinetic undercooling obtained from the steadystate mode predicted by equation (4.6) included in a general undecooling balance at the tip of dendrite allows us to describe the solidification kinetics in glass-forming alloys [50].
where the phase field diffusion coefficient, D φ (n), and the maximum phase field speed, V φ (n), are considered as two parameters which define the relaxation time of the gradient flow of the phase field by With the local equilibrium limit, i.e. with τ φ → 0 (V φ (n) → ∞), equation (5.1) reduces to the linear relation V n = μ(n) G, which is obtained from the parabolic phase field equation and which can exhibit nonlinear behaviour only due to nonlinearity of the function D φ (n) (see discussion in [22]).

(b) Anisotropic functions
By employing cubic harmonics [51,52] and adopting established conventions [53][54][55], the interfacial free energy γ (n) can be written as   where n i (i = x, y, z) are the Cartesian components of the unit normal vector n.
Hoyt et al. [56] defined the values of stiffness for different crystallographic orientations in the solid-liquid interface energy needed for parametrization of γ (n). Following the results of their work, the anisotropic function f (n) and its corresponding stiffness f (n) + d 2 f (n)/dθ 2 can be evaluated around the vanishing angle between the normal n to the interface and the normal to the crystal face, i.e. θ = 0 for different orientations. Tables 2 and 3  (c) Application to growth of nickel crystal (i) Free energy for different crystal orientation These values can be used to compute the interfacial free energies with respect to 100 , 110 and 111 crystal orientation for Ni, i.e. γ 100 , γ 110 and γ 111 , respectively, from the equations shown in table 2. The comparison with the results of works [57,58], and the values of the interfacial free energies obtained by using the average interfacial energy, γ 0 from [59,60] is shown in table 4.

(ii) Velocity-undercooling relationship
In the case of a pure (chemically one component, elemental) system, the driving force, G, can be given by equation (3.12). Therefore, equation ( from which one can identify the expression of coefficients of crystal growth kinetics, μ klm , as where k, l and m are the Miller indices.  in table 2 and table 3, yield a good agreement in comparison with [57,58] and those obtained using the averaged interface energy of [59,60], see table 4. -Using material parameters of Ni from table 5, figure 2 shows very good predictions of MD-data due to Hoyt et al. [47] by the hodograph equation, (3.16), for planar interface, κ 1 = κ 1 = 0, and constant normal velocity V n (t) = const given by equation (5.1), where the driving force, G, is given by equation (3.12). - Table 5 shows the agreement for the obtained values of interface kinetic coefficient, equation (5.8), in comparison with those found in the literature [6,47]. -Using expansion of the interface energy (5.3) and (5.4), we obtained values of kinetic coefficients in three crystallographic directions by the fitting of equations (5.7) and (5.8) to data of MD-simulation. In this respect, one should note that the values of these kinetic coefficients μ klm correlate well with the values obtained by the authors of the work [47] for the same MD-data, see table 5. This confirms the adequacy of the application of the kinetic equation (5.7) together with (5.8) to the description of anisotropic interfaces in a wide range of undercooling.
Note that for each 100 , 110 and 111 orientations with the values of interfacial free energies γ 100 , γ 110 and γ 111 , the melting temperature, T m and the enthalpy of fusion, H f , respectively (table 6), the phase field diffusion D φ,100 , D φ,110 and D φ,111 , and the maximum speeds of phase field propagation, V φ,100 , V φ,110 and V φ,111 , which are considered as parameters, are obtained by fitting MD-data due to Hoyt et al. [47]. These latter are shown in the 10th and 12th rows of

Conclusion
The present work has been devoted to the derivation and analysis of the slow and fast interface propagation described by the hodograph equation following from the phase field model. The obtained solution of the hodograph equation describes the amplitude, width and velocity of the anisotropic interface in the slow and fast regimes of dynamics. We show that the obtained hodograph equation is consistent with classical equations of the Herring-Gibbs-Thomsontype as well as the velocity-dependent Herring equation for crystal growth, Klein-Gordon equation for the front's propagation of the inflation stage of matter, and Born-Infeld equation of the nonlinear electrodynamics. The first benchmarks of the derived hodograph equation of anisotropic interface motion show the consistency of its solution with the steady-state growth of nickel crystals obtained for different crystallographic directions using molecular dynamics simulation of Hoyt et al. [47].
Data accessibility. Electronic supplementary material on asymptotic analysis of the hyperbolic phase field model are attached to the main text of the manuscript.