Table Of Content

Universality classes for unstable crystal growth Sofia Biagi,1,2, Chaouqi Misbah,1,2, and Paolo Politi2,3, ∗ † ‡ 1Université Grenoble 1/CNRS, LIPhy UMR 5588, Grenoble, F-38401, France 2Istituto dei Sistemi Complessi, Consiglio Nazionale delle Ricerche, Via Madonna del Piano 10, 50019 Sesto Fiorentino, Italy 3INFN Sezione di Firenze, via G. Sansone 1, 50019 Sesto Fiorentino, Italy (Dated: June13, 2014) 4 Universality has been a key concept for the classification of equilibrium critical phenomena, 1 allowing associations among different physical processes and models. When dealing with non- 0 equilibrium problems, however, the distinction in universality classes is not as clear and few are 2 the examples, as phase separation and kinetic roughening, for which universality has allowed to n classify results in a general spirit. Here we focus on an out-of-equilibrium case, unstable crystal u growth,lyinginbetweenphaseorderingandpatternformation. Weconsiderawellestablished2+1 J dimensionalfamily ofcontinuumnonlinearequationsforthelocal heighth(x,t)ofacrystalsurface 2 having the general form ∂th(x,t)=−∇·[j(∇h)+∇(∇2h)]: j(∇h) is an arbitrary function, which 1 is linear for small ∇h, and whose structure expresses instabilities which lead to the formation of pyramid-likestructuresofplanarsizeLandheightH. Ourtaskisthechoiceandcalculation ofthe ] quantitiesthatcanoperateascriticalexponents,togetherwiththediscussionofwhatisrelevantor h not to thedefinition of our universality class. These aims are achieved by means of a perturbative, c multiscale analysis of our model, leading to phase diffusion equations whose diffusion coefficients e m encapsulate all relevant informations on dynamics. We identify two critical exponents: i) the coarsening exponent, n, controlling the increase in time of the typical size of the pattern, L ∼ tn; - ii) the exponent β, controlling the increase in time of the typical slope of the pattern, M ∼ tβ t a whereM ≈H/L. Ourstudyrevealsthatthereareonlytwodifferentuniversalityclasses, according t tothepresence(n=1/3, β=0) ortheabsence (n=1/4, β >0) of faceting. Thesymmetry of the s . pattern,aswell as thesymmetryof thesurface mass current j(∇h) andits precisefunctional form, t a is irrelevant. Ouranalysis seems to support theidea that also space dimensionality is irrelevant. m - d PACSnumbers: 05.70.Ln,81.10.Aj,05.45.-a n o c I. INTRODUCTION: UNIVERSALITY CLASSES field). [ When passing to nonequilibrium processes, the phe- 3 The concept of universality is very useful in physics, nomenologyis muchwider anda classificationin univer- v because it allows to classify seemingly different phenom- sality classes is not as firm. A much studied case is the 3 ena and models. Perhaps, one of the oldest examples so-called “phase separation”. Let us consider a system 6 is the universal form of the Van der Waals equation of undergoing a continuous phase transition (at T = T ) c 2 state (law of corresponding states [1]), which is the sim- when passing from a disordered high temperature phase 4 plest equation describing a change of state and valid for to an ordered low temperature phase. If the tempera- . 1 anyfluid. Aclearformalizationofuniversalitywasfirstly ture T is suddenly decreased (quenching) from T > T i c 0 possible for equilibrium critical phenomena, where the to T < T , the system undergoes an ordering process 4 f c renormalization group theory allows to give a rigorous where the typical size L of ordered regions increases in 1 definition of which parameters are relevant (universal) time, L(t). This process, called coarsening, lasts forever : v andwhicharenot. Forexample,withinimportantclasses (for infinite systems) if the system is globallyat thermo- Xi of ferromagnetic spin models, it is known that relevant dynamic equilibrium. In most cases L(t) increases as a parametersare: the physicaldimension of the space, the power law, L(t) tn, which defines the coarsening ex- r ∼ a dimension of the order parameter and its symmetries, ponent n. Generally speaking, T and T are irrelevant i f the (short/long) range of interaction of the coupling. A parametersanditappearsthatthephysicalspacedimen- universalityclass is uniquely defined by its criticalexpo- sion is also irrelevant (as long as T is finite, T > 0). It c c nents, which describe the behaviour of the order param- appears instead that conservation laws are relevant for eterinproximityofthe criticalpointasafunctioneither the dynamics and it is reasonable to expect that a con- ofthecontrolparameter(e.g. the temperature)orofthe servation law slows down the dynamics and reduces the conjugatefieldoftheorderparameter(e.g. themagnetic coarsening exponent, a known fact at present [2]. The spirit of universality also means that the same model is important for different physical problems: for sofi[email protected] example, phase separation and pattern formation have ∗ [email protected] several similarities. Therefore our study, which focuses † ‡ [email protected] onacertainclassofgrowthequationsforcrystalsexhibit- 2 ingpatternformation,isexpectedtoberelevantforboth strong support regarding its irrelevance. fields. Thisclassofequations,seeEq.(3),hasemergedin thelasttwentyyearsasaprototypicaldescriptionofcrys- tal growth by deposition processes. It has some similarities with well known models as the Cahn-Hilliard equa- II. CRYSTAL GROWTH EQUATION tion and the clock models [3], but its general properties havebeen now fully established,as discussedin the next InthisSectionweshallgiveabriefintroductiontothe section. class of equations we are interested in, mainly address- This equation leads to morphological instability of a ing the qualitative aspects of the dynamics rather than planar surface, with formation of mounds/pyramids out their physical derivation (for a thorough discussion on of the flat front. In general, as previously sketched for the physical background, the reader is referred to [6]). domains in phase separation processes, mounds coarsen, A growing planar crystal surface (growing by molecu- but under some conditions we show that other scenarios lar beam epitaxy, for example) can undergo a morpho- take place. In apparent contrast to some existing litera- logical instability resulting into the formation of three- ture(seeSectionVIIA),weareabletostatethatpattern dimensional mounds or pyramids of linear size L and symmetry is irrelevant and only two universality classes height H. The subsequent morphological evolution may result, depending on whether mound’s slope is constant range from a pattern of constant L and an increasing (faceting)oritisanincreasingfunctionoftime. Thislast H up to a perpetual increase of L in the course of feature is known a priori, from visual inspection of the time (coarsening), with H increasing in concert. An in- surfacecurrentj(m)(seebelow)andallowsthedefinition termediate scenario may also take place in some cases, of a second exponent β 0, describing the behavior of where L(t) increases up to a length L reached at a ≥ max thetypicalmound’sslopeM intime,M(t) tβ. Thetwo giventime,beyondwhichthe moundsizeisfrozen,while ∼ universalityclasses we havefound are therefore given,in mound height keeps growing. This scenario corresponds the caseofconstantslope m∗ (j(m∗)=0), by β =0 and to interrupted coarsening[7]. We are not awareof a sce- n = 1/3, and, in the case of increasing slope, by β > 0 nario where both L and H keep constant in time. and n=1/4. In the case of a perpetual coarsening the generic evo- The idea used to establish these results is based on lution law of L(t) is algebraic with coarsening exponent, the statement that coarsening takes place if the steady- n, defined as L(t) tn. During the coarsening process, ∼ state periodic solutions are unstable against perturba- the typical slope M H/L may either keep constant or tions of the phase of the pattern [4]. More precisely, increase in time, M(≈t) tβ, therefore defining a second ∼ a periodic pattern has a constant wavenumber q which exponent β 0. ≥ acquires a space-time dependence when the pattern is Fromamesoscopicpointofview,thelocalvelocity∂ z t perturbed (we can also define the phase of the pattern): of a surface z(x,t) growing under a deposition flux of if the periodic pattern is perturbed, the wavenumber(as intensity F must have the form 0 wellasthephase)willvaryfromonepointtoanother. If the perturbation grows with time, we say that the pat- ∂tz(x,t)=F0 Jtot, (1) −∇· tern is unstable with respect to wavenumber (or phase) fluctuations. If the periodic pattern is unstable with re- provided that the deposited mass on the surface does spect to phase fluctuations then we expect coarsening to not evaporate and that no holes occur in the growing takeplace. Itwillbeshownthatthephaseofthepattern solids [8]. The total current Jtot is a function of the obeysadiffusionequationandinstabilityissignaledbya slope m = z and higher order spatial derivatives and ∇ negativediffusioncoefficientD. Thisdiffusioncoefficient it accounts for all surface rearrangement processes. Its (actually in two dimensions there are several diffusion simplest form is coefficients, as we shall see) depends on the steady-state patternproperties,andmoreparticularlyonthemodulus Jtot =j(m)+Γ 2m, (2) ∇ of the wavenumber q. By using a dimensional relation, D(q) L2/t, where q = 2π/L, we shall extract the wherej(m)isafunctionoftheslopeonlyanditaccounts | | ≈ fortheexistenceofamasscurrentonaterrace. Atsmall coarsening exponent. slopes j ν z: if the current is uphill (ν > 0), the flat Here we are able to make stronger and more general ≈ ∇ surfaceisdestabilizedatsufficientlylargescales. Thesec- statements with respect to [5], facing a wider range ondterm, 2m,regularizesthedynamicsatshortlength of two-dimensional patterns and stressing on universal ∇ scales and it may have different physical origins [9]. features of unstable crystal growth. This is the focus By performing the substitution z h = z F t and of the present paper: going beyond the details of the → − 0 afterappropriaterescalingofxandt,itispossibletoab- equation and of the physical process and pointing out sorbΓ and ν into the new variables so that the equation what is relevant, slope selection or not, and what is can be written in the form not, the symmetry of the pattern and that of the mass current. Although no complete proof about physical ∂h(x,t) space dimensionality is accomplished, our study reveals = [j( h)+ ( 2h)] Jtot, (3) ∂t −∇· ∇ ∇ ∇ ≡−∇· 3 where j( h)= h+ higher order terms. Interestingnonlineardynamicsappearslaterandperi- ∇ ∇ Someimportantfeaturesofthenonlineardynamicscan odic steady-state solutions play the major role, because be discussed by referring to the one dimensional version relevant informations can be drawn from their stability ofEq.(3),whichhasbeendiscussedatlengthinRef.[4]: The general idea used here is that if coarsening takes place, this means that every steady-state solution is un- ∂th= ∂x[j(hx)+hxxx]. (4) stablewithrespecttowavelengthfluctuationsandthere- − fore the relevant variable to describe this phenomenon In fact, by taking the spatialderivativeof both sides, we is the wavelength, or, more precisely, the phase of the get the generalized Cahn-Hilliard equation, pattern, since in nonlinear systems it is known that the phase is a more appropriate variable to deal with rather ∂ m= ∂ [j(m)+m ] (5) t − xx xx than the wavelength itself [11]. This idea was applied with success to study one dimensional fronts in Ref. [4], where the shape of the potential U(m) = dmj(m) de- wheretheabilityofthesystemtodevelopcoarseningwas termines the type of dynamics [10]: (i) stationary solu- R directly related to steady-state properties, with no need tions, satisfying j(m)+m =0, correspondto periodic xx to perform a forwardtime-dependent calculation. It will “oscillations”withinthepotentialwellU(m);(ii)thereis evenbeshownforseveralexamplesbelow,thatthestabil- coarsening if and only if the wavelength of such station- ityorinstabilityofthepatternagainstphasefluctuations arysolutionsisanincreasingfunctionoftheiramplitude; can be concluded analytically. Even more importantly, in general three scenarios, depicted here above, are pos- our approachprovides the values of exponents n and β. sible: perpetual coarsening, interrupted coarsening, no In order to study stability of the periodic steady-state coarsening; (iii) the slope of emerging mounds is con- h , we seek for solutions ofthe nonlinearequationin the stant if U(m) has maxima at finite m= m , otherwise 0 ∗ ± form (with ε small parameter) slope increases forever. When passing from one to two dimensions, i.e. pass- h=h˜ +εh˜ +..., (7) ing from Eq. (4) to Eq. (3), the equivalence between the 0 1 growth equation and the Cahn-Hilliard equation ceases to be valid [6] (see also Section VII). Furthermore, the and linearize the equation. However, in addition to this surfacecurrentjrequiresspecificationofitsin-planesym- quite standard study of linear stability, the crux of our metry, which adds a new degree of freedom to the prob- method is to introduce a multiscale analysisthat will al- lem. The following Sections present the various dynami- low us to extractthe phase evolutionequation,the anal- cal scenarios where the values of the exponents n and β ysisofwhichwillinformusonthepresenceofcoarsening (universalityclasses)areextractedforthefamilyofmod- or the lack thereof. Therefore, besides the fast variables els defined by Eq. (3). We shall follow a multiscale per- x and t, we introduce slow variables defined as turbative approach,discussed in the next Section, which allowsustowritedownthephasediffusionequationthat X=εx, T =ε2t. (8) describes the evolution of the typical mound size in the courseoftime. Thevariousdynamicalscenarioswillcon- The perturbation parameter ε is a small quantity that stitute the subject of Sec. IV while distinct universality defines the fact that we are looking for long wavelength classes will presented in Sec. VI. A thorough discussion modulation of the pattern, which are the most “danger- of our results will follow in Sec. VII. ous”modes (see [4] for more details). In a multiscale spirit fast and slow variables are treated as if they were independent [12]. As already said, it is convenient to III. THE PHASE DIFFUSION EQUATION work with the phase variables rather than with the spatial variables. For that purpose, we introduce (in two dimensions) two scalarphase variablesϕ andϕ . If the AsalreadyanticipatedinthepreviousSection,theflat 1 2 patternisperfectlyperiodicthenthesevariablesaresim- profile, namely the solution h 0 of Eq. (3), is unsta- ≡ ply given by ble. This is easily shown from a linear stability analysis: setting h = δ exp(ωt+ik x) in Eq. (3) and assuming · ϕ :=q x, ϕ :=q x, (9) δ 1, we obtain the linear spectrum: 1 1 2 2 · · ≪ ω(k)=k2 k4, (6) where qi = ϕi are the basis wave vectors defining the − symmetryof∇thestationaryperiodicpattern. Toaccount where k = k. This result shows that there is a band for perturbations of the periodic lattice, q-vectors are | | of wave-vectors (0 < k < 1) with positive ω, so the cor- notjustconstantsbut havea dependance onslowscales: responding harmonic amplitude increases exponentially q = q(T,X); therefore we introduce for convenience the with time until nonlinearities can no longer be disre- slow phase scales: ψi = εϕi, so that qi = Xψi can be ∇ garded. This instability will result first in a deformed expressed as function of slow variables only. (moreorlessregular)surfaceandduringtheinitialstages According to this approach, various differential oper- the amplitude grows quite rapidly. ators in the model equation have to be substituted as 4 follows: that the right hand side of Eq. (15) is orthogonal to the kernel of the adjoint operator of ) are verified. These ∂ ε2∂ = ε2[(∂ ϕ )∂ +(∂ ϕ )∂ ] L t → T T 1 ϕ1 T 2 ϕ2 conditions have the following form [14]: = ε[(∂ ψ )∂ +(∂ ψ )∂ ], (10) T 1 ϕ1 T 2 ϕ2 ∇→∇0+ε∇X, (11) hvi,gi=0, (20) with ∇0 = q1∂ϕ1 +q2∂ϕ2 and ∇X = (∂X,∂Y). Then where functions v1, v2 verify L†[v] = 0. We therefore expansions(7),(10)and(11)arereportedintothemodel calculate the adjoint of our linear operator from the † L equation (3) which yields (by keeping only terms up to definition v,u = v, u . Given that † hL i h L i order ε, see Appendix A) [u]= u 2( 2u), (21) ε[(∂ ψ )∂ h˜ +(∂ ψ )∂ h˜ ]= L −∇0·J∇0 −∇0 ∇0 T 1 ϕ1 0 T 2 ϕ2 0 = ( 0+ε X) j( 0h˜0)+ε ( 0h˜1+ Xh˜0) L is self-adjoint if and only if the Jacobian matrix J − ∇ ∇ ·{ ∇ J ∇ ∇ is symmetric (see Appendix B). This latter case is def- +∇0(∇20˜h0)+ε[∇X(∇20˜h0)+∇0(∇21h˜0+∇20h˜1)]}, initely the most common one, since we find that it is (12) assured by all the explicit forms of j used in the litera- ture. We also stress that a symmetric means that the to be studied order by order. current derives from a potential, J J= δ /δm (see tot Zeroth order – The zeroth-order defines stationary so- − F SectionVII for further details). lutions ˜h as the unperturbed ones: 0 For the sake of completeness, we must keep in mind 0=∇0·[j(∇0h˜0)+∇0(∇20h˜0)]=∇0·(J0)tot ≡N[h˜0], cthaabtilitthyeosfymthmeecturyrrepnrotpmerettyhofodr:Jaipshnaosteadliiffmuistioinnaepqpulai-- (13) tion could be derived formally without having a linear where is a nonlinear operator acting on h˜ . Explicit solutionNs ˜h are in general not available, the0only basic self-adjoint operator. However, for non adjoint opera- 0 tors, the solutions of [v] = 0 can be obtained, in gen- information being that h˜ enjoys periodicity properties L† 0 eral, only numerically [15], even if examples to get them in ϕ , ϕ . Focusing our analysis on high symmetry sub- 1 2 analytically in 1d do exist [4]. strates, for which h = 0, a stronger condition can be h i According to the above discussion, if is symmetric imposed: than = and v = ∂ h˜ (ϕ ). It is nJow possible to L L† i ϕi 0 i rewrite g, see Eq. (16), as follows (see Appendix C for (J ) =0. (14) 0 tot more details): First order – At first order we obtain a linear and in- homogeneous equation for h˜1: g ≡(∂Tψ1)∂1h0+(∂Tψ2)∂2h0−(ψα)βγcαβγ, (22) L[h˜1]=g(h˜0,ψ1,ψ2), (15) where ∂α ≡∂ϕα and h0 ≡h˜0 for ease of notation, and where ∂˜h cα =q ∂ 0 +2q q ∂ ∂ ∂ h˜ g ≡(∂Tψ1)∂ϕ1h˜0+(∂Tψ2)∂ϕ2h˜0+∇0·[J(∇Xh˜0) − βγ δν δ"Jνγ∂qαβ# jγ lβ α l j 0 +∇X(∇20h˜0)+∇0(∇21h˜0)] (16) +3 2q ∂ ∂˜h0 +δ 2∂ h˜ , (23) is a function of stationary solutions h˜ , while ∇0 νβ ν∂qαγ βγ∇0 α 0 0 [h˜ ] [ ( h˜ )+ ( 2h˜ )] (17) with qij =(qi)j as the j-th component of the i-th wave- L 1 ≡−∇0· J ∇0 1 ∇0 ∇0 1 vector. Moreover, the compact notation h ∂ ˜h will j ≡ ϕj 0 is the Fréchet derivative of , defined as be adopted from now on. N Byusingtheaboveexpressionforginthetwosolvabil- [h˜ +εh˜ ]= [h˜ ]+ε [h˜ ]. (18) 0 1 0 0 1 ity conditions (20) we obtain the phase diffusion equa- N N L tions (i=1,2): By virtue of translational invariance of with respect to space variables, it follows that [h˜ (Nϕ +∆ )] must 0 i i ∂ψ vanish as well. In the limit ∆i →0Nwe get ∂Tψi = ∂Xβ∂αXγD˜βiαγ, α,β,γ =1,2 (24) [h˜ (ϕ +∆ )]= [h˜ (ϕ )]+∆ [∂ h˜ (ϕ )]=0, (19) N 0 i i N 0 i iL ϕi 0 i where repeatedindices are to be summed overaccording which alsoimplies that [∂ h˜ (ϕ)]=0. Therefore,since to Einstein’s convention. The diffusion coefficients have ϕ 0 [h˜ ] = 0 has nontrivialLsolutions (∂ h˜ ), the Fredholm the following expressions: 1 ϕ 0 L alternative theorem [13] can be used for Eq.(15). Such h ,cα h ,h h ,cα h ,h theorem guarantees solutions for Eq. (15) if and only if D˜1α = h 1 βγih 2 2i−h 2 βγih 1 2i (25) the so called solvability conditions (expressing the fact βγ " h1,h1 h2,h2 h1,h2 2 # h ih i−h i 5 symmetry of h0(x) Θ p invariances Dβiαγ independentDβiαγsisinfacttwelve;forrhombicandrect- oblique no specific no specific 2-fold 12 angular ones parity allows to reduce this number to six; rhombic no specific 1 2-fold, Π 6 2 then,theincreaseddegreeintherotationalinvariancefor rectangular π/2 no specific 2-fold, Π 6 2 thesquareandhexagonalcasesimplies furtherreduction square π/2 1 4-fold, Π 3 2 to, respectively, three and two independent coefficients. hexagonal π/3 1 6-fold, Π 2 2 The 3-fold symmetry shares similarities with the hexag- triangular 2π/3 1 3-fold, Π 2 1 onal pattern (albeit the two symmetries are distinct). It TABLEI.Presentation of the5 twodimensional Bravais lat- turns out that these two symmetries obey formally the tices,classified accordingtorelativeorientationsbetweenthe same diffusion equation, with the same number of inde- twoq-vectors(Θ)andtheirrelativeamplitude(p);moreover, pendent coefficients. paritysymmetry(Π2withrespecttobothspacevariables,Π1 In the next Section we are going to exploit the phase with respect toasingle space variable) androtational invari- diffusion equations for some symmetries and we will re- ances are specified. The last, extra, row refers to the 3-fold port on some far-reaching consequences. In particu- case. An increasing symmetry corresponds to a decreasing lar, we will examine stability of Eq.(24) with respect to numberof independentdiffusion coefficients Diα. βγ phaseperturbations,arelevantinformationregardingthe coarsening problem. andD˜β2αγ 1↔=2D˜β1αγ. Itisconvenienttodefinenewdiffusion coefficients Diα, by regrouping similar derivatives: βγ IV. THE COARSENING CONDITIONS D˜iα β =γ Diα = βγ (26) A coarsening dynamics is signaled by phase instabil- βγ D˜iα +D˜iα β =γ. ( βγ γβ 6 ity, i.e. by a phase which increases exponentially with time[17]. Phasediffusionequations,Eqs.(24),arelinear Therefore,inthemostgeneralcase,Eqs.(24)havetwelve and can be solved assuming independent diffusion coefficients. Their expressions are in generalquite involved except if some symmetry prop- ψ (X,T)=ψ(0)exp(ΩT)exp(iK X) (28) erties of the steady-state solutions h are evoked. Sym- 1,2 1,2 · 0 metry properties will lower the number of independent and imposing a null determinant for the linear system diffusion coefficients. It should be remembered that h0 with unknowns ψ(0). This way, we can write down a is a perfectly periodic in-plane pattern, defined by one 1,2 quadratic equation for Ω of the five known two-dimensional Bravais lattices. Se- lecting one of these patterns for the stationary solution Ω2+f(Diα,K)Ω+g(Diα,K)=0 (29) h means fixing the two q-vectors and the space group βγ βγ 0 symmetry that leave h0 unchanged. It is convenient to and obtain two entire spectra, Ω1,2 = Ω1,2(K), whose list the Bravais lattices in a sort of hierarchy to face at properties depend on the symmetry of h . We present 0 once howthe demandofsymmetry simplifies the expres- herebelowdetailedresultsregardingrectangular,square, sion of diffusion equations. Let us define Θ as the angle hexagonal and triangular symmetries. Appendix E lists betweenthe twoq-vectorsandp,the proportionalitybe- the q-vectors used in these specific cases. The oblique tween their moduli: and rhombic symmetries will not be treated here since they involve quite lengthy expressions. Since we do not q q q 1 2 1 cosΘ:= · , p:= | |. (27) expectanynewspecificityassociatedwiththem(seelater q q q | 1|| 2| | 2| discussion),wedidnotfeelitworthwhiletodwellonthis The proposedorder for the five Bravais lattices is shown issue. inTableI.Inaddition,givenitsconsiderablerelevanceto experiments [16] we also studied the 3-fold case, that is notincludedamongtheBravaislatticesbut,nevertheless, A. The hexagonal and triangular symmetries canbedealtwithusingthesamemethodasfortheother symmetries. The 3-fold case is characterized by Θ = Inthe6-foldandinthe3-foldcasesthespectrumturns 2π/3 and a p=1, while the parity symmetry holds for a outtobeisotropicinKandthetwoeigenvaluesarefound single space variable only. to be: InAppendixDweprovideanexplicittreatmentofthe phasediffusionequationforthehexagonalsymmetryand Ω (K)= D K2, Ω (K)= D K2, (30) 1 22 2 11 − − determine the number of independent coefficients. This serves as a guide for the other symmetries for which we where D1111 ≡ D2222 ≡ D11 and D2121 ≡ D1212 ≡ D22 (see do not reportthe details. Our results are summarizedin Appendix D). Since D22 is positive: Table I: last column reports the number of independent coefficientsDβiαγ correspondingtoeachpatternsymmetry. D = 9q2 h2 >0, (31) For oblique, that is the most general one, the number of 22 h2 h 12i h ϕi 6 the eigenvalue Ω is negative, signaling stability of the AnegativeD wouldsignalinstability. Wewillseelater 1 11 pattern. The other eigenvalue, instead, has no a priori how to determine this sign analytically and how to dis- fixed sign: criminate among different dynamical scenarios. B. Square and rectangular symmetries For these symmetries, the spectrum of eigenvalues is anisotropicandits analysisis, in principle, more compli- 4q7/4 D = ∂ (q5/4 h2 ). (32) cated. Let us first consider the square case spectrum: 11 h2 q h 12i h 1i K2 Ω (K,θ)= (D +D ) (D D )2+4[D2 (D D )2]sin2(θ)+4[(D D )2 D2 ]sin4(θ) , (33) 1,2 − 2 " 11 22 ± 11− 22 12− 11− 22 11− 22 − 12 # q where K =Kcosθ, K =Ksinθ and where we have used the compact notations: D11 D22 D , D11 D22 D , and1D12 D21 2 D . Expression (33) shows that extremal values for Ω (K1,1θ≡) in2t2h≡e (K11,K2)2p≡lane11ar≡e 22 12 ≡ 12 ≡ 1π2 π π 1,2 1 2 along the directions θ = n and θ = +n . Since we are dealing with a 4-fold symmetry, we consider just two 2 4 2 cases, for each of which we distinguish two different eigenvalues: θ =0, Ω0(K)= D K2 , Ω0(K)= D K2; • 1 − 22 2 − 11 θ =π/4, Ωπ/4(K)= (D +D D )K2/2 , Ωπ/4(K)= (D +D +D )K2/2. • 1 − 11 22− 12 2 − 11 22 12 As already seen for hexagonaland triangular symmetries, also in this case one eigenvalue for each couple is always negative, since: q2 D = [ h2 +3 h2 ]>0, (34) 22 h2 h 11i h 12i h ϕi 4q2 D +D D = h2 >0, (35) 11 22− 12 h2 h 11i h ϕi while the sign of the other eigenvalues is not obvious, being determined by that of the following expressions: 1 D = [∂ (q3 h2 )+q3∂ h2 +q2 h2 ], (36) 11 h2 q h 11i qh 12i h 12i h 1i 4 1 1 D +D +D = q3∂ h2 +q2 h2 + q3∂ h2 +2q2 h2 . (37) 11 22 12 h2 2 qh 11i h 11i 2 qh 12i h 12i h 1i(cid:20) (cid:21) In Section V we propose calculations of the diffusion coefficients valid in the weakly nonlinear regime, in order to treat those eigenvalues whose sign has not been easily recognizable. Analogously, for the rectangular case the spectrum takes the following form: K2 Ω (K,s)= (D11+D22)+(D11+D22 D11 D22)s 1,2 − 2 (" 11 11 22 22− 11− 11 # (D11 D22)2+ 2(D11 D22)(D11 D22 D11+D22)+4D12D21 s+ (D11 D22 D11+D22)2 4D12D21 s2 ±vu 11− 11 " 11− 11 22− 22− 11 11 12 12# " 22− 22− 11 11 − 12 12# ) u t (38) wheres sin2(θ). Theextremalvaluesarenowobtained new other directionswe are ableto specify in the weakly ≡ notonlyalongthemaximalsymmetrydirections,namely nonlinear regime (see Appendix F). According to such along θ = nπ and θ = nπ π/2, but also along two approximation, corresponding to steady states of small ± 7 amplitude a, these directions are close to θ =π/4: bothfromthatoftheone-dimensionalproblemandfrom Eq. (44), must be a decreasing function of q. Thus, we sin2(θ)= 1 2√2 q4 p2−1 p 1 O(a2), (39) can state that the nature of the function whose decreas- 2 ± m2(p2+1)2 ≡ 2 ± ing character determines stability depends on the space dimension and on the class of the considered equations where m = q3∂q(a2)/a2 and p is defined by Eq. (27). [20]. The eigenvalues corresponding to the two first extremal The secondimportantconclusionis thatour results in directions are: this sectiondo not depend onthe natureof the currentj entering in Eq. (3). The readercanrefer to Appendix D θ =0 • in order to check this statement in detail for the 6-fold Ω0(K)= D22K2 , Ω01(K)=−D1111K2; symmetry. 2 − 11 θ =π/2 • V. THE DIFFUSION EQUATION IN THE Ωπ/2(K)= D11K2 , 1 − 22 WEAKLY NONLINEAR REGIME Ωπ/2(K)= D22K2. 2 − 22 In this Section our aim is to analyze if coarsening oc- Again, one eigenvaluefor eachcouple is alwaysnegative, curs or not, while the determination of coarsening expo- since: nents will be presented in the next section. In order to q2 determine the dynamical scenarios for our growth equa- D11 = [ h2 +3p2 h2 ]>0, (40) 22 h2 h 11i h 12i tion, we need an evaluation of the signs of appropriate h 1i diffusion coefficients, see Eqs.(D2,36,37,42). This task q2 D22 = [3 h2 +p2 h2 ]>0 (41) can, in general, be performed only numerically by solv- 11 h2 h 12i h 22i h 2i ingforthesteady-statesolutions. However,byrestricting ourselvesto a weakly nonlinearanalysis,someanalytical while the otherhas no obvioussign,asit is fixedby that results can be obtained. To that end we assume that of the following expressions: the amplitude of the stationary solution h is small. We 0 1 havealreadyperformedthegenerallinearanalysisofour D11 = [∂ (q3 h2 )+q3p2∂ h2 +q2p2 h2 ], (42) 11 h2 q h 11i qh 12i h 12i equation, which has resulted into the spectrum (6). In a h 1i weakly nonlinear approachwe can push further this sta- bilityanalysisextractinganapproximatedsolutionforh 0 D22 = 1 [p2∂ (q3 h2 )+q3∂ h2 +q2 h2 ]. (43) in power series of the amplitudes of the Fourier modes. 22 h2 q h 22i qh 12i h 12i Thanks to the periodic character of the stationary solu- h 2i tionwe canexpressh(x,t) with a Fourier seriesthat can It is worth notice that D2121 =D1212 and D1111 =D2222, for be truncated at some order. The small amplitude limit p=1[18]. The other two extremal directions, defined by is legitimate as long as k 1, so that higher harmonics Eq. (39) have to be considered in the weakly nonlinear are stable, ensuring the co→nsistency of the truncation of regime. The reader can find calculations in Appendix F. the series. Here, it suffices to say that, also in the rectangular case, Since the symmetryofthe growingpatternisidentical once the direction has been fixed, the sign of one eigen- orlowerthansubstratesymmetry,anisotropiccurrentis valueisnegativewhilethesignoftheotherisnotevident. the most general one, i.e. compatible with any Bravais As a summary of this section we can highlight two lattice. We consider a generic class of isotropic currents important conclusions. In the hexagonal and triangular j(m,c ,c ) = m(1+c m2 +c m4) so that Eq. (3) be- 2 4 2 4 symmetries,oneeigenvalueis positive(phaseinstability) comes: if the quantity (see Eq. (32)) ∂ h=L˜[h] c [3(h2h +h2h )+h2h +h2h t − 2 x xx y yy y xx x yy A≡q5/4hh212i (44) +4hxhyhxy]−c4[5(h4xhxx+h4yhyy)+h4xhyy+h4yhxx isadecreasingfunctionofthe wavenumberq. Thequan- +6h2xh2y(hxx+hyy)+8hxy(hxh3y+h3xhy)], tity dependsonlyonthepropertiesofthe steady-state (45) A solutions. Thus, determining whether coarsening occurs or not can be decided on the inspection of steady-state where L˜[h] = ( h) 4h is the linear part. We −∇· ∇ −∇ solutions only. This result generalizes our previous one- focus here on the square and hexagonal symmetries for dimensionalstudytotwodimensions[4],wherewefound which the determination of the steady-state solution h0 that coarsening occurs if h2 (which is nothing but the is relatively simple. h 0i amplitude of the pattern) is a decreasing function of q. For a square symmetry, adopting the wave vector In two dimensions [19] we had previously found for the directions which are specified in Appendix E, we can time-dependent Ginzburg-Landau equation and for the write h(x,t) = n,man,m(t)eiq(nx+my). The constant Cahn-Hilliard equation that a certain quantity, different term (n = 0,m = 0) is zero because of the condition P 8 h =0. Giventherealcharacterofh(x,t)anditsparity- Consequently, phase equation eigenvalues are positive, h i symmetrywithrespecttobothspacevariables,itfollows implying instability with respect to phase fluctuations, thata =a anda =a =a . Invari- i.e. coarsening (this also holds for c > 0 and c < 0). ∗n,m n, m n,m n,m m, m 2 4 anceunderπ/−2ro−tationsim−pliesthefollowing−additional Instead, in the case where c >0, we can easily see that 2 relations: eigenvalues are negative, meaning no coarsening at all. Finally,inthecasec <0andc >0,wefindinterrupted 2 4 h(x,y)=h( y,x)=h( x, y)=h(y, x), (46) coarsening: the length scale L of the pattern increases − − − − until reaching a certain maximum wavelength λ=2π/q. whichtranslateinto the followingconditions for the har- Tofixidea,andwithoutlossofgenerality,wesetc = 1 2 monic amplitudes: and deduce the two stationary solutions from Eqs. (52−): an,m =am, n =a n, m =a m,n. (47) 10q4 (10q4)2 224c q6ω − − − − x (q,c )= ± − 4 1 a2, (55) Therefore, at first order (n,m= 1): ± 4 p 224c4q6 ≡ 1 ± which both coincide at the maximum reachable length h(x,t)=a (t)(eiq(x+y)+eiq(x y)+c.c.) , (48) 1 − scale, where coarsening is interrupted. For the hexagonal case, we proceed in the same way, with a = a = a = a a , so the same dynam1ic,1s equa−t1i,o1n hol−d1s,−fo1r the1,f−o1ur≡ha1rmonics. This setting a new truncated Fourier series: series can be written in the more convenient form h0 = h(x,t)=a (t)(eiq/2(x+√3y)+eiqx+eiq/2(x √3y)+c.c.), 4a cosϕ cosϕ and diffusion coefficients (36) and (37) 1 − 1 1 2 (56) canbe explicitly calculatedin the smallamplitude limit: with a = a = a = a = a = a a , 1,1 1,0 0, 1 1, 1 1,0 0,1 1 4 d d so the same dynamics−equati−on−holds f−or any of th≡e six D11 ≃ h2 dq(q3a21)+q3dqa21+q2a21 harmonics [21]. Again, reporting Eq.(56) into Eq.(45), h ϕi(cid:20) (cid:21) we obtain the amplitude equation: 8q3 da2 MAX 1, (49) ≃ hh2ϕi dq a˙1 =a1(ω1+c29q4a21+94c4q6a41), (57) 16q3 da2 D +D +D MAX 1, (50) where ω1 := ω(k = q). Because this equation has the 11 22 12 ≃ h2 dq same structure and the same signs in front of each term h ϕi asinEq. (51),thesameconclusionsasabovearereached, with the substitution q = qMAX = 1/√2 in the last pas- namely we have coarsening, no coarsening and inter- sages. Itisnowevidenthowthesignoftheeigenvaluesis rupted coarsening scenarios depending on the signs of directly related to the increasing or decreasingcharacter the coefficients c and c . 2 4 ofthesteadyamplitudea1 asfunctionofthesteadywave The above results can be extended to other symme- lengthλ=2π/q. Reportingexpansion(48)intoEq.(45) tries,including the 3-foldsymmetry, whichdoes not cor- we find the equation obeyed by a1: respond to a Bravais lattice. In that case, the start- ing Fourier series corresponds to a linear combination of a˙1 =a1(ω1+c220q4a21+224c4q6a41), (51) h(x,t)andh(x x0,t),whereh(x,t)isgivenbyEq.(56) and x =(4π/q−,0) [22]. where ω := ω(k = √2q), and then solve for stationary 0 3 1 It is worthnoting that in the limit of weak amplitude, solutions: the coarsening criterion always corresponds to the requirement that the amplitude of the stationary solution a =0, ω +20c q4a2+224c q6a4 =0. (52) 1 1 2 1 4 1 isadecreasingfunctionofthe wavevector. Thisisatriv- ial consequence of the single harmonic approximation, We conclude that the number ofstationary solutions de- where h (x)=a (q) exponential factors. pendsonvaluesofc2 andc4: coarseningoccurrenceisdi- 0 1 × rectly associated with specific forms of j currents, whose expression determines one of three possible scenarios. VI. UNIVERSALITY CLASSES Let’s first consider the case in which c =0. We find: 4 ω In this section our aim is to extract analytically the 1 a = , (53) 1 20c q4 coarsening exponents β and n, defined as in Sec. I. In r− 2 ordertodeterminetheexponentwemakeuseofthetem- sothatstationarysolutionscorrespondingtothebandof poral behavior of the phase, see Eq. (28), whose ampli- wavevectors such that ω >0 exist only if c <0. Using tude increases as ψ(0)eΩT. The relevant time scale of 1 2 Eqs. (49), we now obtain: the phase instability is therefore set by Ω−1. According to Sec. IV, the unstable mode has an eigenvalue of the 8 1 form Ω = K2 , where is a suitable combination of D11 = 5 h2 c , (54) diffusion co−efficDients andDwhose negative sign indicates h ϕi 2 9 instability. Therefore, if L is the typical size of mounds current j producing n β 1 after a time t, we have T =ε2t, K and ≈ εL 1 faceting 0 L2 3 (q) , (58) |D |≈ t 1 1 increasing slope with q =2π/L[23]. 4 2(1+α) It turns out that the coarsening exponent n only depends on one single property (see below) of the current TABLEII.Thevalueofthecoarseningexponentnandofthe j entering the general equation (3), while its symmetry, exponent β (|m(t)| ≈ tβ for large t, for the two universality classes resulting from our study. The exponent α appears in as well as the pattern symmetry, is definitely irrelevant. therelation j(m)≈1/|m|α for large m. Theonlyessentialingredientiswhetherthecurrentleads ornottoaslopeselection. Wefindn=1/3forjcurrents giving rise to mounds that grow with a certain constant direction but remains constant along the perpendicu- slopeandn=1/4otherwise. Letus showmoreprecisely lar one, so that Eq.(62) can be mapped onto a one- these results. dimensional equation: m +(1/r)m +1/mα, neglecting ′′ ′ Let us consider the square symmetry. We have seen the angular dependence for m. Plugging in it a solution that coarsening occurs if at least one of the expressions of the form m Arγ, we find γ =2/(α+1) and finally, given by (36) or (37) is negative. Consider one scalar ≃ calculating again the scalar products in Eqs. (36) and product entering expression (36): (37) for the square symmetry case, we get h2 = 1 λdx λdy 1 ∂2h0 2. (59) L t1/4 (63) h 11i (2π)2 q2 ∂x2 ∼ Z0 Z0 (cid:18) (cid:19) for any value of α. For systems exhibiting slope selection, the current j(m) Wecangainfurtherinsightfromdimensionalconsider- haszerosforfinitevaluesm oftheslope,thereforem = ations: if m L2/α+1, then m tβ =t2n/α+1, therefore ∗ ∗ ≈ ≈ (∂ h ,∂ h ) is constant everywhere but along domain givingβ =1/(2(α+1)). Wehavefoundthatalltheother x 0 y 0 walls,thathaveafinitebutsmallthickness. Let’sdenote symmetries hereby mentioned produce exactly the same the thickness by δ. In the large wavelength limit δ λ exponents, n and β, pointing to the existence of univer- ≪ we can also assume that inside domain wall there is a sality classes. Our results are summarized in Table II. linearspace dependence forthe slope m: for example,in It is worth comparing our results for n with the cor- Eq. (59), ∂h /∂x = m Ax+By, with A and B real responding values for the one-dimensional growth mod- 0 x ≈ constants, whose exact values are unimportant for our els. The models with constantly increasing slope yield purposes. Estimation of Eq. (59), therefore, yields n=1/4[4,26], aswehavealsofoundintwodimensions. The1dmodelwithfaceting,instead,isknowntoproduce 1 1 λ λ a logarithmic coarsening in the absence of noise [4, 27] h2 dy dx A2 c λ3+o(λ3) h 11i≃ (2π)2q2 ≃ 11 and n = 1/3 when noise is present [28]. We conclude, Z0 Zλ−δ (60) with the caveat of noise and in analogy with the case with c a positive constant. Similar considerations of phase separation processes (see Sec.VII.B for a thor- 11 lead to h2 c λ3 and h2 c λ2, where con- oughdiscussion),thatthedimensionofphysicalspace,d, h 12i ≃ 12 h 1,2i ≃ ϕ stants, again, are positive. We straightforwardly obtain seems[29]tobeirrelevantforourclassofgrowthmodels, D = 2qc /c fromEq.(36)and(D +D +D )= Eq. (3). 11 12 ϕ 11 22 12 − 2q(c c )/c , with c c [24], from Eq. (37): both 12 11 ϕ 11 12 − ≥ coefficients are evidently negative[25], and they have the same q-dependence. Using (58) the coarsening exponent VII. DISCUSSION can be easily extracted: A. The context L t1/3. (61) ∼ Thefieldofcrystalgrowthprocessesbyavapourphase Formodelswithoutslopeselection,thecurrentjhasno hasbeenveryactiveinthepasttwenty-fiveyears,involv- zeros. A prototype of this kind of currents is asymptot- ing experiments, simulations, and analytic studies. Gen- ically represented by j(m) 1/mα, α > 1. Exploiting ≃ | | eral references of special interest for the present article Eq. (14) and coupling it with the asympotic expression are three review papers [6, 30, 31] and one book [32]. of the current, we obtain Rigorous results for model equations included in the 2m 1/mα. (62) class (3) studied here are very rare. We are aware of −∇ ≃ | | two exact inequalities, concerning isotropic currents and We now switch to polar coordinates and make the as- which are in agreement with our universality classes: sumption that mound profile changes only along one Kohn&Yan[33]studiedthefacetedcase,findingn 1/3; ≤ 10 Bo&Li [34] studied the increasing slope case, finding giving a few more details on point (i), here above. n 1/4. Another worth mentioning paper is the study Equation(3)is aclass ofgeneralmodels because j(m) ≤ by Watson&Norris [35] of the faceted case with a three- has the only requirement to be linear at small slopes. fold symmetric current: authors find n = 1/3, also in Suchequationcannotcoveranypossiblemodelofgrowth agreement with our findings. byvapordeposition. Forexample,somestudieshavesug- Less rigorous results are often based on approximate gested a higher order linear term, which would produce evaluations of the different terms appearing in the equa- slowercoarsening[36] (iftwo linearterms ofdifferentor- tionandonhowsuchterms dependonscalelengthL. A derarepresent,crossovereffectsareexpected,depending significant amount of effort has been devoted to such an on their relative strength). A second remark is related approach by Golubović and collaborators [31]. While in to the up-down symmetry of the emerging morphology. the “no faceted” case [36] results are in agreement with Eq. (3) is invariant under the transformation h h → − ours, the faceted case is controversial. More precisely, in but such a symmetry is weak or absent in simulations Ref.[37]thesix-foldsymmetriccasegivesn=1/3,while and experiments. Therefore, symmetry breaking terms in Ref. [38] the four-fold case may give either n = 1/3 have been introduced, for example in Ref. [40]. While in or n = 1/4, depending on the details of the current. d=1suchatermseemstobeirrelevant[41],thequestion Siegert too, already a few years before, had claimed to of its relevance in d = 2 is still open. A final comment find a slower coarsening, n = 1/4, when integrating nu- concerns the form of j(m). Whatever its symmetry, we merically a (3)-like equation with square symmetry and have assumed it is analytical at m = 0. However, some faceted morphology [39]. results [42] suggest that step-edge diffusion might pro- Thepeculiarresultn=1/4forsquaresymmetrywould duce a current which is singular at vanishing slope. be related, according to the above mentioned authors, to the existence of two different types of domain walls, pyramid edges and roof tops: in the former case only one component of the slope changes, while in the latter B. Crystal growth vs phase separation case both components change. Roof tops, which are not present in a regular,periodic square lattice, act as dislo- In the previous Section we have discussed how our re- cationsandwouldplay amajorroleintheir simulations, sultscomparewithotherstudiesoncrystalgrowth. Here slowing down coarsening. However, they also claim that weratherfocusondifferencesandsimilaritieswithphase the system would be frozen in the absence of roof tops, separation processes. because the square pattern would be metastable. The last statement is in contrast to our findings, according First of all, the irrelevance of space dimensionality we to which the square pattern is expected to be linearly have found to be valid for our growth model Eq.(3) is unstable. We think that the square case would require also a well-known feature of phase-separation processes, some more analysis to gain more insights on the role of as long as T > 0. In such context, for a scalar order c defects, as well as on the roles of the specific form of the parameter the coarsening exponents are n = 1/2 for a current, of the initial conditions and of the simulation nonconserved order parameter and n = 1/3 for a con- time. For the sake of completeness we also make our served one [2, 4]. However, in d = 1, such figures are readerawarethatthesameconsiderationsjustexpressed found when noise is present [28], otherwise a slow log- forthe squarecasemightbe extendedto the lessstudied arithmic coarsening appears [27]. The one dimensional rectangular case, since a different coarsening exponent growth model with faceting is equivalent to a conserved canbe foundintheliterature(n=2/7)[31]. Again,also phase separationprocess(so calledmodel B of dynamics for this symmetry, coarsening is bound to the presence or Cahn- Hilliard equation). Therefore, it is not sur- ofdislocationsandthe specific formofthe currentmight prisingthatitgiveslogarithmiccoarseningwithoutnoise be relevant (e.g., its derivability from a potential). andn=1/3inthe presenceofnoise. However,the anal- While comparing with previous studies of the same ogy with phase separation processes cannot be pushed classofequationsisstraightforward,comparingwithKi- further, because our growth model has, in 2d, peculiar netic Monte Carlo simulations and with experiments be- proprieties. comes difficult and dangerous, mainly for two reasons: Even if our multiscale approach is applicable, in prin- (i) Is the system (the simulated system or the real sys- ciple, to any nonequilibrium current j(m), we have con- tem) described by a (3)-like equation or different equa- sidered here the case of symmetric Jacobian matrix , J tions are more appropriate? (ii) Does the simulation which means ∂j /∂m =∂j /∂m . This condition is sat- i ℓ ℓ i or the experiment attain large enough times to probe isfiedbyallcrystalgrowthequationsweareawareofand the asymptotic regime? Honestly, the variety of results itisequivalenttosayingthatdynamicscanbecastintoa of simulations/experiments and the above two questions variationalformulation. Infact,Eq.(3)canbe rewritten prevent from giving a clear picture of such results and as: from connecting them to specific model equations. As forexperiments,wereferthereadertoTable2orRef.[6] ∂h δ and to Table 4.2 of Ref. [32]. We close this Section by = F, (64) ∂t ∇· δm