Nano-faceted crystals answer the call for self-assembled, physico-chemically tailored materials, with those arising from a kinetically mediated response to free-energy disequilibria (thermokinetics) holding the greatest promise. The dynamics of slightly undercooled crystal–melt interfaces possessing strongly anisotropic and curvature-dependent surface energy and evolving under attachment–detachment limited kinetics offer a model system for the study of thermokinetic effects. The fundamental non-equilibrium feature of this dynamics is explicated through our discovery of one-dimensional convex and concave translating fronts (solitons) whose constant asymptotic angles provably deviate from the thermodynamically expected Wulff angles in direct proportion to the degree of undercooling. These thermokinetic solitons induce a novel emergent facet dynamics, which is exactly characterized via an original geometric matched-asymptotic analysis. We thereby discover an emergent parabolic symmetry of its coarsening facet ensembles, which naturally implies the universal scaling law for the growth in time t of the characteristic length .
Nano-faceting of material interfaces [1,2] is a paradigmatic, non-equilibrium self-assembly process which arises in a wide variety of physical settings; for example, high-efficiency photoelectrochemical cells yielding solar-energy storage through hydrogen production , and enantiomer-specific heterogeneous catalysts with application to biology . Exotic nanocrystal motifs and unexpected orderings abound here, such as tetrahexahedra from electrochemical  synthesis, truncated Wulff shapes in solution-phase  synthesis, adorbate-induced nano-faceting of single-crystal surfaces , and the self-ordering of epitaxially grown quantum dot arrays , to name but a few. Since the pioneering low-energy electron diffraction work of MacRae , advances of in situ and ex situ instrumentation for nanoscale imaging [10–12] are permitting experimentalists to elucidate the essential physics and chemistry underlying this tantalizing range of nano-faceting phenomena, and the pace of discovery shows no signs of abating. However, theories which provide quantitative prediction of such shapes and shape distributions (morphometrics) are still in their scientific infancy. Since physico-chemical characteristics of nanocrystals are known to exhibit an exquisite sensitivity to their shape , morphometric theories will be crucially important for the scientific development of engineered materials with prescribed properties for applications ranging from bio-sensing  to green nanotechnologies .
Certain planar interfaces of macroscopic single crystals can be thermodynamically induced to spontaneously decompose into fully faceted morphologies by thermal [16,17], chemical [7,18–21] or electrochemical means [13,22]. Experimental observations of spinodally decomposing interfaces [23,24], such as those for thermally quenched Si(111)  or the oxygen-induced faceting of Cu(115) , reveal an ensuing coarsening of the resulting fully faceted interfaces, and the concomitant emergence of scaling regimes, wherein a power law governing the time t evolution of the characteristic length of the interface's morphology appears [16,17], The range of scaling exponents n that have been empirically discerned [25–27] suggest they depend on the dominant mass-transport mechanisms, and possibly even the symmetry group of the crystal's Wulff facets. Furthermore, normalized length distributions of such faceted morphologies within these scaling regimes (e.g. edge lengths) have also been empirically observed to be given by universal distributions .
Much of the non-equilibrium nano-faceting phenomena noted above lie definitively outside the predictive scope of purely energetic or purely kinetic considerations , but rather reflects the interface's kinetically mediated response to free-energy disequilibria [7,13]. The ensuing thermokinetic interfacial dynamics, which encodes a thermo-mechanical balance between the local rates of surface-free-energy/bulk-chemical-potential release and a kinetics-induced entropy production , finds expression in multiple length- and time-scale geometric partial differential equations (PDEs).
The isothermal growth of a nano-faceted crystal into its undercooled melt has been thermodynamically modelled at a continuum level by Herring, who adapted Wulff's assumption that the interfacial energy is (geometrically) non-convex with respect to the unit-normal n to the interface  by also including a dependence on the interfacial curvature . Upon assuming attachment–detachment kinetics to be the dominant mass-transport mechanism, one concludes, upon suitable spatio-temporal scaling, that such a one-dimensional Herring interface , oriented with interfacial-normal n pointing into the melt (figure 1), evolves according to the hyperbolic–parabolic, geometric PDE  1.1 where , ϑ and κ, respectively, denote the normal velocity, local-inclination angle and curvature of , ∂s is the partial derivative with respect to the interfacial arc-length s, while the surface stiffness reflects the anisotropy of the interfacial energy, and ε is a dimensionless constant that encodes the degree of undercooling. We note in passing that the Burton–Cabrera–Frank step-flow models for thin-film deposition on crystal substrates  may possess step-bunching instabilities that lead to mound formation, and the coarse-grained description thereof yields continuum models which qualitatively mimic equation (1.1): see Sec. 4.9.8 of Michely & Krug .
In this article, we analytically characterize the slow facet dynamics that effectively govern the solutions to equation (1.1) in the slightly undercooled regime 0<ε≪1. To determine this effective dynamics, which mathematically emerges from equation (1.1) by letting ϵ↘0, we develop novel geometric matched-asymptotic methods; these significantly extend the analysis  previously applied to the convective Cahn–Hilliard equation [35,36]. The emergent facet dynamics we uncover is highly non-trivial since it does not coincide with simply setting ε=0 in equation (1.1), which only yields a mathematically ill-posed (backward parabolic) equation when .
After deriving equation (1.1) in §2, we proceed in §3 to establish the existence of convex (+) and concave (−) translating fronts that asymptote to angles and (thermokinetic angles), which both depend on the degree of undercooling ε and also deviate from the thermodynamically expected Wulff angles ±ω; a thermokinetic effect. In §4, we develop a matched-asymptotic analysis of equation (1.1) which analytically demonstrates that this thermokinetic effect, , slowly drives the large-scale features of solutions between coarsening events. We thereby discover that the interfacial evolution is captured to leading order by purely Wulff-faceted interfaces slowly evolving under an explicitly derived intrinsic dynamics. Furthermore, this emergent facet dynamics possesses a Peclet length, below which a parabolic-scaling symmetry operates; this is the geometric extension of the emergent Peclet length for the convective Cahn–Hilliard equation . In §5, we turn to an analysis of the coarsening dynamics of solutions to equation (1.1) that originate from a small random perturbation of a thermodynamically unstable, macroscopic facet orientation ϑ=0. We theoretically predict and numerically validate that the time t evolution of the characteristic (effective) facet length is governed by the universal scaling law 1.2 We also empirically verify that the associated one-point statistics of the normalized facet lengths are universally governed by a unique time-invariant distribution.
2. Herring's oriented interface dynamics
The surface free energy of a crystal is classically presumed to depend solely on its local crystallographic orientations , and thereby permits edges to arise between minimum-energy facets of equilibrium crystals. Reflecting the configurational abhorrence of an atom to sitting at such an edge, Herring augmented the Wulff energy with a curvature-dependent term  to enforce, over a prescribed nanoscopic length scale, α, a smooth transition between the surface normals of adjacent equilibrium facets; we henceforth refer to α as the bending length. For a one-dimensional interface (schematically illustrated in figure 1) Herring's model energy takes the form 2.1 where γ0 sets the interfacial energy scale, the dimensionless function encodes the dependence on the interfacial orientation, θ, relative to a prescribed crystallographic axis, and s, ϑ and κ:=∂sϑ denote the arc-length, local-inclination angle relative to a prescribed polar axis, and signed curvature of the interface , respectively. In what follows, we will restrict ourselves, for simplicity, to fourfold symmetric surface energies , . Recall that the sign of the surface stiffness, (′ denotes the derivative), determines the stability of a facet orientation. We assume the orientation θ=0 both is thermodynamically unstable, , and sits between two neighbouring thermodynamically stable Wulff angles, ±ω; i.e. 2.2
Now consider such a one-dimensional single-crystal interface in isothermal contact with its undercooled melt (figure 1). The jump in the specific chemical potential across the interface [[μ]] that is induced by the melt's temperature T being below the melting temperature of the crystal Tm, namely T<Tm, is given by where L is the specific latent heat of fusion. We presume that attachment–detachment kinetics is the sole mass-transport mechanism, and let denote the associated interface mobility—an inverse-velocity scale associated with the interfacial entropy production of solidification . We take the ideal situation in which the melt- and crystal-mass densities, ρm and ρc, coincide; namely, ρm=ρc=ρ. In this case, no convective flow field is induced during solidification/melting of the crystal. Choosing the orientation of the interface to be given by the outward unit normal pointing into the melt, the dynamics of the interfacial evolution is encoded in the thermo-mechanically derived geometric PDE where denotes the normal interfacial velocity, the surface stiffness (′ denotes the derivative), and ∂s denotes the partial derivative with respect to the interfacial arc-length s.
The ratio of the surface energy scale γ0 to the jump in chemical potential per unit volume, ρ[[μ]], sets the characteristic length scale on which the surface energy and the bulk free energy balance, namely 2.3 Its ratio with the purely kinetic velocity scale then yields the intrinsic time scale 2.4 Scaling to the associated dimensionless lengths and times, i.e. and , we arrive at the (dimensionless) form of Herring's oriented interface dynamics 2.5 where the dimensionless parameter ε≠0 measures the ratio of the bending length α to the length scale appearing in equation (2.3), namely 2.6 For copper L=2×105 J kg−1, ρ=8.5×103 kg m−3, γ0=1.75 J m−2 and estimating α≈1 nm, one obtains, for a dimensionless undercooling of 1%, the value ε≈0.01.
3. Thermokinetic angles of translating fronts
We recall that a translating front is a time-evolving curve of the form 3.1 where is a fixed curve, V is a constant vector and denotes the curve at time t obtained by translating the fixed curve by an amount tV; we naturally refer to as the shape of the translating front. Letting denote an arc-length parametrization of , where i and j denote the standard unit basis vectors of a prescribed stationary two-dimensional Cartesian reference frame, we note for future reference that naturally furnishes an arc-length s parametrization of .
We now turn to the analysis of translating fronts which solve Herring's oriented interface dynamics (2.5) when 0<ε≪1: what we henceforth call the slightly driven regime. We focus in particular on identifying those with symmetric unbounded profiles which are either strictly convex or strictly concave, and also asymptotically flat at infinity: 3.2 Such profiles and their asymptotic angle Ω will shortly be shown to depend on ϵ and the sign of their curvature (+ or −). When we need to, we will later use the notation and to denote the asymptotic angle of the corresponding convex profile and the concave profile .
Note that reflection in the Cartesian y-axis, (x,y)↦(−x,y), is a symmetry of equation (2.5) since is invariant under the inversion of the polar angle, namely θ↦−θ. It is therefore natural to seek travelling-front solutions to equation (2.5) which also respect this symmetry, and so we take to be symmetric about the y-axis, and set the velocity V to be parallel to the y-axis, : 3.3 Without loss of generality, we choose that arc-length s of the oriented curve for which at s=0 the corresponding local angle of is aligned with the polar axis θ=0. Denoting the corresponding arc-length parametrizations of the local angle ϑ of by , we have 3.4
Anticipating the appearance of the inner-length scale S:=s/ε, over which a sharp transition in the interfacial angle of occurs, we introduce the scaled vector function 3.5 which parametrizes the corresponding scaled profile . By differentiating equation (3.5) with respect to S, one concludes , which confirms that S serves as an arc-length coordinate on . Just as before, we let denote the corresponding arc-length parametrization of the local angle Θ of 3.6 Recasting equation (3.5) as and differentiating both sides with respect to s we find and therefore also find that . In other words, also provides a parametrization of the local angle ϑ of with respect to the inner-scale S:=s/ε: 3.7 Note therefore that on we have the change of variable formula 3.8
To determine whether the translating front (3.3) with is a solution of (2.5), we need to calculate its curvature κ, normal velocity and surface Laplacian of the curvature : 3.9 We conclude that such a translating front is a solution of equation (2.5) precisely when the arc-length parametrization of the local angle Θ of its scaled profile solves 3.10
By incorporating the asymptotic flatness , namely equation (3.2), into equation (3.7), one finds 3.11 Taking the limit in equation (3.10) and using (3.11), one immediately concludes that exists and its value is given by But the only limiting value consistent with equation (3.11) is 0, so we have . It therefore follows that the front speed is necessarily given by 3.12 Substituting equation (3.12) into equation (3.10) we arrive at the third-order differential equation for 3.13 Recalling that the signed curvature of is given by , we may re-write equation (3.13) as a geometric differential equation for the scaled profile , namely 3.14
In the case where is strictly convex (concave), i.e. ∂sΘ>0 (∂SΘ<0), then its local angle Θ is monotone increasing (decreasing) with respect to its arc-length S. Therefore, the local angle of , namely Θ, may serve as natural parametrization of by associating the local angle Θ at each . Therefore, on the arc-length derivative ∂S is related to the angular derivative ∂Θ through the chain rule by 3.15 Furthermore, on the following arc-length and angular limits are equivalent: 3.16 By using equations (3.15) and (3.16), we may recast the assumed asymptotic flatness of expressed in arc-length terms by equation (3.2) into the angular equivalents 3.17 Furthermore, noting the identity on 3.18 we may re-express the geometric differential equation (3.14) in the angular form 3.19 where .
In summary, given the profile of a convex (concave) travelling front solution of equation (2.5) which is asymptotically flat at infinity, we have found that the angular parametrization of the curvature of the scaled profile , is a fixed-sign solution of the two-point boundary-value problem given by 3.20 Note that the second-order partial differential operator will generically need only two boundary conditions to be well posed, and so the prescription of four boundary conditions here can be expected to impose a restriction (compatibility condition) on Ω. Last, we may naturally view the solutions to equation (3.20) with 0<ε≪1 as regular perturbations in ε of the symmetric curvature profiles defined by 3.21 subject to the same boundary conditions and sign condition.
There is a natural converse to the above whereby, given a function and angle Ω which solve equation (3.20), one may readily construct associated travelling front solutions to equation (2.5). First, let denote the profile parametrized by the function which arises as the unique solution to the first-order initial-value problem 3.22 is, by design, the angular parametrization of , which will also be unbounded if its associated arc-length 3.23 diverges as Θ→Ω. It then follows that the travelling front defined by translating the scaled profile with velocity , i.e. 3.24 is a solution of equation (2.5).
(a) Exact translating fronts and thermokinetic angles and
Here we present exact convex and concave travelling front solutions to equation (2.5) that arise for the particular surface energy 3.25 with constants a>0 and ω∈(−π/2,π/2) given. One readily verifies that this choice of surface energy meets the conditions expressed in equation (2.2) with Wulff angles ±ω. Our approach to finding these exact solutions is to use the converse construction noted around equation (3.22). To that end, we seek a positive (negative) solution to equation (3.20), and Ω+>0 ( and Ω−>0), within the ansatz where Ω+>0 and Ω−>0 denote generic positive and negative angles.
Rather remarkably, one finds that both with Ω=Ω+ and with Ω=Ω− each solve the second-order PDE appearing in equation (3.20) precisely when the corresponding angles Ω± are solutions of the trigonometric-polynomial equation 3.26
The range of ε for which there exist solutions to equation (3.26) is readily determined by considering the graph (figure 2). One readily finds that for each there exist positive (negative) solutions Ω+>0 (Ω−<0) of equation (3.26). Furthermore, it is only the least positive (negative) solution Ω+(Ω+) of equation (3.26) for which with with Ω=Ω−) meet the non-zero condition equation (3.20); namely for all Θ∈(−Ω±,Ω±). For each we therefore let Ω+ε(Ω+ε) denote the least positive (negative) solution of equation (3.26) and introduce the corresponding curvature functions 3.27
Now recall the travelling front solutions to equation (2.5) constructed from solutions to equations (3.20) via equations (3.22) and (3.24). In the case of equations (3.27) with ε chosen in the appropriate range, we are able to solve equations (3.22) exactly, and thereby explicitly characterize the corresponding convex and concave profiles, here denoted by and , respectively. Introducing the angular parametrization and , namely 3.28 we find that the cartesian components and are explicitly given by 3.29 The convex and concave curves and are oriented in the direction of increasing and decreasing Θ, respectively. Now choosing the arc-length S for which Θ=0 at S=0, we may identify the arc-length parametrization of the local angle of , denoted by , through 3.30 This first-order initial-value problem is exactly solvable, and we find 3.31
In conclusion, the scaled oriented profiles , translating with the respective velocity , namely 3.32 are each an exact solution of Herring's equation (2.5). Recall that the positively oriented arc-length parametrizations of the local angle of , namely , follow from that of according to , and hence 3.33 Since the Wulff angles ±ω are determined through a free-energy minimization, the deviation from ±ω noted in (3.33) marks the profiles as truly non-equilibrium profiles, whose asymptotic orientations are thermokinetically determined. Furthermore, an examination of the ε-expansion of these thermokinetic angles , namely 3.34 with , reveals the asymmetric nature of this offset. It reveals that the asymptotes of the convex profiles are steeper than would be expected from thermodynamic considerations, while those of the concave profiles are shallower (figure 3). This thermokinetic offset from Wulff angles, as exhibited in the ε-expansion (3.34), will now be shown to be universally true for the convex and concave travelling front solutions of the Herring equation (2.5) at small driving 0<ε≪1, provided only that the surface energy meets the natural faceting conditions expressed in (2.2).
(b) General asymptotic analysis of and
Here we present an analysis of convex and concave travelling front solutions to equation (2.5) for a generic non-convex surface energy, subject only to the natural thermodynamic conditions appearing in equation (2.2). Given the ε-family of exact solutions of this type which we have explicitly constructed for the particular surface energy equation (3.25), it is natural to assume the existence of such an ε-parametrized family of solutions for the general case equation (2.2), at least for 0<ε≪1.
We assume there exists for each 0<ε≪1 a corresponding unbounded, symmetric, convex (concave) profile, (), that is asymptotically flat in accordance with equation (3.2), with , and such that 3.35 solves Herring's equation (2.5). Introducing the scaled profiles , we find as before that the angular parametrization of the curvatures of the respective profiles , namely 3.36 solve equation (3.20) with , and in particular 3.37 Recalling equation (3.21), we may naturally view the ε-family of solutions to equation (3.37), namely with 0<ε≪1, as regular asymptotic expansion in ε of the convex and concave curvature profiles and that follow from setting ε=0 in equation (3.37), namely 3.38
Viewing equation (3.37) as being of the form , where w and f are symmetric in Θ, and also noting that the only symmetric solutions of the associated homogeneous equation ( f=0) are scalar multiples of , we conclude that 3.39 where is an as yet unknown constant, and is the particular symmetric function 3.40 which is obtained by variation of constants.
Recalling the asymptotic flatness of the curvature expressed in equation (3.21), namely and , we find, on taking the corresponding limits in equation (3.39) and then solving for , that 3.41 Substituting this expression for into equation (3.39), we arrive at 3.42 Now differentiating equation (3.42) w.r.t. Θ to yield 3.43 and also recalling the specific asymptotic conditions of equation (3.21) 3.44 we find, upon taking the corresponding limits in equation (3.43), that 3.45 Directly calculating from equation (3.40), we re-express equation (3.45) as 3.46 where, as before, denotes the angular derivative of .
Setting ε=0 first in (3.46), we find 3.47 from which it follows from equation (2.2) that 3.48 Noting this and setting ε=0 in equation (3.42), we immediately conclude , where 3.49 Viewed as an angular parametrization of a curve, the curvature profile coincides with that of the regularized equilibrium profiles previously identified in . Now recalling equation (3.48), we conclude that the ε-expansion of around ε=0 assumes the form 3.50 In addition, given that , and , we are able to deduce from equation (3.46) and L’Hôpital's rule that 3.51 We therefore deduce that η±=η>0, where the universal Herring offset η is given by 3.52 which when inserted into equation (3.50) provides the complete generalization of equation (3.34): 3.53
4. Emergent dynamics
Figure 4 presents a representative numerical solution of equation (2.5) starting from a thermodynamically unstable facet ϑ=0 subject to a small initial random perturbation, which models the spinodal decomposition of a crystal facet [1,23]. After an initial transient instability, which is naturally characterized through a linear instability analysis of equation (2.5) around the unstable interface ϑ=0, one observes the emergence of a hill valley-type profile. It comprises O(1)-length segments with an orientation approximated by the Wulff angles ±ω to leading order. These essentially linear O(1)-length segments (psuedo-facets) are seen to meet and merge through O(ε)-length rounded profiles which smoothly interpolate between the disparate angles ±ω of the joining segments: this is consistent with the behaviour of solutions to approximate theories that have been derived from a small-angle assumption [34–36]. We also note that the characteristic length of the interface grows (coarsens) in time as psuedo-facets are extinguished. Furthermore, our numerical data suggest that there is a unique coarsening mechanism, namely a concave neighbouring pair of psuedo-facets are extinguished as their lengths simultaneously shrink to zero. We figuratively describe this unique coarsening motif as a vanishing hillock. Last, we see that the leading-order kinematics of our computed solution to equation (2.5), upon factoring out the trivial advection term , namely , reveals a slowly evolving Wulff-faceted interface .
To comprehend the slow dynamics of , we look to develop smooth ε-families of profiles , i=0,1,…, that are framed on and approximately solve equation (2.5) up to a uniform error O(εi+1); i+1 being the residual order of the family . We will find that the construction of an ε-family of approximate solutions to equation (2.5) with the second-order residual O(ε2) necessarily imposes an intrinsic dynamics on , which we infer as the emergent dynamics of equation (2.5). Before proceeding to do so we need to fix some preliminary notation.
Let denote a Wulff-faceted profile with O(1)-length facets evolving in the slow time τ=εt, and let denote its arc-length s parametrization. We sequentially order its facets, between critical events, in the direction of increasing arc-length, adopting the convention that the ith facet will have angle (−1)iω, and denoting the arc-length of the ith facet's left and right edges at time τ by si−1(τ) and si(τ), respectively. It is notationally convenient to introduce the arc-length s parametrization of the indicator function for the ith facet, namely where and 4.1 We may thereby express the interfacial angle θ0 of through 4.2 while its normal velocity assumes the form 4.3 where υi denotes the normal velocity of its ith segment (facet) (figure 5).
Consider the ε-families of profiles whose parametrizations are obtained through normally deforming according to where , and ψ(ξ,τ) is a given normal deformation. We look to choose ψ(ξ,τ) so that the advected front, , which is naturally parametrized by 4.4 satisfies equation (2.5) to within a piecewise O(ε2) residual.
A direct calculation shows that the arc-length and parametric derivatives of necessarily coincide to leading order: 4.5 Now letting denote the local inclination of , one then finds by direct calculation upon equation (4.4) that 4.6 where θ1:=∂ξψ(ξ,τ). We thereby deduce the following ε-expansion of the curvature of , namely: Also, for the normal velocity of , here denoted by , one readily verifies Thus, satisfies equation (2.5) through order ε provided on each interval , 4.7
To construct from a family of uniformly valid approximate solutions to equation (2.5), here denoted , one needs to asymptotically connect the O(1)-length (outer scale) pseudo-facets of the profile through the O(ε)-length-scale (inner scale) translating fronts (3.32). Note that the slowly evolving part of , here denoted by , is obtained by factoring out the trivial advection 4.8 Let denote the angular profile of the multiple-scale, slowly evolving profile . Sufficiently far away from the vertices of , namely |s−si|≫ε, must reflect the outer-scale angular structure of as expressed in equation (4.6), and so 4.9 But in the vicinity of the vertices of , s∼si, where the outer structures of merge, the angular profile of must also be compatible with those of the respective non-equilibrium translating fronts of equation (3.32), with corresponding angles (equation (3.31)) and so 4.10 where #i denotes − or + depending on whether i is even or odd. Comparing the inner and outer expansions, namely equations (4.10) and (4.9), on an intermediate scale, ζ=ε1/2(s−si), while recalling that equation (3.53) provides the exact O(ε2)-expansion of , we arrive at 4.11 Now matching the O(ε) terms here yields , which when combined with equation (4.7) yields the family of two-point boundary value problems exhibited in figure 6. It is solvable for θ1 only if the following compatibility condition for the facet velocities of is met: 4.12 where li:=si−si−1 (the length of the ith facet of ), and the Peclet length . We have thus theoretically discovered the intrinsic facet dynamics which governs between critical events.
In figure 7, we compare numerically computed periodic solutions of equation (2.5) with matched asymptotic profiles, and (4.8), that uniformly approximate solutions to orders O(ε) and O(ε2), respectively. We note that coincides with the thermodynamic fronts previously identified in , while necessarily encodes both the thermokinetic orientational offsets of the translating fronts, exhibited in figure 3 and encoded in equation (3.53), and the concomitant distortion these induce on the intervening ground state orientations θ=±ω.
5. Scaling laws
Numerical simulation of the facet dynamics prescribed by equation (4.12) suggests that the generic singularity (critical event) of any solution involves a concave pair of facets whose lengths contract to zero in concert before simultaneously vanishing at a common point at some critical time t⋆. That this behaviour mimics the vanishing-hillock coarsening motif observed in the numerical simulations of the original dynamics given in equation (2.5), as exhibited in figure 4, not only explains this motif, but provides us with the means to infer how to resolve the facet dynamics at such a critical time t⋆. Namely, we generate the daughter faceted interface by simply introducing a new vertex at the point into which the vanishing facet pair of the parent interface disappears, thereby joining the two adjacent, and now neighbouring, facets. Possessing two fewer facets than its parent, this daughter interface is thereby morphologically coarser. By continuing to evolve this daughter interface according to equation (4.12), while successively incorporating further such coarsening events as they arise, the average facet length 〈l(t)〉 of will continually grow in time; figure 8 exhibits a representative computed solution. This coarsening facet dynamics faithfully models the large-scale features of solutions to the original multi-scale PDE given by equation (2.5).
We now turn to a study of the length statistics of solutions to our derived coarsening facet dynamics. Starting from Wulff-faceted interfaces with sub-Peclet initial lengths, i.e. 0<li(0)≪Lp, our numerical studies show that, after a short transient, the average facet length 〈l(t)〉 displays the universal scaling law 〈l(t)〉∼t1/2. This scaling persists until 〈l(t)〉 enters the Peclet regime, 〈l(t)〉∼Lp, wherein a crossover to logarithmically slow coarsening ensues thereafter. Furthermore, within this sub-Peclet scaling regime, the distribution of the lengths of the solution is found to be scale invariant (figure 8).
Given our indexing convention is such that the ith facet of has slope (−1)iω, it follows that the purely kinematic relationship between the facet length li of and its neighbouring normal velocities assumes the form By combining this with our derived facet dynamics (4.12), we conclude that the Peclet-scaled lengths evolve in the re-scaled time according to the intrinsic dynamics 5.1 In the sub-Peclet regime, this dynamics for the lengths is effectively given by 5.2 Note that equation (5.2) is invariant under the 1-parameter family of spatio-temporal scalings given by We refer to this parabolic scaling as an emergent symmetry since it is not present in the original dynamics equation (2.5) with 0<ε≪1, but rather emerges as ε↘0 within its effective dynamics as expressed by equation (5.2).
Following an argument appearing in , we may conclude from equation (5.2) that the essentially sole critical event of solutions to equation (5.2) will involve concave pairs of lengths simultaneously shrinking to zero at a corresponding critical time ; namely This explains the origin of the vanishing-hillock coarsening motif previously noted in figure 4, and thereby provides a natural in situ experimental test of theory.
Given the universal statistical behaviour of the distribution of facet lengths that we have already empirically observed, it is natural to presume universal behaviour on the expected (average) length of any solution to equation (5.2). So in particular, we may safely assume that any given solution of equation (5.2) whose lengths at the initial time τ=0 are uniformly small, , will adopt, after an appropriate temporal transient, a universal form , i.e. 5.3 Note that, having focused attention on solutions with a smallness condition at τ=0, we have adopted a distinguished, albeit arbitrary, choice of initial time from which to start to observe the evolution of the distribution of lengths.
The parabolic scaling acts on a solution of equation (5.2) according to 5.4 being a symmetry of the facet dynamics (equation (5.2)), it follows that is also a solution of equation (5.2). We may therefore conclude from the universality hypothesis (equation (5.3)) that 5.5 Setting in the parabolic-equivariance relation that is equation (5.5), one thereby deduces the power law 5.6 This explains the empirically observed scaling appearing in figure 8. Empirical data from direct simulations suggest that the universal constant .
Note that to develop a mean-field theory to predict the structure of the empirically determined scaling function exhibited in figure 8 would require an appropriate interpolation between the famous Lifshitz–Slyozov–Wagner (LSW) theory  for the small-phase fraction limit of Ostwald ripening, and the theory of von Smoluchowski  for the time evolution of the number density of aggregating particles: mean-field theories incorporating such interpolated transport and coagulation mechanisms for faceted-interface dynamics have already begun to be explored , and we hope to develop in future a comprehensive theory to predict the length distribution appearing in figure 8.
In conclusion, we have analytically uncovered a parabolic scaling symmetry that emerges from driven faceting crystal–melt interfaces governed by equation (2.5) through a novel matched-asymptotics analysis. Key to our analysis is the discovery of two non-equilibrium defects—convex and concave translating fronts, which display asymptotic orientations and that deviate from the thermodynamic (ground state) Wulff angles ±ω. We thereby predict that while the characteristic facet length of such crystal–melt interfaces is below an emergent Peclet length scale Lp then the temporal power law ensues and the length statistics of the interface's facets assume a scale-invariant form. We anticipate that the further development of the analytical methods presented herein may prove fruitful in efforts to theoretically characterize the emergent dynamics of driven phase-ordering systems with generic ground state manifolds [41–43].
- Received July 22, 2014.
- Accepted November 14, 2014.
- © 2014 The Author(s) Published by the Royal Society. All rights reserved.