Kinetic Monte Carlo Simulations of Quantum Dot Self-Assembly Kinetic Monte Carlo Simulations of Quantum Dot Self-Assembly

In the Stranski–Krastanov growth mode for heteroepitaxial systems, layer-by-layer growth is followed by the formation and growth of three-dimensional (3D) islands. In this paper, we use a kinetic Monte Carlo method to simulate this growth mode behavior. We present a detailed and systematic investigation into the effects of key model parameters including strain, growth temperature, and deposition rate on this phenomenon. We show that increasing the strain lowers the apparent critical thickness that is defined by the onset of 3D island formation. Similarly, increasing the growth temperature lowers the apparent critical thickness, until intermixing, and the resulting relevance of entropic contributions, become more significant. We also report the impact on Stranski– Krastanov growth of more model-specific parameters, such as bond strengths between constituent atoms of the system, and surface energy anisotropies.


Introduction
In this paper we discuss results from a model that uses kinetic Monte Carlo (KMC) simulations to describe the growth of heteroepitaxial systems. Under strain, heteroepitaxy often proceeds via the so-called Stranski-Krastanov (SK) mechanism, where two-dimensional (2D) layer-by-layer growth transitions into three-dimensional (3D) island growth mode after a critical thickness is reached. Strain in these systems results from the difference in lattice constant between the substrate and the material being deposited. Archetypal examples of SK growth in semiconductor systems include the growth of InAs on GaAs [1], and Ge on Si [2,3]. The resulting 3D islands are often referred to as quantum dots (QDs).
QD formation is the result of an instability, and this has been investigated theoretically, from the continuum and mesoscales, to the microscopic (atomic) scale. The basic driving force for this instability is a competition between minimizing the surface energy of the system, which stabilizes layer-by-layer growth, and the strain relief that occurs by forming 3D features [4,5]. The early work of Asaro and Tiller [6] and Grinfeld [7] predicted an instability that forms immediately with a critical thickness of zero, in contrast to the experimentally observed SK growth mode. Spencer, Voorhees, and Davis [8] suggested that while the instability is indeed present immediately, it is not observed until 1 Matthew Abramson and Hunter J. Coleman contributed equally to this work an apparent critical thickness is reached. To explain why SK growth begins with several flat layers, Tersoff [9] introduced the concept of a wetting potential, which one can interpret as a thickness-dependent surface energy of the deposited material.
In a 1+1 dimensional KMC model, Baskaran and Smereka [10] showed that all of these concepts are relevant to describing SK growth properly, but that intermixing also plays a crucial role. Intermixing has two effects: it reduces the effective lattice misfit, and it enhances the role of entropic contributions, which become more important as the temperature increases. A follow-up paper by Schulze and Smereka [11] confirmed that the results observed in Ref. [10] for 1+1 dimensions also hold for a 2+1 dimensional KMC model. However, due to computational limitations, the work of Schulze and Smereka only showed results for a few model parameters. By exploring the effects of all model parameters, in this paper we seek to provide a comprehensive study of this KMC model.
Our renewed interest in this topic has been motivated in part by a some unexpected results in recent experiments, particularly in tensilestrained QD systems. The first concerns the thickness of the 2D wetting layer during SK growth. Typically, once the critical thickness is reached and QDs begin to form, the 2D wetting layer stops growing. layer continues to grow even after the 2D-to −3D SK transition has occurred [12,13]. The second concerns the fact that for Ge QDs grown on InAlAs(111)A, researchers have observed a transition from SK growth to the Volmer-Weber (VW) growth mode as the substrate temperature is increased [14]. The KMC model we describe here will form the basis for future studies aimed at understanding these phenomena.
The organization of this paper is as follows. We begin by describing our 2+1 dimensional KMC model, before showing detailed results of the effect of strain, temperature, and deposition rate for a set of fixed bond-strengths that eliminate any effects of surface energies. We then vary the relative bond strengths between substrate atoms and deposited atoms, which amounts to including the effects of different surface energies. Some key results, such as the formation of QDs upon the introduction of strain, are only weakly dependent on the surface energies.

