Si (001) Shallow

NANO LETTERS

Calculation of Vertical Correlation Probability in Ge/Si(001) Shallow Island Quantum Dot Multilayer Systems

2006 Vol. 6, No. 6 1279-1283

Maxim A. Makeev* Collaboratory for AdVanced Computing and Simulations, Mork Family Department of Chemical Engineering and Materials Science, UniVersity of Southern California, Los Angeles, California 90089-0242

Anupam Madhukar Nanostructure Materials & DeVices Laboratory, Mork Family Department of Chemical Engineering and Materials Science, and Department of Physics, UniVersity of Southern California, Los Angeles, California 90089-0241 Received February 3, 2006; Revised Manuscript Received April 10, 2006

ABSTRACT We investigate the behavior of the island vertical pairing probability in multilayer systems of Ge island quantum dots (QDs) in Si(001). By combining a simple kinetic rate model with our previously reported atomistic simulation results on the nature of the stress field from buried shallow Ge islands having {105}-oriented sidewalls, we derive an analytical expression for correlation probability as a function of the parameters characterizing the multi-QD systems. The approach is based upon continuum mechanochemical potential model, which allows one to introduce necessary elements of the kinetics of island formation in a simple way. We compare the model predictions with available experimental data and find that the model provides a satisfactory description of the coupling probability. The correlation probability behavior as a function of capping layer thickness, Ge island size, interisland distance, and Ge adatom diffusion length is investigated within the framework of the developed model.

The mechanism of vertical correlation in the multilayer systems of buried heteroepitaxial islands has been under extensive examination recently.1,2 The interest is motivated by the potential applications of the lattice-mismatched coherent three-dimensional islands3,4 as quantum dot (QD) building blocks for devices in the fields of photonics and nanoelectronics.5 Moreover, in the multilayer structures of QDs, besides vertical self-organization, the islands can gain increased size and shape uniformity from layer to layer in the resultant island superstructure.6 Utilizing the convenient analytical expression for the stress due to a spherical inclusion in a solid medium, ref 1 modeled the effect of buried island-induced stress field on adatom diffusivity and computed the vertical pair correlation probability as a function of the capping layer thickness. The theoretical model was used to fit the calculated vertical pair correlation probabilities to those obtained experimentally for the InAs/ GaAs QD systems. Furthermore, experimental studies of the QD multilayer structures have been performed for a number of realizations of InGaAs/GaAs,7-9 Ge/Si,10-14 and PbSe/ PbEu1-xTex15 systems. In the lattermost case a hexagonal * Corresponding author. E-mail: [email protected]. 10.1021/nl0602600 CCC: $33.50 Published on Web 05/13/2006

© 2006 American Chemical Society

ordering was observed, which was attributed to the existence of threefold symmetry elastic energy minima, each of which serves as a sink for adatoms.15 A very general understanding of the process of self-organization can be reduced to a simple realization that the buried islands induce inhomogeneous stress profiles on the spacer layer surface and the associated stress gradient leads to the adatom diffusion bias during a subsequent layer deposition.1 A nearly perfect vertical ordering was observed for small spacer thickness values, while nearly random positioning takes place when spacer layer thickness is large and, thus, the stress coupling is weak.1 In ref 1, the notion of small or large spacer layer thickness is in relation to the inverse cubic decay of the stress due to a spherical inclusion. Subsequently, the inverse cubic dependence model of the island-induced stress of ref 1 has been used to analyze the experimentally determined correlation probability in Ge/Si island multilayer structures.10 The qualitative features, observed for probability versus spacer layer thickness behavior for the InAs/GaAs and Ge/Si(001) QD systems, were found to be quite similar to each other, except that the characteristic lengths of probability decay differed significantly: that of

(WL) is positioned on the spacer layer surface and has thickness hW. At the basis of our theoretical approach is the mechanochemical potential model, used previously for the purposes of correlation probability computations in ref 1. We concern ourselves with the problem of second row island formation on the spacer layer surface, when a buried island (acting as a built-in stressor) is present in the system. We assume that the driving force of vertical self-organization is the stress-gradient-driven diffusivity. The chemical potential for Ge adatoms diffusing on the Ge WL surface under the stress gradient can be written as µ(x,y,z ) zs) ) µ0 + (Ω/2EGe)σ2τ + γΩR - Ωσnn

Figure 1. Schematic illustration of the model system under consideration. The notations are explained in the text.

