Analysis of a combined influence of substrate wetting and surface electromigration on a thin film stability and dynamical morphologies
Abstract
A PDEbased model combining surface electromigration and wetting is developed for the analysis of morphological stability of ultrathin solid films. Adatom mobility is assumed anisotropic, and two directions of the electric field (parallel and perpendicular to the surface) are discussed and contrasted. Linear stability analyses of smallslope evolution equations are performed, followed by computations of fully nonlinear parametric evolution equations that permit surface overhangs. The results reveal parameter domains of instability for wetting and nonwetting films and variable electric field strength, nonlinear steadystate solutions in certain cases, and interesting coarsening behavior for strongly wetting films.
Submitted to the special issue ”Nanoscale wetting of solids on solids” of the journal Comptes Rendus Physique (Olivier PierreLouis, Univ. Lyon, Editor)
Keywords: Electromigration ; Wetting ; Dynamics of thin solid films
I Introduction
Since the early 1990’s there has been an interest in understanding and accessing the effects of electromigration on kinetic instabilities of crystal steps St ()UCS () and epitaxial islands KD ()KKHV (), and on morphological stability of epitaxial films BMOPL ()TGM1 (). Among other applications (which were developed primarily at the microscale), electromigration has been also used for fabrication of nanometersized gaps in metallic films. Such gaps are suitable for testing of the conductive properties of single molecules and control of their functionalities VFDMSKBM ()GSPBT (). For instance, Ref. GSPBT () describes fabrication of nanoscale contacts by using electromigration to thin down and finally break the epitaxially grown ultrathin (10 ML) Ag films wetting the Si(001) substrate. The gap between contacts can be cyclically opened and closed by applying electromigration current at 80 K to open the gap, and enabling surface diffusion by annealing to the room temperature, to close it.
Thus for this and other emerging applications at the nanoscale WH (); SV () it seems important to understand and characterize the effects of substate wetting and electromigration that are simulaneously active in the physical system. This paper combines these effects in a model that is based on an evolution equation(s) for the continuous profile of the film surface. The focus is on wetting films with isotropic surface energy, but with anisotropic adatom mobility KD (); SK (), although the model allows any combination of wetting properties and anisotropies. We factor in and discuss the effects on film stability and morphological evolution of the electric field that is either parallel, or perpendicular to the initial planar surface of the film, and do not limit considerations to small deviations from planarity, i.e. the arbitrary surface slopes and even surface overhangs are permitted by the model. Models of wetting appropriate for continuumlevel modeling of the surface diffusionbased dynamics of solid films have been developed and discussed extensively primarily in the context of thin film heteroepitaxy CG ()K2 (). Our analysis is based on one such model, called the twolayer exponential model for the surface energy ORS (); GolovinPRB2004 (); LevinePRB2007 (), GillWang ()BWA () (which is particularly useful when the surface energy is anisotropic), but other models of wetting discussed in Refs. CG ()K2 () can be used instead, and the results are expected to be qualitatively similar. The goal of modeling in this paper is not to match the theoretical results to the experiment GSPBT () and thus help in understanding the experiment, but rather to provide the broad analysis of the interplay of two effects (wetting and electromigration) that, at least to our knowledge, has not been addressed in prior publications.
Ii Problem Statement
We consider a 2D singlecrystal film of unperturbed height with the 1D parametric surface
, where
and are the Cartesian coordinates of a point on a surface, is the time and is the parameter along the surface.
The origin of the Cartesian reference frame is on the substrate, and along the substrate (the direction) the film is
assumed infinite. The axis is in the direction normal
to the substrate and to the initial planar film surface.
Surface marker particles will be used for computations of surface dynamics Trygg (). Thus and ()
in fact represent the coordinates of
a marker particle, which are governed by the two coupled parabolic PDEs Sethian1 (); BKL ():
(1)  
(2) 
Here (and below) the subscripts and denote partial differentiation with respect to these variables, is the normal velocity of the surface which incorporates the physics of the problem, and . Here is the surface arclength and the surface orientation angle, i.e. the one that the unit surface normal makes with the reference crystalline direction (chosen along the axis). If the surface slopes are bounded at all times (surface does not overhang), then it is more convenient to describe surface dynamics by a single evolution PDE for the height function of the film. Eqs. (1), (2) can be easily reduced to such “hequation”, which we will use for analysis; however, Eqs. (1), (2) will be used for most computations in this paper. Similar parametric approach was used recently in Ref. OL () for the computation of a hillandvalley structure coarsening in the presence of material deposition (growth) and strongly anisotropic surface energy.
Assuming that temperature is sufficiently high and surface diffusion is operative, the normal velocity is given by
(3) 
where is the adatoms diffusivity, the atomic volume, the adatoms surface density, the Boltzmann’s factor, the surface chemical potential, the anisotropic adatom mobility, the applied electric field, the effective charge of adatoms, , if the electric field is vertical, or , if the electric field is horizontal, and is used to select either stabilizing, or destabilizing action of the field for the chosen combination of the vertical or horizontal orientation of the field and the mobility . The two values of correspond to two possible field orientations once either the horizontal, or the vertical field direction has been chosen (that is, field directed updown, or leftright). The first term in describes the highT surface diffusion, the second term describes surface diffusion enabled by electromigration M ().
The surface chemical potential is assumed to have the contributions from the surface energy and the surface wetting interaction with the substrate/film interface:
(4) 
where is the surface mean curvature and, in the general case, is the height and orientationdependent, i.e. anisotropic, surface energy. (Note again that here stands for the shortest distance between the substrate and a chosen point on a film surface; this distance is the height of the surface if there is no overhangs.) In this paper we focus on the effects due to anisotropic adatom mobility, thus we will use the simpler isotropic model for the surface energy CG ()K2 ():
(5) 
where is the (constant) energy of a substrate/gas (or vacuum) interface, the characteristic wetting length, and the constant energy of a crystal/gas interface (that is, of the film surface). Eq. (5) is the interpolation between the two energies. In the limit of a thick film, , only the latter energy is retained in this expression (because the intermolecular forces between the substrate and the surface molecules are relatively shortranged), and in the limit of a film of zero thickness only the former energy is retained. Despite that Eq. (5) is phenomenological, it matches surprisingly well the experiments and the abinitio calculations (at least for latticemismatched systems) GillWang (); WKSM (); BWA ().
Finally, we assume the (dimensionless) anisotropic adatom mobility in the form SK ()
(6) 
where is the number of symmetry axes and is the angle between a symmetry direction and the average surface orientation. is a parameter determining the strength of the anisotropy. Throughout the paper we present results either for (isotropic case), or for and . For the latter set of parameters values the graphs of the functions and are shown in Fig. 1.
Remark 1. In the limit (where wetting effect is not operative), the present problem is most closely related to the one analyzed by Schimschak and Krug in Ref. SK (). The essential difference is that these authors determine the electric field from the solution of the potential equation with the appropriate boundary conditions on the material boundaries, including the moving surface KAIN (); B (). Thus their solution for the electric field is nonlocal, unlike the local approximation used in this paper. The local nature of the electric field explains why we did not detect traveling wave solutions, which are the hallmark of Refs. SK (); B ().
Next, we choose as the length scale and as the time scale and write the dimensionless counterparts of Eqs. (1), (2), where we use same notations for dimensionless variables:
(7)  
(8)  
(9) 
Notice that the dimensionless forms of and the parametric expression for coincide with their dimensional forms, and that Eq. (5) has been substituted in Eq. (4). For conciseness, we keep differentiation with respect to the arclength, , in the dimensionless equations, but their most transparent forms for the computations result when the differentiation with respect to is replaced with the differentiation with respect to the parameter , using . In Eqs. (7)  (9) is the surface diffusion parameter, is the strength of the electric field, and is the ratio of substrate to film surface energies. For wetting films , for nonwetting films . Notice that may take on positive or negative values through the parameter .
As was mentioned above, if there is no overhangs then Eqs. (7)(9) can be reduced to the single dimensionless evolution equation for the surface height :
(10) 
where
(11) 
, if the electric field is vertical, or , if the electric field is horizontal. Also is given by Eq. (6), where is replaced by .
Using values: cm/s, cm, erg/cm, cm, erg, cm (0.3 nm, or 1ML), gives . Using (where statcoulombs is the absolute value of the electron charge) and the minimum applied voltage difference Volts acting across the typical distance of the order of the film height, nm, gives erg/cm, which translates to . Applied voltage in surface electromigration experiments at the nanoscale can be as high as 1 V GSPBT (), thus in this study we explored the range of field strengths .
In the following sections we analyze several representative situations.
Iii Vertical electric field
iii.1 Linear stability analysis
The small slope approximation of Eq. (10) reads:
(12) 
where the mobility has been linearized about the flat surface , that is, . The last two terms are the simplest nonlinearities from the expansion of the electromigration flux that involve and . Without loss of generality we will take here and elsewhere in this paper. (See Fig. 1. When mobility is isotropic, i.e. , then , as is seen from Eq. (6).) Notice that the anisotropy of the mobility, , does not have an effect on linear stability, as it enters in the coefficients of the nonlinear terms. Also note that the last three terms can be written in a conservative form similar to the terms in the first line of the equation and to the first term in the second line (and this is how they are implemented in the code).
Introducing small perturbation in Eq. (12) by replacing with (where is the dimensionless unperturbed film height), linearizing in and assuming normal modes for , gives the perturbation growth rate , where is the perturbation wavenumber:
(13) 
Remark 2. When wetting interaction is absent (thick film: ), Eq. (13) reduces to the standard one, , which reflects the stabilizing action of the surface diffusion and either stabilizing (, electric field is in the positive direction), or destabilizing (, electric field is in the negative direction) action of the electric field. Such film is absolutely stable when , but when it is long waveunstable.
iii.1.1 Analysis of Eq. (13)