The model
The KMC model that we use has been described in detail in Ref. [11]. It is a cube-on-cube, bond-counting model. The state of the system is described by a discrete height array ℎ , supplemented by a discrete displacement field, . The total energy of the system consists of the bond energies and a discrete elastic energy , and = + . Transitions from one state to the next occur with rates that depend only on the initial state and the location of the hopping particle: where { } indicates the initial position of the atom making the move and −( 0 + ) is the energy barrier for the transition. We take 0 as a fixed constant, and is the difference between the energies with and without a surface atom at ( , ). Like the total energy, consists of two pieces: one that describes the bonding found in bond-counting schemes without elastic effects and that depends only on ℎ , and one that depends only on the elastic energy, As described in Ref. [15], we shall consider two species of atoms, denoted type 1 and type 2. For most of our simulations we will consider the situation in which atoms of type 2 are deposited onto a substrate of type 1. For this model, for a surface atom at site ( , ) is given by where is the strength of the interaction (bond strength) between atoms of type and type , and (1) , denotes the total number of bonds of type and connecting the atom at site ( , ) and its nearest neighbors. (2) and (3) are defined analogously, but for next-nearest neighbors and next-to-next-nearest neighbors, respectively. We choose This implies that is the energy barrier for the diffusion of a type 2 adatom on a type 1 substrate (ignoring the term for now). The parameters , , and allow one to vary the anisotropy of the crystal. For example, the surface energy per unit area for the (001) facet of material 1 is Finally, is the total elastic contribution to the energy. We account for the elastic interactions using a ball-and-spring model with longitudinal and diagonal springs having spring constants and respectively. The elastic effects arise because the natural bond lengths of materials 1 and 2 are different. The elastic energy of the system is then of the general form where is the respective relative change in spring length. In analogy to , is given by the difference in elastic contributions with and without a surface atom at position ( , ). These elastic energies are obtained by minimizing the elastic energies with respect to all . In later sections we will use the surface roughness , which we define as where ℎ is the average height of the system, and the angular brackets ⟨⟩ signify the spatial average over all lattice sites.

Results and discussion
Here, we provide a systematic study of how the key parameters in our KMC model affect the formation and evolution of quantum dots. In 3.1 we investigate the effect of strain for a system where the bond energies between any two atoms of type 1 or type 2 are all the same. We then explore the effects of temperature 3.2 and deposition rate 3.3, which can usually be controlled experimentally, for example during QD growth by molecular beam epitaxy. In 3.4 we investigate in greater detail the role of different bond energies, which translates into the role of different surface energies for atoms of type 1 and 2. Finally, in 3.5 we examine how the choice of different surface energies for different facets affects the resulting shape and orientation of the QDs.
All results shown in Section 3.1 and Section 3.2 are for = 0.3, = 0.5, = 1.0, bond energies 11 = 12 = 22 = 0.2425 eV, and a deposition rate of 1.0 ML/s (monolayer per second). We used an identical value for the three bond energies in these sections to isolate the effect of different surface energies from the effects of other parameters. We chose this specific value since it is the average of the bond energies 11 = 0.26 eV, 12 = 0.2425 eV, and 22 = 0.225 eV used in Ref. [11], and it corresponds to a typical value for semiconductor materials. Our choice for , , and leads to surface energies 100 = 2007 erg/cm 2 , 110 = 1712 erg/cm 2 , and 111 = 1637 erg/cm 2 . We note that many conclusions discussed below are qualitatively similar given different choices for , , and , and for .

