Intermixing in heteroepitaxial islands: fast, self-consistent calculation of the concentration profile minimizing the elastic energy

We present a novel computational method for finding the concentration profile which minimizes the elastic energy stored in heteroepitaxial islands. Based on a suitable combination of continuum elasticity theory and configurational Monte Carlo, we show that such profiles can be readily found by a simple, yet fully self-consistent, iterative procedure. We apply the method to SiGe/Si islands, considering realistic three-dimensional shapes (pyramids, domes and barns), finding strongly non-uniform distributions of Si and Ge atoms, in qualitative agreement with several experiments. Moreover, our simulated selective-etching profiles display, in some cases, a remarkable resemblance to the experimental ones, opening intriguing questions on the interplay between kinetic, entropic and elastic effects.


Introduction
In a simplified view, the Stranski-Krastanov growth mode observed in several semiconductor lattice-mismatched systems [1]- [4] is characterized by the formation of a thin wetting layer of material A covering a substrate B, followed, upon further deposition, by the formation of three-dimensional (3D) islands. 3D structures appear when the energy cost associated with the creation of extra surfaces is compensated by the higher strain relaxation allowed for by the increased degrees of freedom. Islands, however, do not allow for a complete relaxation of the elastic energy, which displays a strongly non-uniform behavior, islands edges remaining strongly loaded with stress [5]. Unsurprisingly then, A/B intermixing within the islands has been frequently reported at typical growth temperatures. A very recent review can be found in [6].
From the quoted set of experiments, and with the help of theory [7,8], it emerged that the actual concentration profile in 3D islands depends on the growth conditions and it is determined by the complex and not yet fully understood cooperation between elastic-energy relaxation, frozen kinetics 2 , and entropy maximization (enthalpy of mixing has been shown to be rather negligible [7,10]). Since the growth-parameter space cannot be sampled at will in experiments, a promising way for understanding the role played by each of such phenomena is surely given by directly relating the experimentally derived concentration profile with theoretical simulations performed under conditions tuned to emphasize the different contributions. It is precisely with this goal in mind that in this work we present a very fast method for determining the concentration profile which minimizes the elastic energy in realistically shaped 3D heteroepitaxial islands for a given average concentration value. Beside a very nice model allowing for an analytical treatment of the problem [11] (limited however, to the description of islands formed at the very first stages of growth), a very direct method for finding the optimal Ge distribution is given by exploiting atomic-scale Monte Carlo (MC) simulations, the system energy in each configuration being computed through classical empirical potentials [12]- [15]. Minimization can be achieved by iterating a procedure consisting of chemical-species random exchanges, energy optimization, and acceptance based on the Boltzmann statistical weight [12]. The drawback of the atomistic MC approach is given by the computational cost. The volume of typical experimental islands can easily exceed 10 5 nm 3 (particularly if alloying is pronounced), making an atomistic treatment impossible. If it is true that the self-similar behavior of the deformation field allows one to work at a reduced scale [16], two main problems might be encountered. Firstly, the effective number of bulk-like atoms could become comparable with the surface ones, making surface effects too strong. Secondly, the complex, multifaceted geometry of realistic islands like domes [17,18] or barns [19] works against excessive size shrinking since each facet must be created with a physically meaningful extension in order to avoid artifacts (such as the overlap between the strain field coming from edges and/or vertex). In practice, this means working with islands of at least 10 4 -10 5 atoms, making the convergence procedure extremely long. Motivated by this limitation, we have designed a rather simple and very fast way for achieving elastic-energy minimization within a finite element method (FEM) framework. While in the present paper, we shall provide an application to Ge/Si(001), a system for which several experiments analyzing alloying have been reported [6,10], [20]- [24], it is important to 3 emphasize that the physics involved is rather general (see the interesting discussion in [25]), and applies also, for example, to the popular III-V systems [26]. In InAs/GaAs films, moreover, the larger lattice mismatch is expected to provide a stronger driving force for non-homogeneous alloying, so that all the effects here below discussed should be even stronger. Future work will be dedicated to this system.