Wetting films . From Eq. (13) one notices that wetting films are absolutely linearly stable when , but they are long waveunstable when (see Fig. 2). The shortwavelength cutoff wavenumber, the maximum growth rate, and the wavenumber at which the latter occurs are:
(14) and are plotted in Figs. 3 and 4. The film stability decreases with increasing , and this trend saturates around (3 nm). That is, for the stated field strength films of thickness do not ”feel” the stabilizing presence of the substrate, and they are as stable as the films that do not interact with the substrate at all. Of course, increasing field strength makes the film less stable, but increasing makes it more stable, since the substrate energy provides stabilizing effect. In Figs. 3(c) and 4(c) the rate of decrease of and with (i.e., the rate of stabilization) increases fast with decreasing . For at and the entire dispersion curve is below the axis, see the dashed curves in Figs. 3(c) and 4(c), and thus the film is absolutely linearly stable for and . By analysing the condition for longwave instability, , one can get an understanding of why this happens.
Figure 3: Wetting films and vertical electric field: Instability cutoff wavenumber . (a) vs. , ; (b) vs. , (solid line), (dashed line). ; (c) vs. , (solid line), (dashdotted line), (dashed line). . Figure 4: Wetting films and vertical electric field: Maximum growth rate . Cases (a)(c) are as in Fig. 3. Indeed, this condition is equivalent to , and for the typical values of and stated above and for moderately large the righthand side of the latter inequality is negative, thus the longwave instability occurs for any film height . This is similar to the situation when wetting effects are absent, see Remark 2 (but of course the spectrum of unstable wavenumbers and the magnitude of the maximum growth rate are different). There is a critical thickness below which the film is absolutely linearly stable only if , i.e. when is of the order of at least one hundred.