Effect of strain
In Fig. 1, we show how the surface morphology of 4 ML of type 2 atoms deposited on type 1 atoms at = 700 K evolves as the strain is increased from 0 % to 7 %. For strains < 4 %, we observe layer-by-layer growth for these parameters and no QD formation. When the strain reaches 4%, we enter the SK growth mode. Following the growth of a 2D wetting layer of ∼ 2 − 3 ML thickness, a transition to 3D growth occurs and QDs begin to form. We refer to the thickness of the wetting layer at this moment as the apparent critical thickness, consistent with the work of Baskaran and Smereka [10]. In Fig. 2 we plot the time evolution of surface roughness for all simulations shown in Fig. 1. During layer-by-layer growth the roughness oscillates  Surface roughness as a function of coverage for systems with 0%, 2%, 4%, 5%, 6%, and 7% strain at 700 K. The simulations correspond to the ones shown in Fig. 1. The transition from layer-by-layer growth to Stranski-Krastanov growth is indicated by an arrow. between 0 and 0.5, with each period corresponding to the formation and completion of a new atomic layer. These roughness oscillations in our model are analogous to the intensity oscillations of the reflection high-energy electron diffraction (RHEED) signal routinely seen during molecular beam epitaxy [16]. In both cases, one oscillation period corresponds to the growth of a complete monolayer of the crystal. There is a clear signature in the roughness plots (indicated by a dashed red line) that indicates the SK transition from the roughness oscillations corresponding to layer-by-layer 2D growth, to a monotonic increase in due to 3D QD formation. For a given set of simulation parameters, we refer to the deposition amount required to produce this transition as the apparent critical coverage. We note that the formal definition of the apparent critical thickness is slightly different in Ref. [10].
As the strain increases, the apparent critical thickness decreases so that the onset of QD formation takes place at lower deposition (Fig. 3). The formation of QDs lowers the elastic energy, while at the same time it increases the overall surface energy. This strain relief mechanism becomes more important as the strain increases, and therefore occurs at an earlier time [17][18][19]. We also see that as the strain increases, the island radius decreases, the island height increases, and the island areal density is increased. This is also well known from experimental results, for example the early work of Snyder et al. for In Ga 1− As on GaAs(100) [1]. For 7% strain the QDs start forming after the deposition of approximately 0.6 ML. This indicates that the system has transitioned to the VW growth mode, where 3D growth of QDs starts before the completion of the first 2D monolayer.

Effect of temperature
In Fig. 4 we show the surface morphologies after the deposition of 4 ML under 5% strain at different temperatures. The corresponding apparent critical thicknesses are shown in Fig. 5. At = 500 K, we have layer-by-layer growth until ∼ 2.8 ML deposition, and Fig. 4(a) shows the beginning of the QD formation. As the temperature increases, the apparent critical thickness decreases. The explanation is that there is an energetic driving force for the formation of QDs to relieve strain. At low temperatures, QD self-assembly is kinetically limited, but as the temperature increases, we get closer to equilibrium, and we see the formation of QDs at lower coverage [14,20]. This result also agrees with the work of Baskaran and Smereka [10]. For temperatures above = 600 K, the apparent critical thickness stabilizes around 1.3 ML, and does not decrease any further. In fact, our data suggests that the apparent critical thickness increases slightly between = 650 K and = 700 K. At first glance, this is an unexpected result. But an analysis of the intermixing of atoms of type 1 (the substrate atoms) and atoms of type 2 (the deposited atoms) can explain this behavior. At temperatures below = 600 K, there is essentially no intermixing of atoms of type 1 and type 2. All atoms in layer 0 (which we define as the top atomic layer of the original substrate) and below are of type 1. All atoms in layer 1 (i.e., the first atomic layer above the original substrate) and above are of type 2. This can be seen in Fig. 6, where we show the concentration of atoms of type 2 in layer 0. At low temperatures, it is essentially 0, which means that almost all atoms in this layer are of type 1. As the temperature increases above = 600 K, this intermixing fraction increases, and it is around 0.2 and 0.4 at = 650 K and = 700 K, respectively. This increased intermixing is also apparent upon visual inspection of Fig. 4, where we clearly see more type 1 atoms exposed at the top layer as the temperature increases. The observed stabilization of the apparent critical thickness of the wetting layer happens because increased intermixing leads to a reduced effective lattice mismatch. Furthermore, intermixing increases the relevance of entropic contributions that stabilize the thickness of the wetting layer.