Methodology
FEM calculations leading to elastic-energy profiles in 3D Ge/Si (or, in other heteroepitaxial systems) islands can be easily carried out [5] by simply setting an initial strain condition on the island completely analogous to the one obtained in atomistic calculations where the initial geometry is built by positioning the atoms in the diamond lattice determined by the Si lattice parameter. Pure Ge is compressively strained by ∼4% in the three perpendicular directions (in-plane x and y, and out-of-plane z). The initial condition thus reads x x = yy = zz = −0.04 (we follow the convention of describing compression by negative strain values), where x x , yy and zz represent the three diagonal components of the strain tensor. As in atomistic simulations the final elastic-energy distribution can be obtained by allowing the system to relax, so it is done by FEM. At the end of the calculation, the continuum elastic body satisfies the null-stress condition σ · n = 0, where σ is the stress tensor and n is the normal to any free surface. An extension to uniformly alloyed Ge c Si 1−c /Si (0 c 1) systems is straightforward, since it simply requires using an elastic constant linearly interpolated between the Ge and Si values, and imposing the initial strain condition x x = yy = zz = −0.04 × c. In order to treat non-uniform distributions, we propose a particular procedure which has the appealing feature of exploiting only typical operations performed by any FEM code, allowing for a simple implementation. Here, we recall that the core of a FEM calculation consists of dividing the space into elements (typically, tetrahedrons for 3D problems), and by evaluating any function f (x, y, z) considering an interpolation based solely on the values of the function at the nodes of the element containing (x, y, z) (see, e.g. [27]). Here, we shall make use of two different meshes at the same time. A first one, called the concentration mesh, and made of n nodes located at (x i , y i , z i ), i = 1, . . . , n, is used to define a non-uniform concentration profile. Values are assigned at the nodes (c i = c(x i , y i , z i ), i = 1, . . . , n), while a linear interpolation scheme is used for evaluating the concentration at any other point, including the nodes (X j , Y j , Z j ) ( j = 1, . . . , N , N n) of a second, much finer grid (called simply the FEM grid in the following), used for finding the elastic body equilibrium condition. For the systems considered in this work, concentration meshes with n ∼ 400 turned out to be sufficient to accurately describe the concentration profiles in 3D islands (in practice, only 1/8 of the n points were used, see section 3).
Once the initial strain condition is assigned, the minimization of the elastic energy proceeds by an MC iterative procedure. Two vertices of the concentration grid are chosen randomly and their assigned values are randomly changed so that the overall composition of the island remains the same. This constraint can be rather easily fulfilled for linear interpolation leading only to a pre-calculation of necessary weights. The new elastic energy is then calculated by a full FEM calculation, using linear elastic theory. While to this goal we used the commercial Comsol Multiphysics code, the method can obviously be implemented with any FEM solver. Moves (i.e. concentration exchanges) are accepted if the energy is lowered and rejected otherwise. This procedure is repeated until the concentration profile does not change anymore. It is important to emphasize that at each MC step a new solution of the elastic problem is found, Figure 1. Evolution of the concentration profile for a 2D ridge, with average Ge concentration c = 0.6. After ∼3000 MC steps, the energy of the system, plotted in the bottom right panel in the form of ratio between the actual energy and that of the uniform distribution, reaches convergence.
so that the calculation is fully self-consistent: changes in the elastic field accompanying Si/Ge redistribution are correctly described. In order to speed-up convergence, we used a guiding function to bias moves. The bias function is built on the fly, and assigns different statistical weights at the concentration values using information on all previously accepted moves.
A schematic representation of the system evolution during our MC procedure is reported in figure 1 for a 2D case, for which illustration of the system evolution is particularly simple. In the reported example, the Ge island (of height to base aspect ratio 0.2) is actually a ridge, infinitely extended in 1D. The Si substrate is not reported in this figure and in the following ones. The concentration grid was made of 109 points, and the average c was set to 0.6. Experimental Ge and Si lattice constants were used [28]. In the figure, a color map is used to highlight different Ge contents. After ∼1000 steps the system starts converging toward the final situation, where Ge accumulates at the top of the island, whereas Si tends to fill the lower edge region. Notice that the final distribution in figure 1 is nicely symmetric, although no constraints were imposed. Indeed, while first validating our approach, we found out that checking for artificial lack of symmetry was sometimes more indicative than simply looking at the elastic energy convergence. After several trial runs (both in 2D and 3D) we were convinced that our implementation was correct and we started exploiting symmetry to speed up convergence.