Nonwetting films . In this case, for the longwave instability occurs for any film thickness and any electric field strength  which is again similar to the situation when wetting effects are absent, see Remark 2. When the longwave instability (with and numerically similar to those shown in Eq. (14)) occurs when the film height is less than critical, . When the film is absolutely linearly stable. is plotted in Fig. 5. It is clear that increases with increasing either , or , or both parameters.
iii.2 Nonlinear surface dynamics of wetting films
Computations with the smallslope Eq. (12) (where ) were performed for and varying field strengths that satisfy the condition for longwave instability, . Computational domain was with periodic boundary conditions, and the initial condition was the smallamplitude cosineshaped perturbation of the flat surface All computations produced steadystate solutions that have the shape of a vertically stretched cosine curve with a fairly large amplitude. However, neither of these steadystates were confirmed when instead fully nonlinear parametric equations (7)  (9) were computed.^{1}^{1}1Of course, we first carefully checked that values of and from the linear stability analyses are reproduced in the computations of the parametric equations. This is one of the methods that we used for checking that the parametric code is errors free. Fig. 6 shows the computed morphology. Overhangs are clearly visible, and the bottom of the ”pit” flattens out as it approaches the substrate due to increase of the repulsive force, promoting overhangs in the vicinity.
Similarly, when Eq. (12) was computed with the random smallamplitude initial condition on the domain (the mobility was again isotropic), the result was a perpetually coarsening hillandvalley structure. This was again not confirmed in the computations of parametric equations. The late time morphology computed using the parametric equations is shown in Fig. 7(a). The surface develops deep “pockets” whose walls overhang and eventually merge, resulting in tear dropshaped voids trapped in the solid; this can be seen, for instance, on the interval . Fig. 7(b) shows the typical random initial condition.
Iv Horizontal electric field
iv.1 Linear stability analysis
In the case of the horizontal electric field the small slope approximation of Eq. (10) reads:
(15) 
Again here we retained two simplest nonlinearities from the expansion of the electromigration flux term. The coefficient in the thirdtolast term has been set equal to one. This nonlinear term and the last (nonlinear) term in Eq. (15) do not have an effect on linear stability. Unless mobility is isotropic, the secondtolast term of Eq. (15) has an effect on linear stability. Assuming anisotropy, the growth rate is
(16) 
Comparing Eqs. (13) and (16), it is clear that one can obtain stability properties in the horizontal field case by replacing in Eq. (13) with . Since (see Fig. 1) for the chosen parameter values in Eq. (6), then for instability one must have (the electric field is in the positive direction). Thus, for instance, the stability analysis of wetting films in the paragraph 1. of Sec. III.1.1 translates directly into the case of the horizontal field simply by replacing in the formulas by , where . Figs. 3(b) and 4(b) are also valid in the case of the horizontal field if the values along the horizontal axes are understood as values of the product (thus the absolute characteristic values are roughly three times smaller than in the case of the vertical field).
iv.2 Nonlinear surface dynamics of wetting films
iv.2.1 Periodic steadystates from sinusoidal perturbations
Computations with either the smallslope Eq. (15), or the parametric equations (7)(9), of the evolution of a one wavelength ()^{2}^{2}2In this computation is not necessarily equal to ., smallamplitude cosine curveshaped perturbation on a periodic domain resulted in steadystate profiles which have a shape of a vertically stretched cosine curve with a fairly large amplitude. The steadystate profiles often displayed a more sharply peaked bottom and less curved walls than the cosine curve, and the amplitude is significantly smaller in the (fully nonlinear) parametric case. In fact, it can be noticed from Fig. 8 that the film is nowhere close to dewetting the substrate for all tested field strengths, notwithstanding that the stabilizing influence of the substrate is minimal for the chosen parameters values (see the discussion in Sections III.1.1, paragraph 1, and IV.1).
Remark 3. When Eq. (15) is used in a computation, the last term there is at large responsible for the saturation of the surface slope and the existence of the steady state. When this nonlinear term is omitted from the equation, the slope increases until the computation breaks down.
In the computations of the parametric equations, unless is larger than approximately , the surface evolves towards the steadystate mostly by vertical stretching of the initial shape. Perturbations of larger wavelength develop a largeamplitude, hillandvalley type distortions which slowly coarsen into a steadystate, cosine curvetype shape. Fig. 8 shows the amplitude of the steadystate profile, , and the height of the profile at it’s lowest point, , vs. .
It can be seen that, at least for moderate strengths of the electric field, growth of the amplitude is logarithmic in the vicinity of , and for growth is linear. The numerically determined slope of the linear section of curve is 0.208, the one of the curve is 0.217, and the one of the curve is 0.209, which suggests that the slope is insensitive (or very weakly sensitive) to the field strength. The decay of with increasing wavelength mirrors the growth of the amplitude, thus the steadystate surface shape is symmetrically vertically stretched about the equilibrium surface position . All computed steadystates are stable with respect to imposition of random small and largeamplitude point perturbations, which we confirmed by computing the dynamics of such perturbed shapes.
iv.2.2 Coarsening of random initial roughness
We employed fully nonlinear simulations based on parametric equations for determination of the coarsening laws at increasing electric field strength and at variable wetting strength (characterized by values of the parameters and ). Computations were performed on the domain with periodic boundary conditions; all runs were terminated after the surface evolved into a largeamplitude hillandvalley structure with 35 hills. Unless wetting is strong, i.e. is small and is large  the case that is discussed in more detail below  the slopes of the hills are constant 24 during coarsening, except for the short initial period. Figs. 9  11 are the loglog plots of the averaged, over ten realizations, maximum surface amplitude and the averaged mean horizontal distance between valleys (kinks) vs. time.
Remark 4. We also attempted to compute coarsening dynamics resulting from the smallslope Eq. (15). While the hillandvalley structure does emerge and coarsens with time, the characteristic constant hill slope is nearly 90. This suggests that additional nonlinear terms must be retained in Eq. (15) for predictive computations and raises the question what terms must be retained. We had not tried to obtain an adequate nonlinear smallslope model, leaving this agenda to future research.
Of course, increasing the field strength results in faster coarsening, as the times needed for coarsening into a “final” structure are , , and for field strengths and , respectively. These values also point out the decrease of the rate of change of a coarsening rate with the increasing field strength. Other than that, the coarsening laws are similar for the three tested field strengths. Fits to the data in the case are shown in Figs. 9(a,b). At small times coarsening is fast (exponential, see Fig. 9(b)), then it changes to a slower power law with the exponent in the range . (Since the values of the amplitude’s logarithm are negative, we were unable to fit the exponential law to the data in Fig. 9(a), thus we fitted the quadratic, which results in the type law.)
Remarkably, when is decreased from 10 to 6 a very different coarsening behaviors emerge, see Figs. 10(a,b) and 11(a,b). At , after the period of slow powerlaw coarsening at intermediate times, the coarsening accelerates to exponential coarsening or type coarsening at late times (Figs. 10(a,b)). (For reference, the linear stability in the case is shown in Figs. 3(c) and 4(c) by the dashdotted line.) At we did not output a sufficient number of data points in the beginning of the computation, but it is expected that initial coarsening is still faster than the power law. The amplitude starts to decrease towards the end of the simulation, and the kinkkink distance coarsens fast on the entire time interval.^{3}^{3}3We also mention that this accelerated coarsening dynamics is not reproduced by the smallslope Eq. 15, as can be expected from the Remark 4. While there is some evidence of acceleration, this computation breaks down too fast.
It seems certain that such unusual dynamics in the strongly wetting states emerges due to nonlinear ”overdamping” of electromigrationinduced faceting instability by the surfacesubstrate interaction force. Toward this end, in Figs. 12(a,b) the surface shapes and surface slopes for and cases are compared at the time when seven hills formed on the surface. In the former case the hills are steep, have rather uniform height and their slopes are almost straight lines. In the latter case they are more irregular, ”rounded”, and the average height is smaller. It is reasonable to expect that in the opposite case of nonwetting films, the surfacesubstrate interaction will instead “sharpen up” the hills, i.e. increase their slopes and make the surface structures appear more spatially and temporally uniform.
Finally, we remark that the discussed coarsening laws for strong wetting cases also are qualitatively different from the laws governing coarsening in the absence of wetting, but in the presence of deposition, attachmentdetachment, strong surface energy anisotropy, and interface kinetics OL (); SGDNV (). Indeed, only the coarsening exponents shown in Fig. 9 (week wetting) are within the same range (0.10.5) for nearly the whole computational time interval, as are the exponents computed in the cited papers.
V Conclusions
In this paper the effects of electromigration and wetting on thin film morphologies are discussed, based on the continuum model of film surface dynamics. It has been shown that wetting effect modifies significantly the stability properties of the film and the coarsening of electromigrationinduced surface roughness. Also it has been shown that the smallslope evolution equations that were employed in many studies of electromigration effects on surfaces, are often inadequate for the description of stronglynonlinear phases of the dynamics. It is expected that the account of the surface energy anisotropy and the electric field nonlocality (through the solution of the moving boundary value problem for the electric potential) will lead to uncovering of even more complicated behaviours.
Acknowledgements
This paper was inspired by the undergraduate research project “An analysis of the electromigration of atoms along a crystal surface using the Method of Lines Approach” by Kurt Woods (MATH 498 Senior Seminar, WKU Department of Mathematics and Computer Science, Fall 2011).
References
 (1) S. Stoyanov, “Currentinduced step bunching at vicinal surfaces during crystal sublimation”, Surf. Sci. 370, 345 (1997).
 (2) D.J. Liu, J.D. Weeks, and D. Kandel, “Currentinduced step bending instability on vicinal surfaces”, Phys. Rev. Lett. 81, 2743 (1998).
 (3) M. Dufay, J.M. Debierre, and T. Frisch, “Electromigrationinduced step meandering on vicinal surfaces: Nonlinear evolution equation”, Phys. Rev. B 75, 045413 (2007).
 (4) J. Chang, O. PierreLouis, and C. Misbah, “Birth and morphological evolution of step bunches under electromigration”, Phys. Rev. Lett. 96, 195901 (2006).
 (5) O. PierreLouis, “Local electromigration model for crystal surfaces”, Phys. Rev. Lett. 96, 135901 (2006).
 (6) J. Quah and D. Margetis, “Electromigration in macroscopic relaxation of stepped surfaces”, Multiscale Model. and Simul. 8, 667 (2010).
 (7) V. Usov, C.O. Coileain, and I.V. Shvets, “Influence of electromigration field on the step bunching process on Si(111)”, Phys. Rev. B 82, 153301 (2010).
 (8) J. Krug and H.T. Dobbs, “Currentinduced faceting of crystal surfaces”, Phys. Rev. Lett. 73, 1947 (1994).
 (9) M. Schimschak and J. Krug, “Surface electromigration as a moving boundary value problem”, Phys. Rev. Lett. 78, 278 (1997).
 (10) F. Hausser, S. Rasche, and A. Voigt, “The influence of electric fields on nanostructuressimulation and control”, Math. Comput. Simul. 80, 14491457 (2010).
 (11) O. PierreLouis and T.L. Einstein, “Electromigration of single layer clusters”, Phys. Rev. B 62, 13697 (2000).
 (12) F. Hausser, P. Kuhn J. Krug and A. Voigt, “Morphological stability of electromigrationdriven vacancy islands”, Phys. Rev. E 75, 046210 (2007).
 (13) P. Kuhn, J. Krug, F. Hausser, and A. Voigt, “Complex shape evolution of electromigrationdriven singlelayer islands”, Phys. Rev. Lett. 94, 166105 (2005).
 (14) F. Barakat, K. Martens, and O. PierreLouis, “Nonlinear wavelength selection in surface faceting under electromigration”, Phys. Rev. Lett. 109, 056101 (2012).
 (15) D. Maroudas, “Surface morphological response of crystalline solids to mechanical stresses and electric fields”, Surf. Sci. Reports 66, 299346 (2011).
 (16) V. Tomar, M.R. Gungor, and D. Maroudas, “Currentinduced stabilization of surface morphology in stressed solids”, Phys. Rev. Lett. 100, 036106 (2008).
 (17) L. Valladares, L.L. Felix, A.B. Dominguez, T. Mitrelias, F. Sfigakis, S.I. Khondaker, C.H.W. Barnes, and Y. Majima, “Controlled electroplating and electromigration in nickel electrodes for nanogap formation”, Nanotechnology 21, 445304 (2010).
 (18) T. Taychatanapat, K.I. Bolotin, F. Kuemmeth, and D.C. Ralph, “Imaging electromigration during the formation of break junctions”, Nano Lett. 7, 652656 (2007).
 (19) G. Gardinowski, J. Schmeidel, H. Phnur, T. Block, and C. Tegenkamp, “Switchable nanometer contacts: Ultrathin Ag nanostructures on Si(100)”, Appl. Phys. Lett. 89, 063120 (2006).
 (20) Z.J. Wu and P.S. Ho, “Size effect on the electron wind force for electromigration at the top metaldielectric interface in nanoscale interconnects”, Appl. Phys. Lett. 101, 101601 (2012).
 (21) D. Solenov and K.A. Velizhanin, “Adsorbate transport on graphene by electromigration”, Phys. Rev. Lett. 109, 095504 (2012).
 (22) C.H. Chiu and H. Gao, in: S.P. Baker et al. (Eds.), Thin Films: Stresses and Mechanical Properties V, MRS Symposia Proceedings, vol. 356, Materials Research Society, Pittsburgh, 1995, p. 33.
 (23) Z. Suo and Z. Zhang, “Epitaxial films stabilized by longrange forces”, Phys. Rev. B 58, 5116 (1998).
 (24) M. Ortiz, E.A. Repetto, and H. Si, “A continuum model of kinetic roughening and coarsening in thin films”, J. Mech. Phys. Solids 47, 697 (1999).
 (25) YaPu Zhao, “Morphological stability of epitaxial thin elastic films by van der Waals force”, Arch. Appl. Mech. 72, 7784 (2002).
 (26) J.N. Aqua, T. Frisch, and A. Verga, “Ordering of strained islands during surface growth”, Phys. Rev. E 81, 021605 (2010).
 (27) T.O. Ogurtani, A. Celik, and E.E. Oren, “Morphological evolution in a strainedheteroepitaxial solid droplet on a rigid substrate: Dynamical simulations”, J. Appl. Phys. 108, 063527 (2010).
 (28) A.A. Golovin, M.S. Levine, T.V. Savina, and S.H. Davis, “Faceting instability in the presence of wetting interactions: A mechanism for the formation of quantum dots ”, Phys. Rev. B 70, 235342 (2004).
 (29) M.S. Levine, A.A. Golovin, S.H. Davis, and P.W. Voorhees, “Selfassembly of quantum dots in a thin epitaxial film wetting an elastic substrate”, Phys. Rev. B 75, 205312 (2007).
 (30) Y. Pang and R. Huang, “Nonlinear effect of stress and wetting on surface evolution of epitaxial thin films”, Phys. Rev. B 74, 075413 (2006).
 (31) B.J. Spencer, “Asymptotic derivation of the gluedwettinglayer model and contactangle condition for StranskiKrastanow islands”, Phys. Rev. B 59, 2011 (1999).
 (32) S.P.A. Gill and T. Wang, “On the existence of a critical perturbation amplitude for the StranskiKrastanov transition”, Surf. Sci. 602, 3560 (2008).
 (33) M. Khenner, W.T. Tekalign, and M. Levine, “Stability of a strongly anisotropic thin epitaxial film in a wetting interaction with elastic substrate”, Eur. Phys. Lett. 93, 26001 (2011).
 (34) M. Khenner, “Comparative study of a solid film dewetting in an attractive substrate potentials with the exponential and the algebraic decay”, Math. Model. Nat. Phenom. 3, 1629 (2008).
 (35) M. Khenner, “Morphologies and kinetics of a dewetting ultrathin solid film”, Phys. Rev. B 77, 245445 (2008).
 (36) L.G. Wang, P. Kratzer, M. Scheffler, and N. Moll, “Formation and stability of selfassembled coherent islands in highly mismatched heteroepitaxy”, Phys. Rev. Lett. 82, 4042 (1999).
 (37) M.J. Beck, A. van de Walle, and M. Asta, “Surface energetics and structure of the Ge wetting layer on Si(100)”, Phys. Rev. B 70, 205337 (2004).
 (38) G. Tryggvason, B. Bunner, A. Esmaeeli, D. Juric, N. AlRawahi, W. Tauber, J. Han, S. Nas, and Y.J. Jan, “A fronttracking method for the computations of multiphase flow”, J. Comput. Phys. 169, 708 (2001).
 (39) J.A. Sethian, “A review of recent numerical algorithms for hypersurfaces moving with curvaturedependent speed”, J. Diff. Geom. 31, 131 (1990).
 (40) R.C. Brower, D.A. Kessler, J. Koplik, and H. Levine, “ Geometrical models of interface evolution”, Phys. Rev. A 29, 13351342 (1984).
 (41) C. Ograin and J. Lowengrub, “Geometric evolution law for modeling strongly anisotropic thinfilm morphology”, Phys. Rev. E 84, 061606 (2011).
 (42) W.W. Mullins, “Solid surface morphologies governed by capillarity”, in Metal Surfaces: Structure, Energetics and Kinetics, 1766 (1963) (American Society for Metals, Cleveland, OH).
 (43) M. Khenner, A. Averbuch, M. Israeli, and M. Nathan, “Numerical simulation of grain boundary grooving by level set method”, J. Comp. Phys. 170, 764 (2001).
 (44) R.M. Bradley, “Electromigrationinduced propagation of nonlinear surface waves”, Phys. Rev. E 65, 036603 (2002).
 (45) T.V. Savina, A.A. Golovin, S.H. Davis, A.A. Nepomnyashchy, and P.W. Voorhees, “Faceting of a growing crystal surface by surface diffusion”, Phys. Rev. E 67, 021606 (2003).