Generic equations for pattern formation in evolving interfaces

We present a general formalism which allows us to derive the evolution equations describing one-dimensional (1D) and isotropic 2D interfacelike systems, that is based on symmetries, conservation laws, multiple scale arguments, and exploits the relevance of coarsening dynamics. Our approach becomes especially significant in the presence of surface morphological instabilities and allows us to classify the most relevant nonlinear terms in the continuum description of these systems. The formalism applies to systems ranging from eroded nanostructures to macroscopic pattern formation. In particular, we show the validity of the theory for novel experiments on ion plasma erosion.


Introduction
Many surfaces and interfaces, which are very different in nature and that occur at very different lengthscales, produce interfaces which are startlingly similar. This similarity originates in the competition between stabilizing and destabilizing physical mechanisms, that are insensitive to the specific value of the interface height h(x, t) (this function providing the height of a interface above position x on a 1D substrate, at time t). Namely, such mechanisms are invariant under global height translations of the form h(x, t) → h(x, t) + ξ, with ξ an arbitrary constant (shift simmetry [1]). For instance, interesting [(sub)micrometric] surface features develop by growth or erosion using various techniques, such as electrochemical deposition (ECD) [2], chemical vapor deposition (CVD) [3], ion beam sputtering (IBS) [4] or molecular beam epitaxy (MBE) [5], but remarkably also in macroscopic systems, such as aeolian sand dunes [6] or underwater (vortex) ripples in sand [7].
For cases in which dynamical instabilities are irrelevant or absent, the shift symmetry has been argued to lead, in the presence of noise, to scale invariant surface morphologies. This is the generic scale invariance [8] associated with the universality classes [9] of kinetic roughening. However, the situation differs in the presence of morphological instabilities, that are the landmark of pattern formation. In these cases, if the height field is seen as a measure of the amplitude of perturbations around a reference homogeneous state, the conservation law associated with the shift symmetry leads to a large-scale instability [10]. This might prevent a unified description by means of a universal equation, such as the Ginzburg-Landau equation is for the case of shortwavelength instabilities [10,11]. Thus, although the systems mentioned above are seemingly governed by similar equations, a general theory of evolving interfaces in which patterns arise is still lacking. In particular, it seems that for every experimental system a specific theory needs to be developed.
In this paper, we present a general formalism for pattern formation in 1D and isotropic 2D surfaces, that relates with previous approaches in the theory of dynamic scaling of rough interfaces [1,12,13], exploiting systematically the formulation of the interface evolution within the assumption of a (generalized) relaxational dynamics. Through symmetries and multiple scale arguments, our approach helps to understand why some effective equations presenting a morphological instability appear almost ubiquitously in the experimental systems mentioned above. Our method enables classification of the most relevant nonlinear terms for these systems, and stresses the relevance of the phenomenon of coarsening to their continuum description. We have also performed some experiments on ion plasma erosion to explore non-trivial implications of our theory. Moreover, we have applied it to other relevant and diverse fields, ranging from thin film production to surface patterns in macroscopic systems, such as fluid waves or some granular systems.
In principle, we consider that the evolution of h is local and, mathematically, takes the form ∂ t h = G({h}) + η, where {h} stands for h and its spatial derivatives, and η is a noise term that we will neglect hereafter. Previous phenomenological proposals for the functional G have been built upon geometry arguments [12], or upon symmetry and conservation requirements [9]. In the first case, the evolution is governed by the local curvature and its derivatives along the interface. This excludes contributions to G that are due to, e.g., anisotropies [14]. Likewise, the second approach can be misleading. For instance, if one considers non-conserved interface dynamics (namely, if G cannot be derived from a current) under the shift symmetry, one could naively drop terms like 1 2 ∂ 2 x h 2 (because it can be written as ∂ x j) or h∂ 2 x h (because it is not shift invariant). However, G could contain a term proportional to (∂ x h) 2 , namely, the sum of both terms. This maybe innocuous in the presence of generic scale invariance, but not in the presence of instabilities, hence a systematic method to classify the terms in the equation is needed.

