Island Distance in One–Dimensional Epitaxial Growth Harald Kallabis1,2, Paul L. Krapivsky2 and Dietrich E. Wolf1,2 arXiv:cond-mat/9806140v3 [cond-mat.stat-mech] 17 Sep 1998 1 FB 10, Theoretische Physik, Gerhard–Mercator–Universität Duisburg, 47048 Duisburg, Germany 2 Center for Polymer Studies, Boston University, Boston, 02215 MA, USA January 9, 2014 The typical island distance ℓ in submonlayer epitaxial growth depends on the growth conditions via an exponent γ. This exponent is known to depend on the substrate dimensionality, the dimension of the islands, and the size i∗ of the critical nucleus for island formation. In this paper we study the dependence of γ on i∗ in one–dimensional epitaxial growth. We derive that γ = i∗ /(2i∗ + 3) for i∗ ≥ 2 and confirm this result by computer simulations. PACS numbers: 81.15-z, 81.10.Bk, 05.50.+q I. INTRODUCTION γ= Molecular beam epitaxy (MBE) is of outstanding importance for the fabrication of microelectronic devices. The highly controllable experimental conditions make this technique an ideal tool for crystal growth. The thickness of the deposited layers can be controlled down to the atomic scale while growth proceeds layerwise. Typically the growth of multiple layers is of primary interest in most experimental and theoretical studies of MBE. However, deeper insight into the physics of multilayer growth requires a comprehensive understanding of the submonolayer regime (see Ref. [1], for instance). The most important length scale in submonolayer growth is the island distance or diffusion length ℓ. The classical result [2] for this length is that it depends on the diffusion constant D of adatoms and the deposition rate F like ℓ ∼ (D/F )γ . i∗ . 2i∗ + d + df (3) Our motivation for studying the one-dimensional case is the following. Much theoretical work is dedicated to the study of basic mechanisms of epitaxial growth. For computational reasons, the numerical work in this area is often restricted to models in surface dimension d = 1. As mentioned above, such studies rely on a good understanding of the physics in the submonolayer regime. Also, quasi one–dimensional MBE may be realized experimentally on surfaces with extremely anisotropic diffusion constant, like the 2x1–dimer–reconstructed Si(001) surface (see, e.g., Ref. [8]). There, the estimated ratio of the diffusion constants along rows and perpendicular to them is of the order of 1000:1 in typical experiments [9]. A computer simulation study of this growth process leads to the estimate that for a deposition rate of F = 1/600 ML/s and at temperatures below T = 450 K, this anisotropy should be pronounced enough to lead to quasi one–dimensional MBE [10]. However, there is a controversy concerning the relative importance of anisotropic diffusion and anisotropic bonding [11,12]. Another experimental example is the reconstructed surface of Au(100) [13]. This paper is organized as follows. In section II the expression (3) will be derived. The subsequent sections are devoted to the numerical confirmation of this result. In section III we employ the coarse–grained MBE (CGMBE) model [3], which allows for a simple algorithmic implementation of a variable critical nucleus size i∗ . To ensure that the results are independent of the chosen algorithm, we reproduce them with a model that is not coarse grained (MBE model) in section IV. It is shown that the results of the MBE and the CGMBE models compare favourably with each other and with the theoretical prediction. The final section V contains a general discussion of the results. (1) The analysis of rate equations yields that the exponent γ depends on the substrate dimension d, on the size i∗ of the critical nucleus and the fractal dimension df of the islands [2–4]: γ= i∗ . 2i∗ + 3 (2) The critical nucleus is defined such that all islands consisting of i∗ or less atoms can shrink by emitting adatoms, while an island of size i∗ + 1 is stable. Eq. (2) applies if island diffusion and adatom desorption are not important. Eq. (2) has been demonstrated to hold in d = 2 for i∗ = 1, 2, 3, 4, see Refs. [5,6]. In this paper we study the dependence (1) in surface dimension d = 1 for different values of the critical nucleus i∗ . In this case, df = 1. So far, only the case i∗ = 1 [3,7] has been studied in one dimension, to our knowledge, and it is in agreement with (2). Here we show that (2) is no longer valid for i∗ ≥ 2, where 1 Setting t ∝ F −1 gives the result (3). For i∗ = 2 there is a logarithmic correction to (6) [15], II. ANALYTICAL CONSIDERATIONS Most treatments of the monolayer growth processes are based on a mean-field rate equation approach. This is appropriate in the most important two–dimensional case since the basic kinetic mechanism – the two-particle diffusion-controlled aggregation – has an upper critical dimension dc = 2 [14]. In one dimension, however, classic rate equations do not apply and we shall employ a modified rate equation approach [7]. We denote by A and I densities of adatoms and stable immobile islands, respectively. We start with adatoms and write dA A =F − . dt τ DA3 dI = , dt | ln I| whereas (5) remains unchanged. Repeating the same analysis gives ℓ∼ D F 2/7  ln  D F 1/7 . (9) III. NUMERICAL RESULTS The idea of the coarse grained MBE (CGMBE) model, introduced in Ref. [3], successfully applied in Refs. [5,6], and studied in detail in Ref. [16], is to avoid the timeconsuming simulation of the surface diffusion on the atomic scale by resolving it only on a much coarser scale ∆x. As long as ∆x is smaller than the island distance ℓ, one still gets correct information about the surface morphology. Specificially the model is implemented as follows: The substrate is divided into cells of linear extension ∆x, so that a monolayer corresponds to ∆x atoms per cell. The partial filling of a cell at site x is given by a counter m(x) = 0, 1, . . . , ∆x − 1, which means the number of atoms on top of completed monolayers in cell x. The height h(x) denotes the number of complete monolayers. For the complete description of the state of a cell, an edge indicator is(x) is introduced. If is(x) = 1, atoms cannot leave cell x because they are irreversibly incorporated into the island edge present in the cell. This allows for a simple implementation of the concept of the critical nucleus i∗ : if upon diffusion or deposition of an atom into cell x the counter reaches the value m(x) = i∗ , is(x) is set to one. Consequently, m(x) can only increase from then on, until m(x) = ∆x, in which case h(x) is increased by one and m(x) and is(x) are set to zero. As long as is(x) = 0, m(x) can also decrease, when atoms diffuse into neighbouring cells. The dynamics has to take into account that a single adatom move from one cell to a neighbouring one represents the diffusion over a distance of ∆x lattice sites. This happens with frequency νD = D/(∆x)2 for all adatoms in the system. The deposition of one atom happens with frequency νF = F L∆x. Hence the probability for deposition during one time step ∆t = (νF + νD )−1 of the computer simulation is (5) This differs from the two–dimensional situation where the number of different sites visited by an adatom grows linearly with time and the adatom density evolves according to dA/dt = F − DIA. Actually, Eq. (5) contains an additional loss term due to the nucleation of new islands. However, this loss term is asymptotically subdominant, since most adatoms are captured by the existing islands. The rate equation for the island density is derived by noting that islands are formed upon colliding of i∗ + 1 adatoms. Hence ∗ dI = DAi +1 . dt  (4) Here τ is a lifetime, i.e., the time between the deposition of an adatom and its incorporation into a stable √ island. During the time interval τ , the adatom visits Dτ different sites, so the adatom lifetime is determined by √ I Dτ ≈ 1, and Eq. (4) becomes dA = F − DI 2 A. dt (8) (6) The many-particle diffusion-controlled reaction process (i∗ + 1)A → I has an upper critical dimension dc = 2/i∗ (see e.g. [15] and references therein). The fact that the upper critical dimension gets smaller when i∗ increases is simple to appreciate. Indeed, many-particle collisions occur rarely, so the system is well mixed and the meanfield approach should become correct earlier. For i∗ > 2 we have dc < 1, so Eq. (6) applies; the case i∗ = 2 is marginal, so we anticipate logarithmic corrections to the mean-field results. For i∗ = 1 correlations dominate the diffusion process in one dimension, so that (6) has to be modified [14] and leads to γ = 1/4 [7]. Eq.(3) is obtained from (6) and (5) in the following way: From (5) it follows that asymptotically A ≃ (F/D)I −2 . Inserting this into (6) and intergrating we find  i∗ /(2i∗ +3) ∗ F (F t)1/(2i +3) . (7) I∼ D p = νF /(νF + νD ). (10) 1 − p is the probability with which each of the adatoms is allowed to diffuse into a neighboring cell. If no adatoms 2 are present, νD = 0. This algorithm speeds the simulation up by a factor of about (∆x)2 . In the simulations presented in the following, the coarse graining length had the constant value ∆x = 20. The system size was chosen to be L∆x = 2 · 105 and averages were taken over 100 independent runs. model was simulated in systems of size L = 65536. The island distance was measured here as follows. During deposition, the number of islands was calculated for times 0 < F t < 1. An island was defined as a connected region where h(x) ≥ a⊥ , where a⊥ is the vertical lattice constant. For early times, the number of islands increases due to nucleation. At roughly F t ≃ 0.5 the number of islands decreases again since the islands coalesce. This means that the island distance has a minimum which was used for evaluation. 4 10 * i =1 * i =2 * i =3 * i =4 4 10 10 0.4 0.3 γeff 3 10 island distance island distance 3 2 10 0.2 0.1 0.0 0 2 2 10 4 6 log(D/F) 8 10 1 10 4 10 6 8 10 10 10 10 D/F 1 10 FIG. 1. The island distance as a function of D/F for i∗ = 1, 2, 3, 4. The solid lines represent fits to the last three decades of data and have the slopes γ = 0.248, 0.293, 0.329, 0.355, respectively. 10 0 2 10 4 6 10 10 8 10 10 10 D/F FIG. 2. The island distance as a function of D/F for i∗ = 1 (circles) and i∗ = 2 (squares). The solid lines with slopes γ = 0.246 and γ = 0.288 represent fits to the last three decades of data. The inset shows how the local slopes approach the theoretically predicted exponents (dotted lines) for large D/F . In the CGMBE model, nucleation events can be identified easily. Thus the island distance can be determined by dividing the system size by the total number of nucleation events in the first layer. Measurements of this quantity as a function of D/F for different i∗ are shown in Fig. 1. Note that the asymptotic algebraic behaviour ℓ ∼ (D/F )γ is reached with ℓ as function of D/F having a weak negative curvature. In Fig. 2, the measurements of the island distance as a function of D/F for the MBE model with i∗ = 1 and i∗ = 2 are shown. From the effective exponent γeff ≡ log(10 · D/F ) − log(D/F ) shown in the inset, it becomes clear that the asymptotic behaviour ℓ ∼ (D/F )γ is reached only for large values of D/F . The dotted lines indicate the theoretical values from Eq. (2), γ(i∗ = 1) = 1/4, and Eq. (3), γ(i∗ = 2) = 2/7, respectively. Note also that the asymptotic exponent is approached from below. IV. RESULTS WITHOUT COARSE GRAINING In the (non–coarse grained) MBE model atoms are deposited onto the surface of a simple cubic crystal (i.e. a square lattice in d = 1) with a rate of F atoms per unit time and area. Atoms with no lateral neighbors are allowed to diffuse with diffusion constant D. In the simplest variant of the model atoms with lateral neighbors are assumed to be immobile so that, e.g., dimers are immobile and stable. Therefore, i∗ = 1. Extension of the rules from i∗ = 1 to i∗ = 2 is straightforward: Dimers are allowed to split into two adatoms, while trimers are stable. In either case, growth commences with a flat substrate, h(x, 0) = 0 for all sites x. On deposition at x, h(x, t) is increased by one. We neglect barriers to interlayer transport (Ehrlich–Schwoebel barriers [17]). Thus the only parameter of the model is the ratio D/F . This V. DISCUSSION In Fig. 3, the comparison of the simulation results of the CGMBE model (Fig. 1) and the MBE model (Fig. 2) with various analytical formulas is shown. The small difference between γ measured in the CGMBE and the MBE model is not only due to statistical noise, but also due to the asymptotic behaviour of ℓ(D/F ) being reached with negative and positive curvature, respectively. Therefore, even better agreement between the results of the CGMBE and the MBE models for larger D/F can be 3 expected. As a consequence, the values for higher i∗ , obtained with the CGMBE model, can be regarded as reliable. length in one dimension for critical island size i∗ ≥ 2 and have confirmed the result by computer simulations. VI. ACKNOWLEDGEMENTS 0.40 D. E. W. acknowledges support by DFG within SFB 166. H. K. acknowledges support by DAAD within the Hochschulsonderprogramm III. P. L. K. acknowledges NSF grant DMR-9632059 and ARO grant DAAH04-961-0114 for financial support. γ 0.35 0.30 VII. NOTE ADDED IN PROOF 0.25 0.20 1 2 3 Following the suggestion of one of the referees we note that our rate equation approach applies to the general case of surface dimension d and island dimension df . The fractal island dimension can be smaller or larger than the surface dimension which is integer in physically relevant cases. In this general case we use F tℓd = ℓdf with ℓ = I −1/d to obtain 4 * i FIG. 3. The submonolayer exponent γ as a function of i∗ in d = 1. Squares and diamonds are the results of the simulations of the CGMBE and the MBE models, respectively. The solid curve is the theoretical prediction for i∗ ≥ 2, Eq. (3). The dashed curve corresponds to formula (2) and the dotted line are the exponents expected, if islands up to size i∗ can diffuse, Eq. (11). γ= where ν = 2 for d = 1 and ν = 1 for d ≥ 2. This does not apply to the exceptional case d = i∗ = 1 where d < dc = 2/i∗ . There equation (6) should read I˙ = DA2 I and one finds γ = 1/(3 + df ). The agreement of the numerical results with the theoretical prediction (3) is very good. For comparison, we also show Eq. (2), which in one dimension applies only for i∗ = 1. It is interesting to notice that the exponent for i∗ = 2 coincides with the one expected, if islands up to size i∗ can diffuse without decaying [3,18] γ= i∗ i∗ = ∗ . ∗ (d + 2)i + df 3i + 1 i∗ , dν(i∗ + 1) + df (11) [1] H. Kallabis, L. Brendel, J. Krug and D. E. Wolf, Int. J. Mod. Phys. B 11, (1997). [2] G. Zinsmeister, Thin Solid Films 2, 497 (1968); ibid. 4, 363 (1969); ibid. 7, 51 (1971); S. Stoyanov and D. Kashchiev in E. Kaldis (ed.): Current Topics in Material Science, Vol. 7 (North-Holland, Amsterdam, 1981) p. 69; J. A. Venables, G. D. Spiller and M. Hannbrucken, Rep. Prog. Phys. 47, 300 (1984); J. Villain, A. Pimpinelli and D. Wolf, Comments Cond. Mat. Phys. 16, 1 (1992). [3] D. E. Wolf, in: Scale Invariance, Interfaces and NonEquilibrium Dynamics, A. McKane et al. (eds.) (Plenum Press, New York, 1995) pp. 215 – 248. [4] D. E. Wolf, in: Dynamics of Fluctuating Interfaces and Related Phenomena, D. Kim, H. Park and B. Kahng (eds.) (World Scientific, Singapore, 1997) pp. 173 – 205. [5] M. Schroeder and D. E. Wolf, Phys. Rev. Lett. 74, 2062 (1995). [6] H. Jeong, B. Kahng, and D. E. Wolf, Physica A 245, 355 (1997). [7] A. Pimpinelli, J. Villain, D. E. Wolf, Phys. Rev. Lett. 69, 985 (1992). In fact, unstable islands effectively diffuse in our models. Consider for example i∗ = 2, where dimers are unstable. When a dimer decays, the adatoms still remain close to each other, and in one dimension the probability is particularly high that they reunite again at a shifted position. The specific way we implemented the decay of a cluster in our computer simulations makes this tendency very clear: Whenever a dimer has formed it must break up in the next diffusion step, which means that each of the two atoms moves to one of the neighboring cells. With probability 1/2 they will be found in the same cell. If no further adatom was met there, this amounts simply to a displacement of the dimer by ±∆x. This analogy breaks down, however, for i∗ > 2. In conclusion, we have investigated the submonolayer epitaxial growth in a one dimensional model where islands of size up to i∗ are unstable while larger islands are stable and immobile. We have derived the exponent γ which determines the D/F –dependence of the diffusion 4 [8] M. Lagally, Y.–W. Mo, R. Kariotis, B. S. Schwarzentruber, and M. B. Webb, in: Kinetics of Ordering and Growth at Surfaces, M. G. Lagally (ed.) (Plenum, New York, 1990) p. 145. [9] Y. W. Mo, J. Kleiner, M. B. Webb, and M. G. Lagally, Phys. Rev. Lett. 66, 1998 (1991); Surf. Sci. 268, 275 (1992). [10] M. C. Bartelt and J. W. Evans, Europhys. Lett. 21, 99 (1993). [11] S. Clarke, M. R. Wilby and D. D. Vvedensky, Surf. Sci. 255, 91 (1991). [12] C. Pearson, M. Krueger, and E. Ganz, Phys. Rev. Lett. 76, 2306 (1996). [13] S. Günther, E. Kopatzki, M. C. Bartelt, J. W. Evans, and R. J. Behm, Phys. Rev. Lett. 73, 553 (1994). [14] For a recent review, see S. Redner, in Nonequilibrium Statistical Mechanics in One Dimension, ed. V. Privman (New York: Cambridge University Press, 1997), p. 3. [15] P. L. Krapivsky, Phys. Rev. E 49, 3233 (1994). [16] H. Kallabis, Ph.D. dissertation, Universität Duisburg (1997). [17] G. Ehrlich and F. G. Hudda, J. Chem. Phys. 44, 1039 (1966); R. L. Schwoebel and E. J. Shipsey, J. Appl. Phys. 37, 3682 (1966). [18] I. Furman and O. Biham, Phys. Rev. B 55, 7917 (1997). 5