Results: realistic 3D islands
Let us now consider more realistic, 3D applications. We have applied our MC-FEM procedure to the three main shapes observed in Ge/Si(001) molecular beam epitaxy experiments [19], i.e. (105) pyramids (height to base aspect ratio (A.R. = 0.1), multifaceted domes (A.R. ∼ 0.2) and multifaceted barns (A.R. ∼ 0.3). The geometry of barns and domes is described in detail in [29]. Ge and Si elastic anisotropy was taken into account in the calculations, even if we verified that qualitatively similar results are also obtained in the isotropic case. Pyramids, domes and barns can be easily divided into 8 equivalent subregions. Since the 3D FEM solution of the elastic problem is computationally rather demanding, considering that it is required at each MC step, exploiting symmetry is rather important for accelerating convergence. We therefore defined our concentration mesh only in the subregions, propagating every attempted concentration exchange in the rest of the island before calling the FEM solver. In this way, the single elastic-energy calculation is not faster, but the number of MC iterations is reduced since fewer concentration mesh points are needed to reach a satisfactory description of the atomic distribution in the subregions. A further advantage is given by the fact that the solution displays the correct symmetry at each step, so that the portion of the sampled concentration space is reduced.
Let us start the analysis of our results by looking at the behavior displayed by (105) pyramids. Simulations were performed using an average Ge concentration value c = 0.75, in order to reproduce the experimental estimate reported in [19]. In figure 2(a) and (b), we report the final Ge concentration in the island using a 3D perspective view ( figure 2(a)) showing c(x, y, z) at the island facets, and a cross-sectional 2D view, obtained by cutting the island along the [010]-direction ( figure 2(b)). Color maps are used to distinguish different concentrations in the online version. A few representative values are also superimposed in figure 2(b), in order to facilitate the reader of a black and white printed version. Overall, there is a spread in concentration of ∼0.35 between Ge-rich and Ge-poor regions, the latter being located close to the island bottom edges. In figure 2(c) and (d), instead, we plot only the subregions occupied by less than 65% (c) and 80% (d) of Ge atoms. The analysis of the corresponding isocompositional surfaces is interesting for two main reasons. If from one side they give a good idea of the 3D Ge distribution in the island interior, the 0.65 case also allows for a direct comparison with the 6 selective-etching experiments of [19], where Ge-rich (above 0.65) regions are eliminated by a suitable etchant. From figure 2(c), it is evident that the concentration profile minimizing the elastic energy in a (105) pyramid of average c = 0.75 displays almost entirely regions where the Ge content is above 0.65. Material is therefore predicted to resist the etching action only at the very bottom of the island, and along the base edges. Comparing our results with the ones of [19], it is evident that the agreement is only qualitative: in both cases it is confirmed that Si-rich regions are located far from the top of the islands, but in the experiments the four corners appear more rich in Si with respect to the center of the edges, i.e. the reverse of the behavior found in figure 2(c). Notice that, as stressed in section 1, differences between the concentration profile minimizing the elastic energy and the real one must in general be expected, since entropic and kinetic effects are not taken into account. With this respect, we find the comparison between our results and the experimental ones particularly useful in evidencing that the actual profiles do not allow for the maximum possible strain relaxation. A further observation can be made. The main difference between our results and the experimental ones does not seem to stem from different degrees of deviation from a uniform distribution, but only from the actual location of the Si-richer regions. In this respect, we believe that the observed difference should not be attributed to entropic effects. Kinetic effects, as suggested in [19] (see also [30]), could explain the observed distribution. The 0.80 isocompositional surface (figure 2(d)), although not corresponding directly to experimental etching profiles, provides further interesting information on Ge distribution, showing that a Ge-rich region, starting from the top of the island (where pure Ge is found), extends in the center of the island reaching the substrate. A similar behavior is found also for dome islands, for which experimental evidence of this hole was reported. Results for domes are reported in figure 3. In this case, always trying to reproduce experimentally reported values [19,31], we considered an average Ge content c = 0.6. Figures 3(a)-(c) display the same kind of results shown in figure 2 for pyramids. Notice that in this case, differences in concentration exceeding 0.9 are reported between the top and the base of the island. Evidently, the strong compression at the edges, induced by the higher aspect ratio [5], strongly attracts Si atoms which are better accommodated in a low latticeparameter environment. The simulated etching profile (figures 3(c)) confirms this observation, showing that removing material exceeding 0.65 of Ge induces a clear hole in the center of the island. Looking at the composition map obtained by atomistic MC at T = 8 K in [12], one can infer an analogous behavior, thus providing a rather convincing proof of the reliability of our continuum-theory-based approach to non-uniform concentration distribution in heteroepitaxial islands.
The comparison between our results for domes and the corresponding experimental results of [19] shows an agreement much more pronounced with respect to the pyramidal-island case, once it is taken into account that convolution effects produce a rounding effect in the images of [19], and that in our figures proportions in the x-, yand z-directions are maintained. Not only is the presence of the central hole confirmed experimentally, but also the slight tendency toward maximal accumulation along the [110]-directions seems to be in agreement. Experiments in [19] were repeated by raising the growth temperature, therefore enhancing the Si concentration in the island, in order to check possible variations in the etching profile. The higher Si content, indeed, was shown to qualitatively change the morphology of the material resisting the etching action which no longer displayed the central depleted zone. In order to check for this change in behavior, we have run a further simulation for c = 0.4, obtaining the result displayed in For the sake of completeness, we have also investigated the behavior of barn-shaped islands, simulated under the exact same conditions (c = 0.6 and 0.4) considered for domes. From our results, it is quite evident that by raising the island aspect ratio, the tendency toward segregating Si at the island base gets stronger. However, the concentration profile now looks more regular, as a direct consequence of the different geometry: isoconcentration lines are more rounded with respect to the dome case, as appears clearly from both the 3D (figure 4(a)) and the cross-cut images ( figure 4(b)). Also the etching profile (figure 4(c)) looks different, the c = 0.6 result for the barn resembling the c = 0.4 case for the dome (depression, not hole). Finally, for c = 0.4, we predict etching to cause a simple and moderate flattening of the very top of the island.
So far we have focused our attention on the concentration profiles within the islands. Let us now look at the influence of the non-uniform concentration profiles on the elastic energy distribution. Since we have verified that qualitatively similar behavior is displayed by all three geometries, we shall directly illustrate only the dome case. In figure 5, elastic-energy maps are presented for c = 0.6 in the case of a uniform distribution ( figure 5(a)) and for the distribution obtained after MC optimization (corresponding to the concentration distribution displayed in figure 3(a)-(c). Energy is expressed with the usual (for atomistic approaches) meV atom −1 ,   figure 3(b). Energies are reported in meV atom −1 . For the uniform case, the elastic energy at the island edges strongly exceeds 10 meV atom −1 . In order to show the main differences between the two cases, however, it was more convenient to display color maps between 0 and 10 meV atom −1 , representing all points above 10 meV atom −1 with a uniform black color.
meaning that a local atomic distribution of atoms yielding, on average, the concentration c(x, y, z) would have the corresponding elastic energy (per atom). Differences between the uniform and the optimized distributions are rather strong, pointing out the importance of using a self-consistent approach. In the former case (already discussed in several papers, for instance 9 see [5]), the upper portion of the island is almost fully relaxed, whereas at the edges the elastic load is particularly significant. Indeed, upon MC-FEM optimization, it is precisely the strong accumulation at the edges which is relieved. While, obviously, the total energy of the system is decreased (by ∼12%) in passing from the situation illustrated in figure 5(a) to the one of figure 5(b), elastic-energy relaxation at the edges is accompanied by a mild (but easily detectable) increase in the central portion of the island, decaying while moving towards its top. Clearly, the common picture of a base-to-top elastic energy relaxation is modified by the non-uniform concentration profile.

Discussion and outlook
When we decided to get involved in the problem of concentration distributions in heteroepitaxial islands our initial goal was to start tackling the problem by providing a fast method 3 able to search for the profile c min = c(x, y, z) minimizing the elastic energy. Being a well-defined quantity, independent of the growth conditions and representative of the maximum possible level of strain relaxation (for a given island shape and average composition), we found c min a quantity of particular interest. We were aware, however, that entropic contributions at experimental growth temperatures were demonstrated by atomistic simulations [12] to be fully dominant over elastic effects, driving perfect intermixing. Notice that also the very low energy values (much lower than k B T for any T relevant for the experiments) reported in figure 5 indicate the negligible stability of the 0 K c min distribution at finite temperatures. In addition, we were also aware of the limitations of a thermodynamic approach allowing for site exchanges for systems where bulk diffusion is frozen. Therefore, we expected to use c min as a simple starting point to analyze deviations at finite temperatures. The successful comparison between the simulated and the experimental selective-etching profiles, came out unexpectedly. Interestingly, it looks like elastic energy is (at least roughly) minimized while entropy is not allowed to rise. These results appear in contradiction with the conclusions of [10]. There, the authors estimated the entropy values based on the experimentally determined concentration profiles. The mixing entropy S, multiplied by the growth temperature T and estimated on the basis of the simple relation holding for random alloys, turned out to be much larger than the elastic energy E el , leading to the suggestion that strain release is not a significant driving force. This result puzzled us, since their experimental profiles were qualitatively rather similar to ours: a dominating entropic term should have prevented the formation of the observed highly non-homogeneous Si distribution in the island. To clarify this point, we have applied their same T S estimate to our minimum-energy concentration profiles, also finding that the T S contribution outweighs the elastic one. For instance, for a 60% Ge-rich dome, the average T S at the typical growth temperature of 873 K is almost 7 times larger than E el , a value very similar to what was reported in [10]. Entropy, however, played no role at all in determining our concentration profile, the only driving force allowed by our MC-FEM procedure being elastic-energy minimization. Obviously, we are not 3 For an island composed of ∼10 4 atoms, we directly verified that atomistic MC based on typical empirical potentials is ∼50 times slower than MC-FEM in providing the minimum-energy concentration distribution. The computational boost allowed for by MC-FEM is however in many cases much higher. Realistic islands at experimental concentration values are in fact made of more than 10 5 atoms, making the use of atomistic simulations prohibitive. Also if the method is applied at finite temperatures, an application not discussed in this work, the speedup guaranteed by MC-FEM becomes much larger, due to the need of averaging over multiple independent atomistic simulations.
saying that entropy plays no role in the experimental profiles, but simply that E el minimization is indeed important, the comparison between the absolute values of T S and E el being misleading in this sense.
It still remains to be explained how the system is able to approach the minimum-energy configuration given that bulk diffusion is frozen. A model capable of predicting the whole formation process based on an accurate description of the system kinetics seems to be needed. The most promising attempt so far is probably the one reported in [8], even if important ingredients, such as the kinetics of Si/Ge exchanges at islands facets, are still missing.
Finally, we would like to point out to the reader an interesting work that appeared while this manuscript was being prepared for submission. In [32], indeed, an alternative approach aimed at yielding concentration profiles within a continuum framework was presented. Optimal distributions within heteroepitaxial islands were investigated in the presence of an extra term, tuned to mimic the additional driving force for perfect alloying (entropic contribution). While the presence of such a term and some approximations (isotropic materials, same elastic constants for island and substrate, simplified island geometry) do not allow for a direct comparison with the present results, it looks like the methodology presented in [32] also provides a powerful tool to investigate non-uniform distributions avoiding a direct atomistic treatment.