Non-conserved dynamics
We proceed on by studying separately systems with non-conserved and conserved dynamics. As in Refs. [1,13], we assume that the dynamics of the systems is relaxational in a generalized sense (namely, with a non-constant mobility). Although this assumption is not strictly necessary, it clarifies greatly the interpretation of the coefficients of the non-linear terms below. Hence, where Γ (≥ 0) and F are, respectively, a function and a functional of the height field and its spatial derivatives. Hereafter, we restrict {h} to be {h, ∂ x h, ∂ 2 x h, ∂ 3 x h} [15]. If the relaxation rate Γ is a constant, symmetry under global shifts forbids any explicit dependence of F on h, otherwise, what needs to be independent of h is the right hand side of Eq. (1). Consequently, one can write (note our difference in sign convention from [1]): where the height scale s is the generator of the group of symmetry under translations in h, and we have written Close to the instability threshold, we can assume that local slopes are small, |∂ x h| ≪ 1, and we can expand Γ 1 and F 1 in power series. Therefore, Using these expansions in Eq. (1), we obtain a list of terms compatible with the shift symmetry. The resulting equations are quite involved. Hence, we first consider the linear terms. Later we will include the most relevant non-linear terms. Thus, we find where and γ ≡ γ (000) , f ≡ f (000) . The real part of the dispersion relation for single mode solutions of Eq.
a morphological instability arises if ν < 0 and K > 0. The characteristic lengthscale for the ensuing pattern will be roughly given by the wavelength of the mode maximizing ω q , namely, l = 2π 2K/|ν|. When nonlinear terms are incorporated into Eq. (6), this length-scale can grow in a coarsening process [16], that will help us to identify which nonlinear contributions are actually playing a role in the full dynamics. Close to the instability threshold ν = 0, Eq. (8) can be written as with ǫ a small dimensionless parameter, and (010) ). Note the length-scale l diverges as ǫ −1/2 in the limit ǫ → 0, as K is s-independent, according to Eq. (9). Consequently, we can rescale length, time and height as x → x/l, t → t/l z and h → h/l α , and, using the multiple scales method, look for the most relevant nonlinear term. Similarly to the theory of dynamic scaling [9], we obtain it to be λ(s, where (s, s 2 ) is used to stress the polynomial dependence of λ on s. For many experimental situations, however, the hydrodynamic limit may lie beyond practical reach [17,18], in which case the final structure depends crucially on subdominant terms, that must be included in the evolution equation in order to capture essential features of the original nonlinear system. This can be done within our multiple scale framework. For instance, the ratio between (∂ x h) 2 and (∂ x h)(∂ 2 x h) scales as: ∼ l so, in the limit l → ∞, the first one dominates the second one. This calculation be reproduced for every term resulting in the series expansion, with the results that follow for the first correcting terms, in order of relevance: and we have used that (∂ 2 proposed [19] as the generic amplitude equation for weakly non-linear waves, appearing in contexts from Marangoni or Rayleigh convection to periodic waves in binary fluid mixtures (see Refs. in [19]). In our surface context, it describes e.g. dynamics of (1D) ion-sputtered interfaces under oblique incidence, see Refs. in [20]. In general, while the term (∂ x h) 2 dominates asymptotically, it interrupts coarsening [21,22]. This is actually the case for Eq. (16). Consequently, eventual coarsening of the pattern must be due to other terms in Eqs. (10). E.g., if one adds the λ 2 term in Eq. (10) to the right hand side of Eq. (16), (interrupted) coarsening indeed occurs [23]. Further details and examples are provided below in Section 3.1.

Conserved dynamics
In the case of conserved dynamics for the interface, we can write [1] The assumption of (generalized) relaxational dynamics is now stronger than in the nonconserved one, the operator between Γ and δF/δh constraining the form of the possible terms in a multiple scale expansion (see below). In order to make the problem as general as possible, we need to include the non-relaxational current, j nr , into the right hand side of Eq. (17). The origin of this type of current has been widely discussed e.g. in the context of mound formation in homoepitaxial growth [24,5]. Assuming j nr to depend only on derivatives of h in the form and using Eqs. (4)-(5) as above, Eq. (17) takes the form where all the coefficients in Eq. (19) depend on γ (ijk) , f (ijk) and J (ijk) . Thus, An interesting result derived from the above equations is that, whenever j nr = 0, the condition for unstable growth reads ǫ ∝ s 2 so that power counting (with l ∼ ǫ −1/2 , t ∼ l z and h ∼ l α ) provides terms which give a divergent scaling of the local slope, namely, |∂ x h| scales with a positive power of l. This means that the small slope approximation breaks down, and the expansions above are no longer valid, in analogy e.g. with the so-called superrough scaling in the context of surface kinetic roughening [9], occurring specifically in the case of conserved dynamics. In these conditions, detailed knowledge of system specifics is needed and a strongly nonlinear analysis must be performed in order to obtain the correct interface equation, as in [25].

Results
Thus far, we have derived the most general and relevant interface equations compatible with the shift symmetry. We proceed on by testing the validity of the theory and, more importantly, its capability to identify the right terms of the evolution equation in each context, based on experimental or theoretical evidences, through our procedure that accounts for conservation, symmetries and coarsening, in systems where a surface instability is present. To illustrate this method of analysis, we discuss several 1D examples of systems with lengthscales ranging from microscopic to macroscopic, for each type of dynamics, both non-conserved and conserved. Moreover, some of the examples correspon to isotropic 2D surfaces.

Application to one-dimensional systems
To begin with, we consider non-conserved growth of nanometric sized patterns by Ion Beam Sputtering (IBS). In this technique, the surface of a solid target is bombarded isotropically with a flux of energetic ions. Experimentally, a typical length scale develops at early times, that is seen to grow later (coarsening). The system is not conserved due to the mass loss as a consequence of the bombarding process. Finally, the experimentally obtained surfaces do not develop steep slopes (so a small slope approximation is well justified). With these ingredients in mind, our theory predicts that the evolution equation for these systems must contain the term proportional to λ due to the nonconserved nature of the experiment; and terms λ 2 , λ 101 and λ 210 due to the presence of coarsening. Hence, Consistently with the above expectation, a hydrodynamic theory of erosion has been reported [21] which couples a fast field describing the density of mobile particles at the surface, with a slow one which corresponds to the observed experimental surface. This theory allows to relate experimental magnitudes with corresponding equation parameters and leads to, precisely, Eq. (26). At intermediate length scales (microns to millimeters), we can take growth by ECD [2], where the surface of an aggregate grows by incorporation of cations from a solution. The assumption of locality above is reasonable for a small sticking probability, i.e., if the cations do not attach to the first visited surface site but are, rather, allowed to diffuse further until aggregation occurs. The cation flux at any surface point is isotropic, hence the system is reflection (x → −x) symmetric. Finally, the characteristic length of the unstable structures observed experimentally does not grow in time, i.e., there is no coarsening. These data suggest that surface dynamics for ECD is described by the Kuramoto-Sivashinsky (KS) equation, namely, V = β = 0, due to the reflection symmetry, and that all the terms in Eqs. (10) are negligible, compatible with the absence of coarsening. The same scenario is observed in surface growth by CVD [3]. Indeed, the KS equation has been derived explicitly for both ECD and CVD from constitutive equations, in agreement with our results [17].
Our final example of non-conserved interface growth is taken from surface instabilities in fluids due to surface tension gradients (see [26] for a comprehensive reference), an standard example being provided e.g. by Marangini-Bénard convection [27]. In this system, a thin fluid layer is subject to an external transverse temperature gradient, there being non-conserved fluxes that implement the dissipation in the system. When heating from below a convective instability sets in whose spatial characteristics do not coarsen in time. This instability leads to travelling waves on the free surface of the system that break the x → −x reflection symmetry. An asymptotic study of the corresponding hydrodynamic problem leads to the equation [27,26] in agreement with the general expectation from our approach under the constraints considered.
As an important example of conserved interface dynamics at nanometric scales, we consider the growth by Molecular Beam Epitaxy (MBE) of surfaces that are vicinal to a high symmetry crystalline orientation [28]. In this type of experiments, the molecules incorporated from the chamber arrive isotropically at the surface (and hence reflection symmetry is fullfilled), the relevant one-dimensional interface being given by the location of the step separating one terrace from an adjacent one. Under conditions in which the rate of adatom evaporation from the terrace is vanishingly small, the dynamics is conserved. It is precisely in this case that the small slope approximation has been seen to break down [25], as for Eq. (17) for j nr = 0. However, if adatom desorption occurs on the terraces, the system becomes non-conserved and the relevant interfacial description is provided by the Kuramoto-Sivashinsky equation [29], as expected from our above arguments.
Finally, an example of conserved interface growth in a macroscopic system is provided by the formation of macroscopic ripples in aeolian sand dunes [12]. In this case the wind direction breaks the reflection symmetry and, surface dynamics at the sand bed being conserved due to the conservation in the number of sand grains, one would expect Eq. (19) to hold as the relevant continuum equation. Indeed, this has been shown to be the case within the so-called hydrodynamic approach to pattern formation in these systems [12].
The results in this section, as well as some other experimental systems that can be analyzed along similar lines, are summarized in Table 1, where we provide references to the experimental observations and to theoretical derivations compatible with our conclusions. We provide the reader with broader references when connection with specific experiments is less unambiguous.

Application to novel experiments
The examples in last subsection show how, based on general considerations, one can construct a general equation in good agreement with current specific theories. Notwithstanding, we want to emphasize the wide range of generality of our method by applying it to novel experiments on ion plasma erosion. Thus, a Si(100) wafer (2 inches diameter) was immersed in an argon plasma, which was confined magnetically, at a pressure of 5 × 10 −3 mbar. The only parameter that was changed in the experiments was the immersion (i.e. sputtering) time. After the sputtering process, the central part of the wafer was analyzed by Atomic Force Microscopy operating in tapping mode with a silicon cantilever. Under these conditions the flux of incoming ions on the central part of the wafer is mainly isotropic. We plot two snapshots of the surface taken at 3h (Fig. 2a) and 6h (Fig. 2b). Moreover, we have measured a non-zero coarsening exponent n = 0.53 ± 0.02 in ℓ ∼ t n , ℓ being the mean cell diameter. The system being isotropic, we can generalize our results above to this case. Out of the most relevant subdominant terms in Eq. (10) (once generalized to two dimensions) compatible with isotropy, λ 2 [30,22] and λ 210 [31,32] are known in the literature to induce coarsening behavior when taken each of them in combination with λ 100 and a KS dispersion relation. Extensive numerical simulations [23] show that the term λ 101 also induces coarsening, but the term λ 400 does not. Hence, thus far we are left with an equation of the following form: ¿From the experimental data we can extract some quantitative information about the coefficients in Eq. (28). Thus, a closer inspection of the experimental surface cross section ( Fig. 1) reveals that cell shape can be approximated by which is a solution of Eq. (28) in 1D with λ 2 = λ 101 = 0 (and, consequently, approximately valid for large cells in the isotropic 2D case). This suggests that λ 2 and λ 101 are negligible since large values for them would lead, rather, to parabolic cells [21]. Consequently, Eq. (28) reduces effectively to the so-called convective Cahn-Hilliard equation [32]. For short times (linear regime), the typical lengthscale allows determination of the K/ν value, while remaining parameters are estimated by trial and error in order to improve the resemblance with the experimental morphologies. Results of such numerical simulations are compared with experiments in Fig. 2 showing an excellent agreement between them. Thus, the small panels (upper insets) display two snapshots of the numerical integration of Eq. (30) and the large panels the equivalent results from experiments.

Discussion and conclusions
Although table I is not meant to be exhaustive, it leads to some general observations. Thus, as anticipated above, the occurrence of coarsening is an important property in guiding the identification of the appropriate nonlinear terms of a given system. This phenomenon was associated with the seeming inability to describe patterns that are stable both in wavelength and amplitude by means of local height equations [7,16]. However, there may exist now counterexamples [21,22] for this shortcoming, at least for an order range that is much larger than the single cell size. As a consequence, the task of determining the necessary and sufficient conditions for coarsening [16], as well as the possible values (universality classes) for the coarsening exponent n, remains to be completed, lying outside the scope of this paper. Note that, once coarsening occurs, the fact that it interrupts [16] or not, and the value of exponent n, both depend on the type of leading and subleading nonlinearities that appear in the evolution equation. Hence, these provide additional criteria in order to propose a continuum description for a given physical system, although it can nevertheless be the case that two different subleading terms lead to the same coarsening exponent value, such as is the case for λ 2 and λ 210 , both of them yielding n = 1/2 [33,31]. Moreover, system dimensionality and the (possibly related) relevance of noise may influence the value of the coarsening exponent in the corresponding universality classes.
In summary, we have provided a systematic method to obtain and classify the relevant nonlinear terms in interfacial systems which present morphological instabilities and the shift symmetry. Besides, we have applied our method to some reported systems in the literature as well as to novel experiments on ion plasma erosion performed for this purpose.