the InAs/GaAs systems being smaller as compared to the Ge/Si systems. In the latter case, the pairing probability was found to persist for lengths nearly three times larger. Recently, the behavior of the stress field for the {105}faceted Ge island systems was uncovered via atomistic simulations utilizing Stillinger-Weber model of silicon and germanium.16,17 It was found that, for shallow {105}-faceted Ge/Si(001) QDs, the stress dependence on the spacer layer thickness follows inverse linear dependence in contrast to the 1/hsp3 behavior for spherical inclusions. However, the much steeper {101}-faceted InAs islands in GaAs show nearly 1/hsp3 behavior. This calls for a reexamination of the Ge/Si island vertical correlation data.10,14 It should be emphasized, however, that besides the stress field distribution a number of other factors can also contribute to the mechanism of vertical coupling, including various kinetic rates.1 In this work, we develop a simple model for the correlation probability by combining essential elements of the mechanochemical potential consideration1 and atomistic simulation results reported previously by us on stress behavior with spacer layer thickness and lateral stress field variations.16,17 The main virtue of our approach is that it allows for investigation of the buried island’s geometry effects on the stress distribution, with geometrical factors being incorporated in the model through the scaling relation derived in refs 16 and 17. Consequently, our model can be used to study the pairing probability in the multilayer Ge/Si(001) QD systems as a function of nearly all of the relevant experimental parameters. This is unlike the models based upon reducing the extended geometry of an inclusion to the effective pointlike force dipole.1 A schematic illustration of the model system geometry is shown in Figure 1. Note that the model we consider is a 2D model (x-z), with Si(001) substrate being along the x axis. The unit cell (double unit cell is shown in Figure 1) consists of Ge buried island of lateral size d and vertical height h, covered with a spacer layer of thickness hsp. The unit cell is repeated periodically in the x direction and, thus, the interisland distance in the systems is l. A Ge wetting layer 1280

(1)

where µ0 is the surface chemical potential for an unstrained surface, στ is the tangential part of the surface stress, R is the local surface curvature, σnn is the normal surface stress, Ω is the atomic volume of Ge atoms, and zs ) hsp+ hW. The normal stress at a free surface is negligibly small (as also confirmed by our atomistic simulation studies, reported in refs 16 and 17) and can be neglected. The surface curvature was the focus of our attention in ref 17. It was shown that the contribution from this term can also be neglected. Within the framework of our model, we consider the island growth as a process of stress-gradient-driven motion of adatoms to the region -d/2 < x < d/2 on the Ge WL surface. The probability that an adatom will be driven to this region can be computed as K)

D µ(ls,zs) - µ(0,zs) kBT ls

(2)

where D is the temperature-dependent coefficient of surface diffusion of Ge adatoms, kB is the Boltzmann constant, T is growth temperature, and ls is the spacer-layer-thicknessdependent range of the surface stress field. Following ref 1, we proceed further by invoking the classical model for surface adatom concentration during deposition.18 At the borders of a free-standing island, the adatom concentration can be computed as follows n((d/2) ) FLD

sinh Q K cosh Q + (LD/τ) sinh Q

(3)

where F is the incident atom flux, τ is the effective diffusion time, LD ) (Dτ)1/2 is the diffusion length, and Q ) (l d)/2LD. We assume that the range of the stress field is larger than half of the distance between the islands and, thus, use l/2 as the average stress field range.1 The amount of material arriving to the surface region with extent of the island lateral size per unit time, due to both driven adatom diffusivity and deposition from the vapor phase, and the net material amount arriving to the spacer layer surface region per unit time can be computed, respectively, as M1 ) Fd + Kn(d/2) + Kn(-d/2)

(4)

Nano Lett., Vol. 6, No. 6, 2006

and MT ) Fl

(5)

The correlation probability is defined as PV ) M1/MT. After some simple algebra, the following expression for the correlation probability can be obtained readily PV )

d + ls

(LD/ls) sinh Q ls kBTEGe 1 cosh Q + sinh Q LD Ω Σσ

[( )( ) ]

(6)

where we have defined Σσ ) δσ2τ (l/2,zs) - δσ2τ (0,zs)

(7)

The effective stress at the surface of the system, δστ, is defined as a difference between the Ge WL associated stress and the stress induced by the buried island at the Si(001) spacer layer surface, that is δστ(x,zs) ) σGe - στ(x,zs)