Effect of deposition rate
For homoepitaxial growth, increasing the temperature produces many trends that are similar to what one observes when reducing the atomic deposition rate. The reason is that the adatom diffusion length, and correspondingly the density of islands (in the submonolayer regime), depend on the ratio ∕ , which can be increased or decreased by changing either the diffusion constant or the deposition rate . In fact, in the absence of any additional processes such as edge diffusion or adatom detachment from island edges, increasing or lowering leads to identical morphologies. However, once adatom detachment is allowed as an additional process, this equality no longer holds. Let us assume for simplicity that only singly bonded atoms detach from an island edge at a rate , which is dependent on the temperature , but is independent of the deposition rate . One can now change and to keep ∕ constant, but at the same time, ∕ cannot be kept constant.  During homoepitaxial growth, the differences in morphology that emerge upon varying vs. are small for a wide range of parameters, In contrast, for heteroepitaxial growth, these morphological differences are more significant and lead to qualitative changes that result in QD formation. Fig. 7 shows the morphologies of 4 ML under 6% strain deposited at four different rates. Raising the deposition rate leads to a higher density of smaller QDs due to a reduction in adatom diffusion length. This outcome is similar to what we observe when we reduce the growth temperature (see Section 3.2). This result is also in agreement with results obtained from continuum models for strained island growth by Aqua and co-workers [21,22].
However, when we look at the effect on the apparent critical thickness (c.f. Fig. 8), we see a significant qualitative difference. Above we have shown that a lower temperature leads to prolonged layer-bylayer growth, and a larger apparent critical thickness. In Fig. 8 we see the opposite effect, where raising the deposition rate beyond 10 ML/s results in the system immediately forming QDs. The reason is that lowering the temperature reduces ∕ but lowers ∕ even more, and so reduces the relative importance of detachment. Atoms are less likely to detach from an island edge and move to a higher atomic layer, thus stabilizing and prolonging layer-by-layer growth. In contrast, when we increase , the relative importance of ∕ vs.
∕ remains constant. Furthermore, since a higher value of leads to smaller islands, it is now easier for atoms that detach from an island edge to hop up to a higher layer, destabilizing the layer-by-layer growth.

Effect of surface energies
In this section we discuss the effect on surface morphology when we vary the bond energies 11 , 12 , and 22 . We present results for the cases where we keep 12 fixed, and increase (decrease) 11 while decreasing (increasing) 22 by the same amount. Table 1 summarizes the values used for 11 , 12 , and 22 for the five cases studied here. We acknowledge that this discussion is by no means complete. One could, for example, also choose to keep 11 fixed, and then vary 12 and 22 . Fig. 9 shows surface morphologies after the deposition of 4 ML under 5% strain for cases (a)-(e) in Table 1. In the preceding sections we showed results for 11 = 22 , which is the case shown in Fig. 9(c). As we increase 11 (panels (b) and (a)), the most noticeable effect is   intermixing as defined earlier is shown for all five cases. For cases (b) and (a) there is essentially no intermixing. With increasing 11 , the substrate gets stabilized, the apparent critical thickness decreases, and QDs form a little sooner (c.f. Fig. 10). If we instead decrease 11 , we see from Figs. 9(d)-(e) that intermixing is enhanced. The substrate gets de-stabilized, and type 1 atoms float to the top. This can also be seen from the larger intermixing fraction shown in Fig. 11. The increased intermixing in Figs. 9(d)-(e) also lowers the effective misfit, and for the parameters discussed here, it appears that the surface begins to roughen very early without the clear formation of QDs. We note that the intermixing fraction in Fig. 9(e) is slightly smaller than in Fig. 9(d), and not larger. This happens because we increase the bond strength between two atoms of type 2 when we lower the bond strength between atoms of type 1, and once the bond between atoms of type 2 is too large, the intermixing is reduced. We also show the measured apparent critical thickness for Fig. 9(d)-(e) in Fig. 10. In Section 3.2 we showed that enhanced intermixing at higher temperature can stabilize layer-by-layer growth due to the increased importance of entropic contributions. In this section we have demonstrated that the detailed interplay between different atomic interactions is also an important factor influencing intermixing and the associated stability of the thin film.

Effect of facet anisotropies
As described in Section 2, the parameters , , and define the anisotropy of the surface energy for different facets. In this section we show how our model can capture important details such as surface anisotropies. For all results discussed so far, we have chosen = 0.3, = 0.5, and = 1.0, which leads to surface energies 100 = 2007 erg/cm 2 , 110 = 1712 erg/cm 2 , and 111 = 1637 erg/cm 2 . The (110) and (111) facets are the preferred facets (with a ratio 110 ∕ 111 = 1.05). Indeed in all results shown above, these were the dominant side facets for the QDs, while the (100) facet was least favorable. The (100) facet appears only at the top of each QD, and disappears as the QDs get taller (see for example Fig. 1). We note that in many figures we see more (110) than (111) facets, which we presume is because of kinetic limitations.
If we now set the next-next-nearest neighbor interactions to zero (i.e., = 0.0), these surface energies change to 100 = 733 erg/cm 2 , facet now becomes energetically more favorable than the (110) facet. More importantly, the ratio 110 ∕ 111 = 1.13, so that the (111) facet is now even more preferred energetically compared to the (110) facet. In Fig. 12 we show results of the surface morphologies with = 1.0 and with = 0.0. The main effect is indeed that for = 0.0 the (111) facet starts to dominate. In addition, since the (100) facet has now a lower surface energy than the (110) facet, the (110) side facets are being replaced by (100) (i.e., vertical) side facets, and the tops of the QDs do not get smaller as they grow. This suggests that if one were able experimentally to control the different surface energies (for example using surfactants), one might be able to manipulate the orientation and details of QD shape.

