First-principles calculation of the effect of strain on the diffusion of Ge adatoms on Si and Ge (001) surfaces ∗ A. van de Walle, M. Asta, and P. W. Voorhees Materials Science & Engineering Department, Northwestern University, Evanston, IL 60208-3108 3 (Dated: February 2, 2008) 0 First-principles calculations are used to calculate the strain dependencies of the binding and 0 diffusion-activationenergiesforGeadatomsonbothSi(001)andGe(001)surfaces. Ourcalculations 2 revealthatthebindingandactivationenergiesonastrainedGe(001)surfaceincreaseanddecrease, n respectively, by 0.21 eV and 0.12 eV per percent compressive strain. For a growth temperature a of 600 ◦C, these strain-dependencies give rise to a 16-fold increase in adatom density and a 5-fold J decrease in adatom diffusivity in the region of compressive strain surrounding a Ge island with a 3 characteristic size of 10 nm. 1 ] In heteroepitaxial growth of lattice-mismatched thin icsinGe/Si(001). Thissystemrepresentsoneofthemost h films, the Stranski-Krastanov (SK) growth mode has widely studied examples of QD formation induced by c e been widely investigated as a basis for self-assembling SK growth,17,18,19,20,21,22 yet to date the effect of strain m arraysof coherent nanostructuredislands, commonly re- upon Ge adatom binding energies and diffusion rates on ferred to as quantum dots (QDs).1 The need for highly Ge wetting layers remains unstudied. We employ first- - t monodisperseQDarraysinsemiconductoroptoelectronic principlescalculationsto computethe straindependence a t deviceapplicationshasmotivatedextensiveresearchinto ofGeadatomenergeticsonSi(001)andGe(001)surfaces. s the microscopic mechanisms influencing the evolution of Our calculations were set up as follows. We con- . at island size distributions during SK growth. While the sider adatom diffusion in the dilute limit (as opposed m stabilizing thermodynamic effects associated with the to dimer diffusion23). The initial atomic configuration elastic interactions between strained islands have been for all of our calculations is taken to be the well-known - d investigated in some detail (e.g., Refs. 2,3, and refer- minimumenergyc(4×2)reconstruction(depicted inthe n ences cited therein), the role of various proposed kinetic upper left corner of Fig. 1). The adatom calculations o mechanisms on the development of island size distribu- are performed using a supercell geometry where an ar- c tions remains less clear. tificial three-dimensional periodicity is imposed on the [ Within a self-consistent mean-field rate theory, Ko- system. This artificial periodicity is harmless, provided 1 duvely and Zangwill4 demonstrated that with decreas- that convergence with respect to the distance between v ing island-island separation, a strain-mediated decrease the periodicimagesis achieved. Followingearlierstudies 5 in the barrier for adatom-island detachment leads to a of adatom diffusion of similar systems,24,25,26 our super- 9 reduction in the mean island size, with an associated cell consists of two repetitions of the surface unit cell, 1 1 narrowing of the size distribution. These findings are depicted in Fig. 1 along the dimer rows direction (i.e. 0 qualitatively consistent with experimental observations5 a 4 × 4 supercell). Twelve layers of atoms were used, 3 in InAs/GaAs. Madhukar6 and Penev, Kratzer and separated by 6 layers of vacuum. By varying the thick- 0 Scheffler7,8 have considered the growth of InAs QDs on ness of the slab and the width of the vacuum separating / t GaAs, where increasing island size leads to a buildup of the periodic images perpendicular to the surface, it was a compressiveelasticstrainsinthesurroundingsubstrate.9 verified that the selected system size provides energies m Within simplified models of diffusion-limited growth withanaccuracyofthe orderof0.03eV.Singleadatoms - these authorsdemonstratedthat, inthe InAs/GaAshet- were placed on each of the two slab surfaces, which has d n eroepitaxialsystem,thestrain-dependenceoftheparam- the effect of doubling the adatom’s energetic contribu- o etersgoverningadatomdiffusiongivesrisetoareduction tion, therebyimprovingthe accuracyof the calculations. c in the flux of adatoms reaching larger islands relative to Allcalculationswereperformedusingtheab initio total- v: small ones, leading to a reduced rate of coarsening and energy and molecular-dynamics program VASP (Vienna i an associated narrowing of the size distribution. ab initio simulation package) developed at the Institut X While the potentially important consequences for is- fu¨r Materialphysik of the Universit¨at Wien27,28, which r land growth kinetics arising from strain dependencies implements Vanderbilt ultra-soft29 pseudo-potentials30. a in adatom binding and migration energies have been The energycutoff for the plane-wavebasis set wasset to clearly demonstrated, attempts to determine the mag- 150eVanda 2×2×1mesh ofk-points wasusedfor the nitude of these effects in specific systems have been un- 4×4 supercell. All atoms were allowed to relax. dertakeninrelativelyfewsystems.8,10,11,12,13,14 Theseef- On the Si (001) surface, the diffusion of a Ge adatom fects are most readily investigated using first-principles alongthe dimer rowsis knownto be of the orderof1000 calculations8,10,11,15,16, as they are extremely difficult to times faster than across dimer rows31 and we therefore isolate experimentally. The purpose of the present work focus solely on the diffusion along dimer rows. Since a is to investigate the effect of strain-dependent adatom typicalsurface consistsof terraceswhere the directionof ◦ diffusion and binding energies upon island growthkinet- thedimer rowschangesby90 ateachmonoatomicstep, 2 FIG.2: Geometriesofthebindingsite(left)andsaddlepoint (right)configurationsforaGeadatom(indicatedbyanarrow) diffusing on theGe (001) surfaces. FIG. 1: Geometries of thebindingsite and saddlepoint con- figurations for a Ge adatom (filled circle) diffusing on the Si (001) and Ge (001) surfaces. Only the two topmost atomic layers are shown. The size of each circle reflects the atom’s proximitytotheobserver. Inthelower left quadrant,crosses indicatetheadatompositionforthe6othercandidatebinding sites that were considered. fast diffusion in any direction is possible at the meso- scopic level, even though the fast diffusion is unidirec- tional at the microscopic level. Earlier studies24,25,26,32 haveunambiguouslydeterminedthelocationofthebind- ing site and the activated state of the Ge adatom on the Si(001)surface(seeFig. 1)aswellasthepreciseconfigu- rationofthesurfacedimersinthevicinityoftheadatom. FIG. 3: Strain-dependence of the binding energies E (left), b Since experimental and computational evidence is thesaddlepointenergiesEa+Eb (middle)andtheactivation scarcer in the case of the Ge adatom on the Ge (001) barriersEa (right)foraGeadatomdiffusingonSi(001)(top) surface,33 we considered 7 possible binding sites and andGe(001)(bottom)surfaces. Theseplotsareobtainedbya found the minimum energy site to be as shown in Fig- polynomial fit tothecalculated energies as afunction strain. ures 1 and 2. This site remains energetically favored for Let a, aSi and aGe respectively denote the substrate, the Si, values ofthe substrate lattice parameterupto 3%larger and theGe lattice parameters. than the Si lattice parameter. Beyond that threshold, a binding site analogous to the binding site of Ge on Si (001) becomes favorable. We shall neglect this alternate of the surface. The resulting energy versus strain re- binding site, as the substrate lattice parameter required lationships are plotted in Fig. 3. Our results reveal to stabilize it falls outside the range sampled during Ge four important observations. First, linear approxima- on Si (001) heteroepitaxy. tions to the strain-dependence of the binding energies We identifiedthe diffusionpathusingthe nudgedelas- andactivation barriers11 must clearly be used with care. tic bandmethod34,35 andrefinedthe positionof the sad- Second, the Ge adatom binding energy does not nec- dle point (see Figures 1 and 2) using the quasi-Newton essarily exhibit a minimum when the lattice parameter algorithm.28 We focused on the diffusion along the val- of the substrate matches the lattice parameter of bulk leys between the dimer rows, because the binding sites Ge. Third, the sign of the change of these energies locatedatopthedimerrowhadabindingenergyatleast under strain are highly system-dependent: calculations 0.3 eV larger than the saddle point energy for diffusion onthe structurallysimilarIn/GaAs(001)heteroepitaxial within the valleys, strongly suggesting that any other system8foundthestrain-dependencesofthebinding(Eb) saddle points located on the dimer rows would also have and saddle-point energies (Eb+Ea) the be of a sign op- a higher energy than the saddle point we analyzed. The positetotheonesintheGe/Ge(001)system. Thesignof nearly one-dimensional diffusion and the adatom’s pref- the strain-dependence of the activation barrier (Ea) for erence for sites located in the valley between the dimer Ge on both Si (001) and Ge (001), however, agrees with rows agrees qualitatively with the results of earlier cal- earlierstudiesinIn/GaAs(001)8andSi/Si(001)11,12,13,23. culations based on semi-empirical potentials,33 although Finally, the magnitude of the effect of strain on bind- thepreciselocationofthebinding sitesandofthesaddle ing and saddle point energies found in the case of Ge on points differ. Ge(001)isbyfarthelargest,relativetoanyothersystem For each type of surface (Si or Ge), the energies of studied so far.8,10,11,12,13,23 three geometries were calculated (the free surface, the In order to quantify the importance of these results in saddle point and the binding site configurations) at var- the context of quantum dot growth, we calculated the ious levels of biaxial strain imposed parallelto the plane strainfieldinthevicinityofaGe(105)-terminatedpyra- 3 Scheffler,8 employing a simple surface diffusion-limited one-dimensional model of adatom diffusion between two neighboringislands. Weexpandslightlyontheirworkby including boundary conditions at the island edges that impose local equilibrium between adatoms and islands. For simplicity, we consider all entropic prefactors to be strain-independent. We perform a stability analysis by considering two islands of equal size located at x = 0 and x = l and calculate the changes in the adatom flux toward each island as the size of one is perturbed. Let F ,F denote the flux towards each island and consider 0 l the (dimensionless) asymmetry in the flux towards each island, F = (F −F )/(2φl), where φ is the adatom de- l 0 position rate and l is the distance separating the two islands. Solving the above diffusion problem leads to FIG. 4: Strain field experienced by the substrate surface in the vicinity of a Ge island (dotted lines) of a typical size (10 K eβEi,l −eβEi,0 M nm)terminatedby(105)facets. Thehomogenouscomponent F =− − 1, (1) M φl(cid:0) l (cid:1) M of the in-plane strain relative to the Si lattice parameter is 0 0 plotted ((ǫxx+ǫyy)/2, for a surface normal to thez axis). where 1 l l n midal island with a characteristic size of 10 nm.36 For Mn = ln+1 Z (cid:18)x− 2(cid:19) eβ(Ea(x)+Eb(x))dx. (2) 0 the calculation of strain fields around coherent strained islands, both finite-element techniques (e.g., Ref. 9) as andwhereE (x)andE (x)are,respectively,theactiva- a b well as approximate analytical methods37 have been ap- tion barrier and binding energy as a function of position pliedpreviously. Inthe presentwork,elasticstrainfields x; E andE are,respectively,the energiesofanatom i,0 i,l have been derived from the results of the linear stabil- bound to the islands located at x = 0 and x = l; K is ity analysis of Spencer, Voorhees and Davis.38 In this a constant incorporating all entropic prefactors and β is analysis,the linearizedstrainfields arising from Fourier- reciprocal temperature (k T)−1. B mode perturbations in the surface height were derived Our analysis focuses on the diffusion of Ge adatoms within isotropic elasticity theory. By summing the re- on the Ge (001) surface, strained to match the lattice sulting expressions over the Fourier amplitudes describ- constant of the Si substrate, since a Ge wetting layer ing the shape function of a Ge pyramid, we calculated coveringthe SisubstrateisknowntoforminSKgrowth. the in-plane strain field (shown in Fig. 4) at the surface We assume that adatoms are insensitive to the presence of a Ge wetting layer (three atomic-planes in thickness) of Si under the Ge wetting layer, that no substantial Si- on Si(001). In Fig. 4 the strain is referenced to the Si Geinterdiffusionoccursandthatsurfacereconstructions lattice constant, so that zero corresponds to a wetting are unaffected by epitaxial strain. Accounting for the layer epitaxially strained on the Si substrate. In agree- latter (e.g., through the inclusion of a strain-dependent mentwithpreviouscalculationsforrelatedsystems,8,9,37 concentration of “missing dimers”39) may provide an- the strain field surrounding the island is compressive in other source of strain-dependent diffusion deserving fur- nature. ther consideration. In Fig. 4, the strain is seen to range in magni- Under the above assumptions, the quantity E (x)+ a tude from zero to roughly −1% in the vicinity of the E (x) decreases under a compressive strain, as shown b Ge(105)-terminated island. Over this range of compres- in the middle panels of Fig. 3. When the two islands sive strains, the binding energy of a Ge adatom on a are identical, F = 0, since E = E and M = 0 (as i,l i,0 1 Ge surface varies by about 0.21 eV, while the activa- E (x)+E (x)isthensymmetricwithrespecttox= l). a b 2 tion barrier varies by about 0.12 eV, as seen in Fig. 3. Now, consider an increase in the size of the island at ◦ At a typical deposition temperature of 600 C, these x=l. The first term of Equation (1) describes the ther- variations correspond, respectively, to a 16-fold varia- modynamic drivingforcefor coarsening: eβEi,l decreases tion in equilibrium adatom density and a 5-fold varia- duetocapillarity(alargerislandhasasmallersurfaceto tion in adatom mobility. In contrast, the correspond- volume ratio). This term thus causes F to become pos- ing strain-dependences on the Si (001) surface are much itive, increasing the flux towards the larger island. This less pronounced. Of course, the strain-dependence of term is inversely proportional to φ, demonstrating that the entropic prefactors8 could very well modify those increasing the deposition rate would act to reduce this order-of-magnitude estimates based solely on the strain- natural coarsening effect. dependence of the energetic contributions. The second term of Equation (1), which is indepen- To further explore the implications of our calculated dent of deposition flux, quantifies the effect of strain- results, we follow the analysis of Penev, Kratzer and dependent adatom binding and migration energies. A 4 slightincreasein the size ofthe islandat x=l induces a qualitative effect of accelerating the natural coarsening compressive strain in the wetting layer in the vicinity of rate of larger islands, in contrast with earlier results ob- that island, thus inducing a decrease in the saddle point tained for the related InAs on GaAs system.8 The large energyE (x)+E (x)inthesameareaandcausingM to magnitudeofthestraindependenceofadatomenergetics a b 1 become negative. The result is a relative increase in the onGe(001)obtainedinthepresentcalculationsalsocon- adatom flux towards the larger island. In other words, trasts with earlier findings in other systems8,10,11,12,13,23 this analysis suggeststhat the strain-dependencies ofGe and warrants further consideration of these effects in adatom binding and migration energies plotted in Fig. 3 more detailed kinetic models to further elucidate their act to accelerate the rate of coarsening of larger islands. effect on island growth. Insummary,first-principlescalculationshavebeenem- ployed to compute the strain dependence of Ge adatom binding energies and activation energies for diffusion on Acknowledgments both Si(001) and Ge(001) surfaces. For a growth tem- ◦ perature of 600 C, these strain-dependencies give rise to a 16-fold increase in adatom density and a 5-fold de- This work was supported by the NSF under programs crease in adatom diffusivity in the region of compressive DMR-0102794 and NSF-MRSEC DMR-00706097, using strain surrounding a Ge island with a characteristic size computer resources provided by the National Partner- of 10 nm. Within a simplified model of diffusion-limited ship for Advanced Computational Infrastructure at the growth, these strain dependencies are found to have the University of Michigan. ∗ [email protected]; http://www.mit.edu/~avdw/ 21 Y.-W.Mo,D.E.Savage,B.S.Swartzentruber,andM.G. 1 D.J.EagleshamandM.Cerullo,Phys.Rev.Lett.64,1943 Lagally, Phys.Rev. Lett. 65, 1020 (1990). (1990). 22 J.Tersoff,C.Teichert,andM.G.Lagally,Phys.Rev.Lett. 2 V.A.ShchukinandD.Bimbert,Rev.Mod.Phys.71,1125 76, 1675 (1996). (1999). 23 E. Zoethout, O. Gu¨rlu¨, H. J. W. Zandlivet, and 3 I. Daruka and A. L. Barab´asi, Phys. Rev. Lett. 79, 3708 B. Poelsema, Surf. Sci. 452, 247 (2000). (1997). 24 V. Milman, D. E. Jesson, S. J. Pennycook, M. C. Payne, 4 H.M.KoduvelyandA.Zangwill, Phys.Rev.B60,R2204 M. H. Lee, and I. Stich, Phys.Rev.B 50, 2663 (1994). (1999). 25 V.Milman, S.J.Pennycook,andD.E.Jesson,ThinSolid 5 N. P. Kobayashi, T. R. Ramachandran, P. Chen, and Films 272, 375 (1996). A. Madhukar,Appl.Phys. Lett. 68, 3299 (1996). 26 V. Milman, Int.J. of Quan. Chem. 61, 719 (1997). 6 A. Madhukar,J. Cryst. Growth 163, 149 (1996). 27 G. Kresse and J. Furthmu¨ller, Phys. Rev. B 54, 11169 7 P.Kratzer, E.Penev,and M.Scheffler,Appl.Phys.A 75, (1996). 79 (2002). 28 G. Kresse and J. Furthmu¨ller, Comput. Mat. Sci. 6, 15 8 E. Penev, P. Kratzer, and M. Scheffler, Phys. Rev. B 64, (1996). 085401 (2001). 29 D. Vanderbilt,Phys.Rev. B 41, 7892 (1990). 9 N. Moll, M. Scheffler, and E. Pehlke, Phys. Rev. B 58, 30 J. C. Phillips and L. Kleinman, Phys. Rev. 116, 287 4566 (1998). (1959). 10 C. Ratsch, A.P. Seitsonen, , and M. Scheffler, Phys. Rev. 31 Y. W.Mo and M. G. Lagally, Surf. Sci. 248, 313 (1991). B 55, 6750 (1997). 32 G.M.Dalpian,A.Fazzio,andA.J.R.daSilva,Phys.Rev. 11 D.J.Shu,F.Liu,andX.G.Gong,Phys.Rev.B64,245410 B 63, 205303 (2001). (2001). 33 K. Mea, Thin Solid Films 395, 235 (2001). 12 C. Roland and G. H. Gilmer, Phys. Rev. B 46, 13428 34 H. J. G. Mills and G. K. Schenter, Surf. Sci. 324, 305 (1992). (1995). 13 H. Spjut and D. A.Faux, Surf.Sci. 306, 233 (1994). 35 G. M. H. Jonsson and K. W. Jacobsen, in Classical 14 M. Schroederand D. E. Wolf, Surf. Sci. 375, 129 (1997). andQuantum DynamicsinCondensed Phase Simulations, 15 E. Kaxiras, Comput. Mat. Sci. 6, 158 (1996). editedbyB.J.Berne,G.Ciccotti,andD.F.Coker(World 16 A. P. Smith, J. K. Wiggs, H. Jonsson, H. Yan,L. R. Cor- Scientific, 1998). rales, P.Nachtigall, and K.D.Jordan, Journal OfChemi- 36 G. Medeiros-Ribeiro, A. M. Bratkovski, T. I. Kamins, cal Physics 102, 1044 (1995). D. A. A. Ohlberg, and R. S. Williams, Science 279, 353 17 F.M.Ross,J.Tersoff,andR.M.Tromp,Phys.Rev.Lett. (1998). 80, 984 (1998). 37 J. Tersoff and R. M. Tromp, Phys. Rev. Lett. 70, 2782 18 F.M.Ross,R.M.Tromp,andM.C.Reuter,Science286, (1993). 1931 (1999). 38 B. J. Spencer, P. W. Voorhees, and S. H. Davis, J. Appl. 19 R. S. Williams, G. Medeiros-Ribeiro, T. I. Kamins, and Phys. 73, 4955 (1993). D.A.A.Ohlberg,Annu.Rev.Phys.Chem.51,527(2000), 39 F. Liu, F. Wu, and M. G. Lagally, Chem. Rev. 97, 1045 and references therein. (1997). 20 J.A.Floro, E.Chason,R.D.Twesten, R.Q.Hwang, and L. B. Freund,Phys. Rev.Lett. 79, 3946 (1997).