(8)

The above expressions contain geometrical factors and the effective tangential stress distribution along the x direction, which can be deduced readily from our simulation data reported previously.16,17 As we have shown in refs 16 and 17, the hydrostatic stress (which is approximately tangential at the surface), can be written in the following form Tr(σ) ) Psim(d simhsim)

1 zsim s

Figure 2. Vertical correlation probability is shown as a function of the spacer layer thickness as extracted from experimental studies: solid squares correspond to data from Figure 4 of ref 10; solid circles correspond to data from Figures 2 and 3 of ref 21. Dashed and solid lines represent our model predictions for the two experimental data sets. Dotted and dash-dotted lines represent the fitting to spherical inclusion model for the same sets of parameters: dotted line corresponds to data from Figure 4 of ref 10; dashdotted line corresponds to data from Figures 2 and 3 of ref 21. The two values of diffusion length are the same as those obtained for the best 1/zs fit.

By combining the sets of eqs 6-8 and eq 10, and computing numerical factors involved, we arrive at the following expression for the coupling probability d PV ) + l

(9)

cosh Q +

2(LD/l) sinh Q 0.08 × Tzs l LD d × h × f(l,zs)

[( )(

)]

(11) sinh Q

where where the subscript “sim” corresponds to parameter values used in our simulations.16,17 Usage of the simulation data allows one to obtain the stress field from any Ge {105}faceted buried island of pyramidal geometry in Si(001) as

στ(x ) 0,zs) ≈ σsim τ

( )( ) zsim dh s d simhsim zs

(10)

where the scaling of the stress field with surface area of the buried island is used.16,17 Note that the simulation data have been shown previously to provide a very satisfactory quantitative agreement with the analytical models, when such models are used appropriately.19 In the following, we use numerical value 0.1 eV/Ω for σsim τ , which corresponds to the spacer layer thickness of 41 ML. The buried island geometry (hsim, dsim) is exactly as in refs 16 and 17. The lateral variations of the stress field (the x dependence) can be approximated as P1(cos(θ)), where θ is the angle between normal to the capping layer surface and a vector that connects island center and corresponding point on the surface, where the stress field is computed. Nano Lett., Vol. 6, No. 6, 2006

f(l,zs) ≈

x(l2/4z2s + 1) x(l2/4z2s + 1) - 1

(12)

In the above expressions, temperature is measured in Kelvin and all lengths are in nanometers. The diffusion length is a fitting parameter because it cannot be extracted from the results of experimental studies easily. However, some estimates of reasonable range of variation of this quantity can be made based upon the findings of ref 20, where a lower bound of 700 ( 100 nm for T ) 700 K was reported. We use this value to roughly define the range of LD variations in our model. Consequently, knowing the temperature, average lateral and vertical sizes of an island, and island density (or interisland distance) one can compute the correlation probability readily using eqs 11-12. Moreover, experimental data reported previously10,14,21 allow for a direct comparison of our model predictions with pairing probability behavior observed in real systems. The normalized vertical correlation probability, P ) (PV - PR)/(1 - PR), is shown in Figure 2 as a function of the 1281