Application to III-V growth
As mentioned in the introduction, our renewed interest in this model has been motivated in part by some recent experiments on III-V QD self-assembly [12][13][14]. We would like to discuss briefly why a cubic KMC model is relevant to such systems, and how one could include the specific effects of the presence of the group V species. In almost all cases, molecular beam epitaxy of III-V materials takes place under group V-rich conditions, and growth is limited by the kinetics of the group III species [23]. An excess of the group V species is provided to stabilize the III-V surface at the growth temperature, and so whenever a group III atom attaches to an island or step edge, there is a group V atom available to rapidly bond to it. We can therefore think of the cubic atom used in our model as the group III element. We can capture the details of the bonding of the group III atom to the system (and to the group V atoms) by effective parameters for the cubic atom in the model. These effective parameters can be tuned to a specific system, for example by matching RHEED patterns [24].
The other important effect of the presence of group V atoms is their impact on the surface morphology. For example, a combined theoretical and experimental study of InAs(001) showed that the stable surface reconstruction changes from 2(2 × 4) to 2(2 × 4) as the As overpressure is reduced [25]. The mobility of adatoms and the surface diffusion anisotropy is determined by the surface reconstruction. One must therefore modify the parameters of the KMC model presented here to account for the different reconstructions that result as we vary group V overpressure.

Conclusions
In this paper we have systematically studied the effects of strain, temperature and deposition rate on the apparent critical thickness of the wetting layer and formation of QDs during SK growth. These are parameters that can be experimentally controlled. In general, increasing the temperature lowers the apparent critical thickness, while also increasing intermixing between atoms that are deposited and atoms in Fig. 13. Apparent critical thickness for 6.5% strain showing the Stranski-Krastanov transition from layer-by-layer growth to 3D QD formation as a function of temperature. The apparent critical thickness drops below 1.0 ML for ≥ 500 K, indicating a transition from SK growth at lower temperatures to VW growth. the substrate. In a certain intermediate regime, we find that raising the temperature increases intermixing just enough that increased entropic contributions stabilize the wetting layer, and therefore reverse the trend of a decreasing apparent critical thickness. The deeper understanding of intermixing and entropic contributions we have gained during this study will play a key role in future investigations of the growth of GaAs QDs on InAlAs(111)A, where experiments unexpectedly show that the wetting layer continues to grow even after QD formation has begun [12,13].
Of the model parameters that we have investigated, strain i.e., the lattice mismatch between the substrate and the deposited material, is the most influential. As strain increases, we find that the apparent critical thickness decreases, leading to a transition in the growth mode from layer-by-layer, to SK, and ultimately VW growth. The preliminary results in Fig. 13 suggest that for a strain of 6.5% we can find a transition from SK growth (with an apparent critical thickness above 1.0 ML) to VW growth (with an apparent critical thickness below 1.0 ML). These results contrast with what we saw for 5 % strain (Fig. 5) where the apparent critical thickness also decreased with temperature, but remained above 1.0 ML, indicating SK growth even at the highest temperatures we considered. A more detailed study of the effect of the bond energies and intermixing on the SK-to-VW transition, and possible relations of our chosen model parameters to real III-V materials, will be subject of future work. We therefore believe that this model is well suited to exploring experimental observations of this SK-to-VW growth mode transition in the Ge on InAlAs(111)A QD system [14]. Finally, we have investigated how QD formation is affected by other model parameters such as bond strength and bond anisotropy, that are more difficult to control experimentally.

Declaration of competing interest
The authors declare that they have no known competing financial interests or personal relationships that could have appeared to influence the work reported in this paper.

Data availability
Data will be made available on request.