Stability of a strongly anisotropic thin epitaxial film in a wetting interaction with elastic substrate Mikhail Khenner,1 Wondimu T. Tekalign,2 and Margo S. Levine3 1Department of Mathematics and Computer Science, Western Kentucky University, Bowling Green, KY 42101 2School of Mathematical Sciences, Rochester Institute of Technology, Rochester, NY 14623 3Department of Radiology at Massachusetts General Hospital and Harvard Medical School, Boston, MA 02114 (Dated: January 4, 2011) 1 The linear dispersion relation for surface perturbations, as derived by Levine et al., Phys. Rev. 1 B 75, 205312 (2007) is extended to include a smooth surface energy anisotropy function with a 0 2 variableanisotropystrength(fromweaktostrong,suchthatsharpcornersandslightlycurvedfacets occur on the corresponding Wulff shape). Through detailed parametric studies it is shown that a n combination of a wetting interaction and strong anisotropy, and even a wetting interaction alone a results in complicated linear stability characteristics of strained and unstrained solid films. J PACS: 68.55.J, Morphology of films; 81.15.Aa, Theory and models of film growth; 81.16.Dn, Self- 3 assembly. ] l l I. INTRODUCTION former may be present even when the latter is absent. a h For wetting (non-wetting) films, the solution of the elas- - ticfree-boundaryproblemwithboundaryconditionsthat s Studiesofthemorphologicalinstabilitiesofathinsolid include wetting stress terms, results in additional desta- e filmsareafirststeptowardsunderstandingcomplexphe- m bilizing (stabilizing) contributions in the dispersionrela- nomena such as the formation of a three-dimensional tion. Some stability characteristics have been analyzed . nanoscale islands in strained alloy heteroepitaxy. Such t inRef. [9]withinthe frameworkoflongwaveapproxima- a studies became common after the pioneering works of m tion, where in addition the surface energy is assumed Asaro and Tiller [1], Grinfeld [2] and Srolovitz [3]. The isotropic. This communication extends that work by - classicalAsaro-Tiller-Grinfeldinstabilityisoneofanuni- d adding strong anisotropy and considering not only wet- axiallystressedsolidfilmonarigidinfinitesubstrate. Its n ting films, but also non-wetting films. Here we recognize o variants for the single-component and alloy films on the that the film thickness and the wetting length are two c rigid as well as on the deformable substrates have been independent characteristiclengths,i.e. theformerlength [ studiedandthisresearchcontinues. Reviewsofworkson is determined by film deposition, while the latter one is the single-component films have been published, see for 5 determined by the molecular structure and properties of v instance Ref. [4]. the film-substrate interface. Since wetting length may 6 Film-substrate wetting interaction is a relatively new be, and normally is, less than the deposited film thick- 0 concept in the field of research on morphological insta- ness,theperturbationwavelengthsmaybecomparableto 8 bility and evolution. When surface slopes are not very the film thickness but still much larger than the wetting 5 . large,this additionalwetting energycan be considereda length. In this case the long wavelength approximation 3 function of the local film thickness h, but not the slopes may hold with respect to the wetting length, but not 0 of h [5–7]. In Refs. [8–10, 12] and others it has been with respects to the film thickness. In Sec. IV we show 0 1 shown that wetting interaction damps long-wave insta- that this approach reveals linear stability features that, : bility modes in a certain range of film thickness, thus we believe, went unnoticed in prior publications. v changing the instability spectrum from long-wave type i X to short-wave type. The latter mode of instability is r more relevant to the process of formation of island ar- a rays [13]. In Ref. [9] it is recognized that in the pres- II. PROBLEM STATEMENT enceofwettinginteraction,theboundaryconditionsthat describe stress balance at the film free surface and at Following Refs. [9, 17], we consider a dislocation-free, the film-substrate interface must be augmented by wet- one-dimensional, single-crystal, epitaxially strained thin ting stress terms - that are proportional to the rate of solid film in a wetting interaction with a solid, semi- change of the surface energy with h. (Wetting stress infinite elastic substrate. The film surface z = h(x,t) is called conjoining (or disjoining) pressure when study- evolves due to surface diffusion. This evolution is de- ing dynamics of thin liquid films on substrates [14, 15]. scribed by This pressure is partially responsible for so-called spin- odal instability, which typically leads to film dewetting (rupture); for discussions of spinodal instability, see for 2 −1/2 ∂h ΩDN ∂ ∂h ∂M instanceRef. [16]andreferencestherein.) Wettingstress = 1+ , (1) andlattice-mismatchstresshavedifferentorigin,andthe ∂t kT ∂x (cid:18)∂x(cid:19) ! ∂x 2 where Ω is the atomic volume, D is the adatoms diffu- k is the perturbation wavenumber and h is the uniform 0 sivity, N is the adatomssurface density, kT is the Boltz- thickness of unperturbed planar film): mann factor, and M is the surface chemical potential ω(h ,k,µ,ǫ,ǫ )=Aǫ2(µ+A h k)k3− [18]. Thelatterhascontributionsfromtheelasticenergy 0 γ 1 0 in the film, the anisotropicsurface energy,anda wetting Bǫ[µ(h −1)+h (B h −A )k]k3aexp(−h )+ 0 0 1 0 1 0 interaction: F (Λ−(G+Λ)exp(−h ))k4−∆k6− 0 ∂2γ κ3 ∂2κ (cid:2)ak2{exp(−h0)−B2aexp(−2h0)} .(5) M=E(h)+Ω γ+ κ−δ + + ∂θ2 2 ∂S2 (cid:20)(cid:18) (cid:19) (cid:18) (cid:19) ℓ has been chosen as the length scale and ℓ2/D as(cid:3)the ∂γ ∂2γ time scale. Here ǫ is the misfit strain in the film, and −sign(h ) cosθ ,(2) x ∂h ∂h∂θ µ = µ /µ is the ratio of the film shear modulus to the (cid:18) (cid:19) (cid:21) f s substrate shear modulus. Other parameters are: where θ is the angle that the unit surface normal makes with the reference crystalline direction, say [01] (cho- 8NΩ2α (1+ν )2µ 4NΩ2γ α (1+ν )ν s f f 0 s f f sen along the z-axis, which is normal to the substrate), A= kTℓα2 , B = kTℓ2α2 , γ(h,θ) is the height- and orientation-dependent surface f f (6) energy, κ is the curvature of the surface, S is the ar- clength along the surface and hx is surface slope (note, 2C1 νfC2−αf βfγ0 ∂/∂S = (cosθ)∂/∂x = 1+h2 −1/2∂/∂x). The term A1 = α α , B1 = 2α α ν , B2 = α µ ℓ, (7) x f s f s f f f proportionaltothesmallpositiveparameterδ isthereg- (cid:0) (cid:1) ularizationthatisrequiredinviewofill-posednessofEq. (1)forstronganisotropy,thatiswhenǫγ ≥1/(m2−1)in C1 =αf+αfβsµ−α2sµ2, C2 =4αf+3αfβsµ−4α2sµ2, Eq. (4) below [19–22]. Note also that the mixed deriva- (8) tiveterminEq. (2)isnonlinearandthusithasnoimpact α =2 1−ν , β =1−2ν , (9) on linear stability. f(s) f(s) f(s) f(s) Inthiscommunicationweconsiderthetwo-layerexpo- (cid:0) (cid:1) nential model for the surface energy [5, 7]: NΩ2γ δ 0 F = , ∆= , (10) kTℓ2 γ ℓ2 γ(h,θ)=γ (θ)+(γ −γ (θ))exp(−h/ℓ), (3) 0 t s t whereγ isthesurfaceenergyofthesubstratewhenthere Λ= m2−1 ǫ −1, G=γ /γ , a=G−1−ǫ . (11) s γ s 0 γ is no film, ℓ is characteristic wetting length, and γ (θ) is t the anisotropic surface energy of a thick film: In Eq(cid:0)s. (6)-((cid:1)9) ν is Poisson’s ratio. Note coupling of wetting interaction and misfit strain through the term γ (θ)=γ (1+ǫ cosmθ), ǫ ≥0. (4) proportionaltoǫ(the secondlineofEq. (5)). Thisterm, t 0 γ γ responsible for breaking symmetry between compressive Here γ is the mean surface energy, ǫ is the strength and tensile stress states, drops out of the growth rate 0 γ of anisotropy and m is the integer parameter specifying in the absence of wetting interactions, h0 → ∞. In the anisotropy type (i.e., four-fold, six-fold, etc.) By com- square brackets of this term, −µ and −A1h0k are the parison with experimental and ab initio computational contributions from the wetting stress; see also Refs. [9, studies the two-layer exponential model has been shown 25]. Another contribution from the wetting stress is the the most accurate to-date [23, 24]. In the absence of term proportional to exp(−2h0) in the square brackets anisotropy, γ = γ = const., and δ vanishes. δ is taken in the last line of the equation. t 0 zero also at weak anisotropy, ǫ < 1/(m2−1). In real- Ourgoalistoelucidatetherolesofanisotropy,wetting γ ity, the maximum of γ (θ) might occur at θ = β, where interaction and wetting stress and to characterize film t the non-zero angle β is a misorientation from the refer- stability in the space of dimensionless parameters h0, k, ence direction. Without significant loss of generality we µ, ǫ and ǫγ. Other material parameters will be fixed to assume β =0 in Eq. (4). theirmostcharacteristicvalues. We choosethefollowing The expression for the elastic energy E(h) in Eq. (2) values: D = 1.5×10−6 cm2/s, N = 1015 cm−2, Ω = is derived in Ref. [17] without accounting for wetting 2 × 10−23 cm3, kT = 1.12× 10−13 erg, γ0 = 2 × 103 interaction. Wetting interaction is considered in several erg/cm2, νf = 0.198, νs = 0.217, µf = 1012 erg/cm3, papers,includingRefs. [8–12,23,25,27](inRef. [27]the δ = 5×10−12 erg, and ℓ = 3×10−8 cm. The value of wetting effect arises not from the dependence of surface the characteristic wetting length is of the order of 1 ML energy on thickness, but from the thickness-dependent thickness for InAs or Ge film [23]. We assume strong elastic energy, which cannot be calculated from linear anisotropy, i.e. ǫ > 1/ m2−1 ≡ ǫ(c) and thus Λ > γ γ elasticitytheory). Tothis end, bycombiningexpressions 0. For strained films considered in Sec. IV, we choose (cid:0) (cid:1) derived by us in Refs. [9, 25] we state the dimensionless m = 32 as the most characteristic value [26]. However, lineargrowthratein the longwave limit, kh ≪1 (where as far as the effect of anisotropy on linear stability is of 0 3 interest, similar results are obtained for other common 0.20 86 4 2 3 579 values such as m=4 or m=6. That is, choosing larger m simply means that smaller values of ǫ are above the 0.18 γ critical value ǫ(c). Wetting films require a >0 [25], thus γ we choose γ = 2γ , ǫ(c) < ǫ < 1. For analysis of non- 0.16 s 0 γ γ wetting films (a<0), we choose γ =γ /2, ǫ >ǫ(c). It s 0 γ γ 0.14 isclearthatwettingstressterms(pointedoutabove)are Γ Ε destabilizing (stabilizing)in wetting (non-wetting)films. 0.12 0.10 III. FILMS WITH WETTING INTERACTION AND ZERO MISFIT STRAIN AND WETTING 0.08 STRESS When misfit strain and wetting stress are not present, 0.0 0.2 0.4 0.6 0.8 Eq. (5) reduces to: k ω(h0,k,ǫγ)=F (Λ−(G+Λ)exp(−h0))k4 FIG.1. Contourplot ofthecritical thicknessh(c1) forstrong 0 (cid:2) −∆k6−ak2exp(−h0) . (1) anisotropy, ǫγ >ǫ(γc) =1/15. First, we consider wetting films. It follows th(cid:3)at the perturbations with the wavenumbers larger than k = c destabilized. However note that models of Refs. [23, 27] Λ/∆cannotdestabilize a film of any thickness. (Here, do not allow staightforward separation of the effects of k isnotthecustomarycut-offwavenumber,butisdeter- pc surfaceenergyandmismatchstress,andthus our results minedfromtheconditionω <0foranyh .) However,in 0 can’t be easily compared to these papers. the opposite case k <k only the films of thickness that c Next, we consider non-wetting films. One example is less than the critical, h(c1), are stable: 0 of such material system may be the energetically-driven dewetting of silicon-on-insulator [28, 29]. Repeating the Λk2−∆k4 h <h(c1) =−ln . (2) analysisandreferringto the criticalvalues shownabove, 0 0 a+(G+Λ)k2 it follows that film of any thickness is stable with re- spect to perturbations with wavenumbers larger than With∆=25/9correspondingtothematerialparameters stated above, m=4, and ǫγ =0.1, we obtain kc =0.42. max kc,kc(u) ,wherekc(u) = −a/(G+Λ). Ifkc <k < Taking typical k = 0.1kc in Eq. (2) gives h0(c1) = 6.94, kc(u),(cid:16)then film(cid:17)is stable if h0 >ph0(c1) and unstable other- whichtranslatestothedimensionalvalueof7ML.Fig. 1 wise. Ifk(u) <k <k ,thenthefilmisstableifh <h(c1) showsthecontourplotofh(c1)(k,ǫ ). Itcanbeseenthat c c 0 0 0 γ and unstable otherwise. Finally, if k < min k ,k(u) , stronger anisotropy decreases h(c1). We notice also that c c strong anisotropy destabilizes (0that is, the contribution then the film of any thickness is unstable. Wit(cid:16)h G=0(cid:17).5 proportionalto k4 inthe squarebracketis positive)only and ǫγ =0.1, kc(u) =0.77>kc, and therefore the second relatively thick films, such that possibility, k(u) < k < k , must be dismissed. Typically, c c the firstscenario(k <k <k(u))holds,andthus thereis Λ c c h >h(c2) =−ln . (3) a critical thickness below which the film is unstable [10]. 0 0 G+Λ Results similar to shown above for wetting and For the chosen values, h(c2) = 1.6 ML. Such threshold- non-wetting films can be obtained (numerically) with 0 non-zero wetting stress, since the negative exponent type influence of strong anisotropy is distinctly different exp(−2h )decaysfastcomparedtothe termsin Eq. (1) from the simplified model in which wetting interaction 0 that are proportional to exp(−h ). is absent. The latter model can be obtained by taking 0 h → ∞ in Eq. (1), and thus this equation becomes 0 ω(k,ǫ ) = F Λk4−∆k6 , from which it is clear that γ strong anisotropy has destabilizing influence on a film IV. WETTING FILMS WITH NON-ZERO (cid:0) (cid:1) MISFIT STRAIN AND WETTING STRESS of arbitrary thickness. These findings to some extent echo Refs. [23, 27], where the existence of the critical perturbationamplitude that is necessaryto destabilize a The situation presented in this section is common for film in the presence of a cusp in the surface energy γ(θ) Stranski-Krastanovgrowth of epitaxial thin films. (which is the case below the roughening temperature), As we pointed out in Sec. II, in the presence of mis- has been demonstrated. Thus if a film is thin, critical fit strain and wetting interaction, Eq. (5) contains the amplitude may be unattainable and the film will not be term that is proportional to the first power of misfit 4 strain. Whether this term is destabilizing or stabilizing lized by short-wavelength perturbations, u>0.02, above (for, say, ǫ > 0) depends on the sign of the expression some critical value of the misfit parameter ǫ. Increasing f(k,µ,h )=µ(h −1)+h (B h −A )k. Onlyforsuf- film thickness in the isotropic case to values as large as 0 0 0 1 0 1 ficiently small k and large µ this is positive. Then the 50ℓ only makes the domain of stability to shrink. How- second term in Eq. (5) is stabilizing as shown in Fig. 2. ever, for larger film thickness in the strongly anisotropic Increasing h makes the domain of stabilization smaller. case (Fig. 3(c,d)) two stability domains emerge sepa- 0 As µ is in the range 0.5 - 1.0 for a typical heteroepitax- rated by the domain of instability. The splitting of a ial semiconductor system, the coupling of misfit strain single domain into two domains occurs at r = 0.5. The and wetting interaction has stabilizing effect on an ul- size of stability domains decreases with increasing film trathin film of thickness of the order of several wetting thickness. Overall, the film is less stable with increasing lengths, for longwave perturbations. Note that the stan- anisotropy (as expected). Note that instability in Fig. dard ǫ2-term is always destabilizing for all perturbation 3(c,d)ispresentforsomeuevenwhenmisfitiszero. Re- wavelengths [17]. sponsible for this is the combined destabilizing effect of anisotropy and wetting stress, which together overweigh 5 the stabilizing effect of the wetting layer; also see Sec. III. Similar behaviour is observed for increasing µ while keepingthicknessfixed. We alsonotice thatonly r=0.1 4 casecanbe(probably)capturedbylongwaveapproxima- (a) tion, as kh = 2πru = 0.38 ∼ 1 for r = 0.1,u = 0.6, and is even larger for other values of r in Fig. 3. In order to 3 Μ (b) 1.0 1.0 2 0.5 U 0.5 U S Ε 0.0 S Ε 0.0 -0.5 U 1 -0.5 U -1.0 0.00.10.20.30.40.50.6 -1.0 0.0 0.1 0.2 0.3 0.4 0.5 0.6 u u 0 (a) (b) 0.00 0.05 0.10 0.15 0.20 k 1.0 1.0 0.5 U 0.5 U FIG. 2. Zero level curve of f(k,µ,h0) = µ(h0−1) + Ε 0.0 S Ε 0.0 S h0(B1h0−A1)k. (a) h0 =2, (b) h0 =3. To the left of each -0.5 U -0.5 U curve, the symmetry-breaking term in the longwave growth -1.0 -1.0 rate(5)(secondline)isstabilizing,totheright-destabilizing 0.00.10.20.30.40.50.6 0.00.10.20.30.40.50.6 (when a,ǫ>0). u u (c) (d) In order to demonstrate some effects of arbitrary re- lation between wetting length, film thickness and the perturbation wavelength, in conjunction with strong FIG. 3. Neutralstability curves. (a): ǫγ =0. r=0.1 (solid), anisotropy,we use next the full dimensional growthrate r = 1 (dash), r = 3 (dot). (b-d): ǫγ = 0.01. (b) r = 0.1, expression involving hyperbolic functions of kh, where (c) r =0.518, (d) r = 1 (dashed), r =3 (solid). Domains of k, h are now the dimensional wavenumber and mean surface stability (instability) are marked byS (U). thickness, respectively. (Available on request from au- thors.) (Eq. (5) emerges upon expansion of this growth characterizethe horizontalspacingbetweentwostability rate in powers of small dimensionless parameter kh, re- domains in Fig. 3(c,d), in Fig. 4 we plot the neutral taining the dominant terms of the expansion (longwave stability curve corresponding to the level ǫ = 0. It can approximation), and non-dimensionalization.) As is Eq. be seen that for all reasonable r this spacing does not 5, the full growth rate is quadratic in ǫ, allowing one to exceed 0.3. For comparison, the case of slightly larger explicitly determine the boundaries of neutral stability, anisotropy is also shown. ω =0, in the r−ǫ oru−ǫ planes. Here r andu (dimen- To summarize, we considered all combinations of sionless) are defined by h=rℓ, k =u(2π/ℓ). wetting interaction (through the exponential two-layer Fig. 3showsneutralstabilitycurves,inu−ǫplane,for model),lattice-mismatchandwettingstrains,andstrong ǫγ =0 and 0.01. (For value m=32 used in this Section, anisotropy. Our results demonstrate complicated linear ǫ(c) = 0.001.) For all three values of a film thickness in stability of ultrathin films (h ∼ 1 : 5 wetting lengths). γ the former (isotropic) case, and for the smallest value in In particular, we show that extremely thin (h ∼ 1 : 2 the latter (strongly anisotropic)case,the film is destabi- wetting lengths), unstressed wetting films are not desta- 5 3.0 they are very thin and at any reasonable mismatch stress level. Our final remark concerns two-dimensional 2.5 surfaces and corresponding surface energy anisotropies γ (h ,h ) of the generic form (4) (where contribution t x y (a) in the y-direction is additive, as is commonly assumed). 2.0 We conjecture that, if the surface orientation (of a thick film) is still one of the high-symmetry crystallographic r1.5 (b) orientations, such as [001] or [111], then the effect of such in-plane anisotropy is nonlinear and thus the lat- 1.0 ter anisotropy will not affect the results. This can be qualitatively understood, for instance, by following the 0.5 analysis leading to Eq. (15) in Ref. [8] while accounting for the nonlinear nature of the mixed derivative term in Eq. (2) and the form of Eq. (3). Due to complexity 0.0 of formulation and derivation, the exact proof is beyond 0.00 0.05 0.10 0.15 0.20 0.25 0.30 the scope of this note. The results are also unchanged if u the surface is two-dimensionalbut in-plane anisotropyis zero. FIG.4. Neutralstability curves(seetext). (a) ǫγ =0.01, (b) ǫγ =0.014. Film is stable below each curve. ACKNOWLEDGMENTS bilized by arbitrarily strong anisotropy. Anisotropic, stressed wetting films are destabilized by any level of M.K. acknowledges the support of WKU Faculty mismatch stress, but only in the narrow range of per- ScholarshipCouncilviagrants10-7016and10-7054,and turbation wavenumbers. Such films can remain sta- thanks BrianJ.Spencer (University atBuffalo)for com- ble withrespecttoshort-wavelengthperturbationswhen ments on the first draft of the manuscript. [1] R.J. Asaro and W.A. Tiller, Metall. Trans. 3, 1789 [16] A. Sharma and R. Khanna, Phys. Rev. Lett. 81, 3463 (1972). (1998). [2] M.A. Grinfeld, Sov. Phys. Dokl. 31, 831 (1987). [17] B.J. Spencer, P.W. Voorhees, and S.H. Davis, J. Appl. [3] D.J. Srolovitz, Acta Metall. 37, 621 (1989). Phys. 73, 4955 (1993). [4] H. Gao and W.D. Nix, Annu. Rev. Mater. Sci. 29, 173 [18] W.W. Mullins, J. Appl. Phys. 28(3), 333 (1957); J. (1999);V.A.ShchukinandD.Bimberg,Rev.Mod.Phys. Appl. Phys. 30, 77 (1959). 71, 1125 (1999). [19] C. Herring, Phys. Rev. 82, 87 (1951). [5] C.-H.ChiuandH.Gao,in: S.P.Bakeretal.(Eds.),Thin [20] S. Angenent and M.E. Gurtin, Arch. Rational Mech. Films: StressesandMechanicalPropertiesV,MRSSym- Anal. 108, 323 (1989). posia Proceedings, vol. 356, Materials Research Society, [21] A.DiCarlo, M.E. Gurtin,andP.Podio-Guidugli, SIAM Pittsburgh, 1995, p. 33. J. Appl. Math. 52, 1111 (1992). [6] Z. Suoand Z. Zhang, Phys. Rev. B 58, 5116 (1998). [22] B.J. Spencer,Phys. Rev. E 69, 011603 (2004). [7] M.Ortiz,E.A.Repetto,andH.Si,J.Mech.Phys.Solids [23] S.P.A. Gill and T. Wang, Surf. Sci. 602, 3560 (2008). 47, 697 (1999). [24] M.J. Beck,A.vandeWalle,andM.Asta,Phys. Rev. B [8] A.A.Golovin,M.S.Levine,T.V.Savina,andS.H.Davis, 70, 205337 (2004). Phys. Rev. B 70, 235342 (2004). [25] M. Khenner, Phys. Rev. B 77, 165414 (2008); Math. [9] M.S. Levine, A.A. Golovin, S.H. Davis, and P.W. Model. Nat. Phenom. 3, 16 (2008). Voorhees, Phys. Rev. B 75, 205312 (2007). [26] J. Tersoff, B.J. Spencer, A. Rastelli, and H. von Kanel, [10] Y.PangandR.Huang, Phys.Rev.B 74,075413(2006). Phys. Rev. Lett. 89(19), (2002). [11] W.T.TekalignandB.J.Spencer,J.Appl.Phys. 96,5505 [27] H.R.EisenbergandD.Kandel,Phys.Rev.Lett. 85,1286 (2004). (2000); Phys. Rev. B 66, 155429 (2002). [12] M.D. Korzec and P.L. Evans, Physica D 239, 465 [28] P. Sutter, W. Ernst, Y.S. Choi, and E. Sutter, Appl. (2010). Phys. Lett. 88, 141924 (2006). [13] C.-h. Chiu, Phys. Rev. B 69, 165413 (2004). [29] D.T. Danielson, D.K. Sparacin, J. Michel, and L.C. [14] V.M.Starov,M.G.Velarde,andC.J.Radke,Wettingand Kimerling, J. Appl. Phys. 100, 083507 (2006). Spreading Dynamics (CRC, Boca Raton, 2007). [30] A.A. Golovin, S.H. Davis, and A.A. Nepomnyashchy, [15] J.Israelachvili, Intermolecular and Surface Forces (Aca- Physica D 122, 202 (1998). demic, London, 1991).