spacer layer thickness for two different buried island systems, parametrized as explained below (solid line and dashed line). The average random overlap probability is computed as PR ) d/l. Also shown in the figure are the experimental data from ref 10 (solid squares) and ref 21 (solid circles). As one can observe, the two computed probability functions compare well with the experimental data, if the diffusion lengths are chosen as LD ) 550 nm (T ) 893 K) and LD ) 1100 nm (T ) 973 K) for the two corresponding experimental data sets. The other parameters are taken as reported in the corresponding experimental studies. Two characteristic features pertaining to our model should be emphasized. For small values of the spacer layer thickness, the model underestimates the correlation probability. However, it overestimates the experimental data slightly, in the region where nearly random positioning is expected. In general, the model captures the general features of behavior well, avoiding any unphysical effects in the fitting procedure. Note that, in experimental studies, the correlation probability was obtained using TEM,14 Raman scattering, and X-ray diffraction methods.21,22 The comparison of the results obtained via these methods is given in refs 21 and 22. The largest difference in the results, taking into account the error bars, can be estimated as ∼10%. It is also worth noting that differences between model predictions and real systems behavior are expected because during capping layer deposition at elevated temperatures buried islands can undergo an uncontrollable shape transformation from perfect to truncated pyramid geometry. This transformation leads to corresponding changes in the nature of stress behavior, with the inverse power law exponent increasing from ∼1 (pyramid) to ∼3 (truncated pyramid). Thus, the trend observed by us is fully consistent with experimental observations. For completeness, we have also used the spherical inclusion approximation (eq 5 of ref 1) to fit the experimental data. The results of fitting to this 1/z3s dependence (shown in Figure 2 as dotted and dash-dotted lines) demonstrate that the spherical inclusion approximation does not work well in the case of {105}-faceted shallow Ge islands. In fact, if one attempts to obtain the closest agreement by varying LD, the value of diffusion length is severely underestimated.20 Indeed, while the functional form of the stress behavior in the case of truncated pyramid might be similar to the case of spherical inclusion, the strength of the field is quite different. In general, the 2D model we consider should not be expected to provide a complete quantitative description of the phenomenon, although, as we have shown, it is quite adequate in reproducing the major qualitative features thereof. The behavior of the correlation probability as a function of the spacer layer thickness is shown in Figure 3 for six different values of the diffusion length. The lateral and vertical sizes of the buried island are 60 and 6 nm, respectively. The island facet orientation is always kept at {105}. The distance between two islands is ∼125 nm. These values are chosen to reproduce the experimental data, reported in refs 14, 21, and 22. As one can observe, the probability function has three different regimes of behavior. For small values of zs, it decays slowly and, after a certain 1282

Figure 3. Correlation probability is shown as a function of the spacer layer thickness for different values of the diffusion length, LD. The island sizes and the interisland distance are as described in the text.

Figure 4. Vertical correlation probability is plotted as a function of the island size, h, for different values of the diffusion length, LD. The spacer layer thickness is fixed at hsp ) 50 nm.

(other parameters dependent) value of the spacer layer thickness is exceeded, a rapid fall off is observed. For large spacer layer thickness, the probability function approaches zero asymptotically. It should be emphasized that the behavior of the probability is dependent strongly on the adatom diffusivity. An increase in LD leads to considerably slower decay of paring probability and extends the plateau (slow decay regime), observed for small values of the spacer layer thickness. Figure 4 shows the correlation probability as a function of the vertical size, h, of the buried island, computed for four different diffusion lengths. The lateral size of the island is chosen to keep the {105} facet orientation fixed. In the whole range of h variations, the probability of correlation increases with increasing island height. As can be seen in the figure, in the intermediate range of island sizes, LD plays a decisive role in shaping the correlation probability. However, for large island sizes, the variations of the diffusion length from 500 to 1200 nm do not change the probability Nano Lett., Vol. 6, No. 6, 2006

In summary, using the mechanochemical potential-based model for driven adatom diffusivity, we have studied the behavior of the correlation probability in the multilayer systems of Ge/Si(001) islands. By incorporating our previously reported simulation results into the model, we have derived an analytical expression, which can be used to predict the behavior of the correlation probability for systems with arbitrary Ge island sizes and densities, provided the diffusion length for Ge adatom migration on the surface of the sample can be estimated accurately. The expression was used to investigate the behavior of the vertical correlation probability as a function of the Si spacer layer thickness, Ge island size, and the spacer-layer-thickness-dependent stress field ranges.

Figure 5. Vertical correlation probability is shown as a function of the interisland separation distance for different values of the diffusion length, LD. The spacer layer thickness is hsp ) 50 nm.

function much. Indeed, an increase in the island size leads to effective increase in the spacer surface stress and, consequently, leads to enhancement of the correlation between buried and newly formed islands, provided adatoms on the surface have sufficient mobility to arrive to the corresponding regions on the spacer surface within finite time intervals. The effect of the interisland separation, l, on the vertical coupling probability is illustrated in Figure 5. Note that the impact of l should be considered in relation to ls, the lateral range of the stress field at the spacer layer surface. This range, defined as the tensile region extent, for the Ge island sizes and spacer thicknesses considered here, was estimated as ∼497 nm. Therefore, we concentrate on the regime, where the stress field range exceeds half the interisland distance. In this regime, the pairing probability is controlled by two factors. First is the purely statistical factor d/l. The second is the contribution of the interisland-separation-distancedependent stress-gradient-driven diffusivity, which enhances the pairing probability. This enhancement is also controlled by adatom thermal diffusivity through LD. Note that monotonic dependence on the diffusion length is a very encouraging fact here, because LD is the fitting parameter of the model. Indeed, given the observed behavior, the model is expected to provide robust qualitative results, even if the diffusion length cannot be obtained with a proper precision. For large interisland distances (l/2 > ls), the stress field range is a constant for a given island geometry, while the interisland distance becomes irrelevant for driven diffusivity in this regime. In this case, the observed suppression of the coupling probability is due to the statistical factor and can be healed by increased thermal diffusivity at a constant stress field gradient.

Nano Lett., Vol. 6, No. 6, 2006

Acknowledgment. This research was supported by the Defense University Research Initiative in Nanotechnology (DURINT) program through DARPA/AFOSR grant no. F49620-01-1-0474. References (1) Xie, Q.; Madhukar, A.; Chen, P.; Kobayashi, N. P. Phys. ReV. Lett. 1995, 75, 2542. (2) Shchukin, V. A.; Ledentsov, N. N.; Bimberg, D. Epitaxy of Nanostructures; Springer-Verlag: Germany, 2003. (3) Guha, S.; Madhukar, A.; Rajkumar, K. C. Appl. Phys. Lett. 1990, 57, 2110. (4) Eaglesham, D. J.; Cerullo, M. Phys. ReV. Lett. 1990, 64, 1943. (5) Nano-Optoelectronics; Grundmann, M., Ed.; Springer-Verlag: Germany, 2002. (6) (a) Tersoff, J.; Teichert, C.; Lagally, M. G. Phys. ReV. Lett. 1996, 76, 1675. (b) Teichert, C.; Lagally, M. G.; Peticolas, L. J.; Bean, J. C.; Tersoff, J. Phys. ReV. B 1996, 53, 16334. (7) Solomon, G. S.; Trezza, J. A.; Marshal, A. F.; Harris, J. S., Jr. Phys. ReV. Lett. 1996, 76, 952. (8) Rouvimov, S.; Liliental-Weber, Z.; Swider, W.; Washburn, J.; Weber, E. R.; Sasaki, A.; Wakahara, A.; Furkawa, Y.; Abe, T.; Noda, S. J. Electron. Mater. 1998, 27, 427. (9) Pal, D.; Towe, E.; Chen, S. Appl. Phys. Lett. 2001, 78, 4133. (10) Rahmati, B.; Jager, W.; Trinkaus, H.; Loo, R.; Vescan, L.; Luth, H. Appl. Phys. A: Solids Surf. 1996, 62, 575. (11) Mateeva, E.; Sutter, P.; Bean, J. C.; Lagally, M. G. Appl. Phys. Lett. 1997, 71, 3233. (12) Schmidt, O. G.; Kienzle, O.; Hao, Y.; Eberl, K.; Ernst, F. Appl. Phys. Lett. 1999, 74, 1272. (13) Le Thanh, V.; Yam, V.; Boucaud, P.; Fortuna, F.; Ulysse, C.; Bouchier, D.; Vervoort, L.; Lourtioz, J.-M. Phys. ReV. B 1999, 60, 5851. (14) Kienzle, O.; Ernst, F.; Ruhle, M.; Schmidt, O. G.; Eberl, K. Appl. Phys. Lett. 1998, 74, 269. (15) Springholz, G.; Pinczolits, M.; Mayer, P.; Holy, V.; Bauer, G.; Kang, H. H.; Salamanca-Riba, L. Phys. ReV. Lett. 2000, 84, 4669. (16) Makeev, M. A.; Madhukar, A. Phys. ReV. Lett. 2001, 86, 5542. (17) Makeev, M. A.; Yu, W.; Madhukar, A. Phys. ReV. B 2003, 68, 195301. (18) Schwoebel, R. L. J. Appl. Phys. 1969, 40, 614. (19) Makeev, M. A.; Madhukar, A. Phys. ReV. B 2003, 67, 073201. (20) Berrie, C. L.; Liu, B.; Leone, S. R. Appl. Surf. Sci. 2001, 175-176, 69. (21) Cazayous, M.; Groenen, J.; Huntzinger, J. R.; Mlayah, A.; Schmidt, O. G. Phys. ReV. B 2001, 64, 033306. (22) Cazayous, M.; Groenen, J.; Huntzinger, J. R.; Mlayah, A.; Schmidt, O. G. Mater. Sci. Eng., B 2002, 88, 173.

NL0602600

1283