Strain Relaxation of Si/SiGe Heterostructures by a Geometric Monte Carlo Approach

Solid‐state‐based quantum technologies, such as electronic spin‐qubits, constitute a leading approach to the realization of quantum computation. Electronic spin‐qubits hosted in semiconductor heterostructures demand the highest crystalline quality, specifically with respect to the structure and formation of in‐plane misfit dislocations (MDs). Here, the formation mechanisms of MD networks in such Si/SiGe heterostructures are investigated. For this purpose, strained Si layers on relaxed Si0.7Ge0.3 with thicknesses above the critical thickness are grown by molecular beam epitaxy to allow for the formation of MDs. MDs are mapped out by electron channeling contrast imaging and the experimental results are compared to MD networks calculated by Monte Carlo simulations. The simulations show that the ratio between the threading dislocation density (TDD) and the average length of MDs is a constant for a given set of parameters, here 2.9 for sufficiently large dislocation statistics and simultaneous extension of MDs. Further, the fractal dimension of the MD network determined by the box‐counting method as an alternative quantity to characterize MD networks is introduced. This allows to infer the formation mechanism of MDs in real devices and hence the quantification of resulting average relaxation within the Si‐quantum well.


Introduction
For quantum computers, the limited fidelity of real qubits imposes an overhead of error corrections requiring in a network of thousands, or even millions, individual qubits. [1] In this context, gate-defined quantum dots in SiGe heterostructures are one of the most promising systems for the large-scale integration of spin qubits since they take advantage of state-of-the-art fabrication processes of existing Si-based transistor technology. [2][3][4] Furthermore, for Si and Ge, it is possible to create nuclear spin-free heterostructures by isotopic enrichment. [5] However, the fabrication of SiGe heterostructures with sufficiently low defect densities and lateral homogeneity is a remaining challenge in material science. [6] In the SiGe material system, the necessary quantum well to confine electrons in the vertical direction can be established using strained Si grown on Si 1Àx Ge x with a typical value of Ge concentration of x ¼ 0.3. Here, the tensile biaxial strain lifts part of the conduction band degeneracy, i.e., valley degree of freedom, which ultimately allows probing of spin quantum states of individual electrons with sufficiently long spin coherence times, a requirement for quantum computing. [7] Growing of such SiGe single crystals cannot be achieved by bulk crystal growth methods such as the Czochralski technique. [8] Hence, Si 0.7 Ge 0.3 buffers must be epitaxially grown on Si. [9] However, the inherent lattice mismatch between Si and SiGe leads to a range of defects, such as compositional fluctuations and dislocations. [10,11] These ultimately lead to strain relaxation in quantum well layers above a certain critical thickness, leading to local changes in the band splitting, and ultimately decoherence of the spin states. [12] Hence, controlling, understanding, and reducing these defects, in particular, the relaxation due to in-plane misfit dislocations (MD), is compulsory for developing a functional SiGe-based quantum computer.
Here, we show that the comparison of purely geometric Monte Carlo simulations with experimental observations of MDs provides insight into the kinetics of formation and formation mechanism of MDs.

Formation of MDs in Si/SiGe
The Si quantum well to host the electronic spin qubits is grown on SiGe virtual substrates. The required SiGe virtual substrates are epitaxially grown on Si (001) and inherently suffer the presence of threading dislocations (TDs) with a density between 10 5 and 10 7 cm À2 . [13,14] Moreover, they exhibit compositional and sinusoidal thickness undulations aligned along 110 directions, so-called cross-hatch patterns. [10] DOI: 10.1002/pssr.202200398 Solid-state-based quantum technologies, such as electronic spin-qubits, constitute a leading approach to the realization of quantum computation. Electronic spin-qubits hosted in semiconductor heterostructures demand the highest crystalline quality, specifically with respect to the structure and formation of in-plane misfit dislocations (MDs). Here, the formation mechanisms of MD networks in such Si/SiGe heterostructures are investigated. For this purpose, strained Si layers on relaxed Si 0.7 Ge 0.3 with thicknesses above the critical thickness are grown by molecular beam epitaxy to allow for the formation of MDs. MDs are mapped out by electron channeling contrast imaging and the experimental results are compared to MD networks calculated by Monte Carlo simulations. The simulations show that the ratio between the threading dislocation density (TDD) and the average length of MDs is a constant for a given set of parameters, here 2.9 for sufficiently large dislocation statistics and simultaneous extension of MDs. Further, the fractal dimension of the MD network determined by the box-counting method as an alternative quantity to characterize MD networks is introduced. This allows to infer the formation mechanism of MDs in real devices and hence the quantification of resulting average relaxation within the Si-quantum well.
The TDs will extend under the presence of a mismatchinduced strain of the tensile strained Si quantum well and form MDs at the interface between the quantum well and the SiGe virtual substrate, once the thickness of the Si-layer exceeds the critical thickness as described by the Matthews-Blakeslee criterion, as is depicted in Figure 1b. [15] The resulting formation of MD segments leads to partial relaxation of the Si layer. This relaxation is laterally inhomogeneous and must be avoided for strained quantum well layers.
According to the theory by Matthews and Blakeslee, the elongation of TDs to form MD segments is a result of a balance between a driving force originating from the mismatch-induced strain -the Peach-Koehler force acting on the TD and line tension of the forming MD. While this line tension is independent of the dislocation length, the Peach-Koehler force increases linearly as the TD length increases with the layer thickness. Hence, above a certain critical thickness h crit an MD is forming, elongating the dislocation until either the end of the crystal is reached, or the MD segment elongation is blocked by an obstacle. Such an obstacle may be, for instance, MD lying on its path. The blocking of MD elongation is related to another critical thickness, as elaborated elsewhere, and associated with the Freund criterion. [16] The thicknesses of Si quantum well layers on SiGe are typically around 10 nm and consequently above the critical thickness by Matthews and Blakeslee and below the critical thickness according to Freund. [14] Hence, one observes the formation of MDs by elongations of TDs and the blocking of MDs by MDs. With the MDs being oriented along the [110] and ½110 directions─the two 110 directions in the (001) plane lead to a characteristic appearance of the MD networks in strained Si. The quantification of MDs with respect to pre-existing TDs, as well as the formulation of statistical criteria to characterize distribution patterns of MDs, constitutes the focus of this work.

Basic Definitions Related to MD Networks
Before we describe dislocation networks in SiGe heterostructures, we introduce some basic definitions. The threading dislocation density (TDD) is defined as the total threading dislocation length (TDL) divided by the crystal volume V which contains them, yielding TDD ¼ TDL V . This is approximately equal to taking a random plane through this volume and diving the number of intersecting dislocations N by the area A of the plane, which yields TDD ¼ N A . Analogously, the misfit dislocation density (MDD) can be defined as the total misfit dislocation length (MDL) divided by the area of the heterointerface A, yielding MDD ¼ MDL A . The MDD can be approximated by taking a random straight line through the heterointerface and dividing the number of crossing dislocations N by the length of the line L, so that MDD ¼ N L . Throughout this work, we simplify and assume that TDD ¼ N A and MDD ¼ N L , which is generally valid for large enough statistics. Based on these definitions, the average spacing between dislocations can be calculated. The average threading dislocation spacing (TDS) can be written as and the misfit dislocation spacing (MDS) can be described as

Relationship Between TDD and Relaxation
Here we consider the relaxation of the strained Si layer by MD formation. Essential for describing the relaxation is the average MDS and hence, the MDD. More specifically, the misfit accommodated by dislocation can be described by simple geometric dislocation configurations and Equation (2) as where b is the Burgers vector of Si. [15] The average MDS can not be independent of the TDD, and it is natural to ask what is the intrinsic relationship between MDS and TDD, and consequently between relaxation and TDD. The details of this relationship depend on the specific dislocation network, but some general statements can be formulated based on simple assumptions as demonstrated below.

Growth
Strained 28 Si layers were grown by molecular beam epitaxy (MBE) on Si 0.7 Ge 0.3 virtual substrates produced by the www.advancedsciencenews.com www.pss-rapid.com Lawrence Semiconductor Research Lab. The influence of interface contamination between Si and SiGe was avoided by growing a 30 nm Si 0.7 Ge 0.3 on the virtual substrate prior to Si growth. The Si was grown above the critical thickness at 350°C followed by annealing in the MBE chamber at 600°C for 10 min to trigger the formation of MDs and while ensuring the blocking of MDs by MDs.

Electron Channeling Contrast Imaging
Electron Channeling Contrast Imaging (ECCI) experiments were conducted in a Thermo Fisher Scientific's Apreo S scanning electron microscope (SEM) at a beam energy of 10 keV and a current of 3.2 nA. The SEM T1 in-lens backscattered electron detector, which is located inside the final lens, was used to acquire the micrographs at a working distance of 3 mm in electromagnetic immersion mode.
Since the ECCI contrast is only visible in a limited area, several images were stitched together to obtain a larger field of view. Tilt and rotation of the sample holder were utilized to ensure the sample was always oriented under the selected diffraction condition.

Monte Carlo Simulations
Here we apply a geometric Monte Carlo approach in two dimensions to investigate the MD networks in Si/SiGe heterostructures. [14] In this regard, the term "geometric" refers to calculations not involving any real physical interaction between dislocations. More specifically, assume a random distribution of TDs serving as the starting point for MD growth, in general at random times and in one of the four 110 directions in the (001) interface. Furthermore, we assume that MDs are blocked by other MDs, as described by the Freund criterion. Any further dislocation interactions such as by long-range elastic strain fields are neglected.
The simulation itself is scale-invariant and the only intrinsic length scale is given by the TDS and the finite size of the simulated area. In the case of this work, we normalize all length measures to the average TDS according to Equation (1). Normalizing the length scales to the TDS, one can easily show that the average misfit dislocation length (MDL avg ) is equal to the ratio TDS/MDS.
To appreciate the significance of this, we recall an interesting result from a previous investigation, namely the constancy, or statistical robustness, of the ratio TDS/MDS. [14] The MDD is defined as the total length of MDs divided by the interface area. If N is the number of dislocations, one can write the total MDL as the N ⋅ MDL avg and the MDD as The number of dislocations typically is varied between 1000 and 10 000 per simulation. The space of discretized elements consists typically of a 2 12 Â 2 12 grid, which turned out to be a good trade-off between convergence and computational cost.
Boundary conditions have to be considered but turned out to play an inferior role if the number of dislocations is large enough, i.e., TDS is sufficiently small compared to the size of the simulated area. Two natural choices are periodic boundary conditions and non-periodic boundary conditions (dislocations end at the simulation boundary). Both yielded similar dislocation structures and all results shown throughout this work are based on non-periodic boundary conditions.
The factors which strongly influence the MD structure are the initial distribution of TDs and the exact model of how MDs grow, in particular starting time and velocity of MD growth. Natural special cases are simultaneous MD growth, all at the same speed, or subsequent growth with a, i.e., Poisson distribution of growth velocities.

Results
The MD structure of a 30 nm thick strained 28 Si layer grown on Si 0.7 Ge 0.3 was imaged by ECCI as depicted in Figure 2a. The MD segments are visible as fine white lines. The TDs originating from the SiGe substrate are characterized by a short bending of the lines exemplary marked by the black arrow in Figure 2a. The TD segments going through the strained Si layer cannot be resolved by ECCI.
The visibility of the MD segments depends considerably on the sample tilt and is a trade-off between the field of view and the contrast of the MDs. Hence, to investigate a significant portion of the MD network a total of more than 20 ECCI measurements were stitched together with high overlap, as depicted in Figure 2b. The red box marks the ECCI image in Figure 2a. The thereby identifiable MD network is marked by black lines in Figure 2b. The TDs are marked by black triangles in Figure 2b. Every MD segment arises from a TD while the end is blocked by another MD. The long-ranged modulations of the ECCI contrast are caused by the cross-hatch pattern of the SiGe substrates.
The observed TDD in this area of 160 Â 133 μm 2 (minus the two corner areas) is 1.8 Â 10 5 cm À2 , with a TDS of 23.6 μm according to Equation (1). The total length of MD segments in the respective area is 1957 μm and consequently, the MDD is 975 cm À1 . By simply applying Equation (2) one obtains an MDS of 10.3 μm and a TDS/MDS of 2.3. As will be discussed below, the MDS is strongly influenced by MD formation conditions. Comparison with simulations can therefore provide insights into the mechanisms of MD growth.

MD Networks by Monte Carlo Approach
Calculating the MD network arising from a random distribution of TDs requires some assumptions, in particular when MD segments form over time. Not surprisingly, we find that the distribution in start-of-growth times strongly influences the final MD structure. This opens an opportunity to investigate the nature of MD formation in real samples. Hence, several models for MD formation are applied in the simulations and compared to experimental results.
First, two extreme choices are investigated, namely all MDs growing simultaneously with the same velocity, and each MD growing after each other with arbitrarily chosen starting times (but the same growth velocity for all). The MD networks arising from simultaneous and non-simultaneous MD formation for these two scenarios are shown in Figure 3a,b, respectively. A total of 10 4 TDs were randomly placed for these simulations. The MD networks are apparently very different. The network in Figure 3a seems more "homogeneous" and consists of MD segments with a relatively narrow distribution of lengths (see Figure). In contrast, the MD network in Figure 3b seems very "inhomogeneous" displaying some dominant and long MD-segments initiated by the MDs that grow first. Consequently, the statistical distribution of MD-segments, in this case, appears broader.
The average MD length for Figure 3a,b is 1.95 and 2.52 TDS, respectively. Consequently, the ratio TDS/MDS is 1.95 and 2.52 for Figure 3a,b, respectively.

The Ratio Between TDS and MDS
A series of 300 independent simulations with varying random TD distributions has between conducted, where MDs grow either simultaneously or one after another (non-simultaneous). The corresponding TDS/MDS ratios can be seen in Figure 4a as black and red points, respectively. The graph shows the ratio between TDS and MDS of the final MD network for each simulation. In contrast to experimental methods, with our simulations MD networks arising from a very high number of TDs are accessible. This allows us to find the limit of the ratio of TDS/MDS for a very large number of TDs, and furthermore, predict the influence of dislocation number on this ratio. A striking result is the  www.advancedsciencenews.com www.pss-rapid.com constancy of this ratio for sufficiently high TDD, which we already reported elsewhere. [14] However, the exact value of this ratio depends strongly on how the MDs actually grow and provides means to investigate MD formation in real samples. Going towards sparse TDD, the TDS/MDS ratio continuously decreases. This can be understood as a finite size effect of the simulation and one can derive an analytic expression assuming that the MDL of individual dislocations is completely independent of another where ðTDS=MDSÞ ∞ is the ratio of TDS/MDS for infinitely many TDs. Applying this equation to the data in Figure 4a for simultaneous MD formation with a ðTDS=MDSÞ ∞ ¼ 1.95 yields an excellent fit for the complete range of TDs. However, the model does not fit well with the non-simultaneous MD formation, since the TDS/MDS ratio converges much slower than N À 1 2 to ðTDS=MDSÞ ∞ . We found empirically that the same Equation (4) works quite well for an effective scaling with N À 1 3 as shown here which is plotted in a solid red line in Figure 4a. The respective value of is ðTDS=MDSÞ ∞ ¼ 2.9.

The Relaxation of the Quantum Well
Using the result TDS=MDS ¼ const., one can infer the dependency of the average relaxation of the strained Si layer as a function of TDD. First, we can write MDS ¼ TDS const: and together with Equation (3), we can write relaxation ¼ b 2 ⋅ const: TDS , which together with Equation (1) results in Hence, by knowing the value of ðTDS=MDSÞ ∞ -which is constant for MD networks arising from strained Si layer grown below the critical thickness according to Matthews and Blackeslee-one can calculate the misfit relaxation accommodated by dislocations. Using the same simulations as in Figure 4a, the corresponding relaxation was calculated for a series of several hundred simulations for simultaneous and non-simultaneous MD formation with an increasing number of TDs, as is depicted in Figure 4b. The respective Equation (6) is plotted for ðTDS=MDSÞ ∞ ¼ 1.95 and ðTDS=MDSÞ ∞ ¼ 2.9 in black and red, respectively. Here, the simulations ranged from 10 2 to 10 4 TDs and the simulation space was denoted with a physical dimension to reflect reasonable TDDs occurring in SiGe virtual substrate, similarly the relaxation is given in units of %. The relaxation by MDs is with a maximum of 0.02% seemingly low, but still substantial compared to the total strain in Si quantum well layers on Si 0.7 Ge 0.3 of around 1%. Especially, considering the inhomogeneous nature of the MD distribution.

Comparison of the Experimental and the Simulated MD Network
From comparing the Monte Carlo MD networks with the experimental data, one can hope to obtain insight into the nature of the formation of the MD network. Comparing the experimental MD network of Figure 2 with the simulated MD networks in Figure 3, suggests a closer resemblance of experimental data with the MD network arising from non-simultaneous MD formation. This is also suggested by the ratio of TDS/MDS, which is 2.3 for the experimental MD network and 2.52 for the result of the Monte Carlo simulation depicted in Figure 3b.
For a more detailed comparison, we limit the simulations to a comparable number of the experimentally accessible number of  (4) and (5), respectively. (b) shows the respective average relaxation as a function of threading dislocation density (TDD), and the solid lines are a simple model for relaxation as by Equation (6), due to MDS for around TDS/MDS ¼ 2 and TDS/MDS ¼ 3, depicted in black and red, respectively.
www.advancedsciencenews.com www.pss-rapid.com MDs. The influence of the boundary conditions of Monte Carlo simulations was minimized by calculating the MD network according to the Experimental part on a four times larger area and selecting the center square. Figure 5a depicts the determined MD network from experimental data shown in Figure 2b (rotated by 90°) together with a simulation for a comparable number of TDs with simultaneous (Figure 5b), and non-simultaneous MD formation (Figure 5c), respectively. Consistent with the average MDL or TDS/MDS ratio, it seems that the dislocation network of the measured sample in Figure 5a resembles the network from non-simultaneous MD formation in Figure 5c and is rather distinct from the simultaneous MD formation shown in Figure 5b. We would like to define a quantitative measure capturing the apparent differences in MD distributions and to distinguish the MD formation mechanism. One could consider the standard deviation of MDLs and MDSs, but MDS is generally not normally distributed, and instead show log-normal distribution. [17] Hence, the standard deviation is not a good statistical measure and other parameters of statistical distributions have to be considered. Hence, we analyzed several approaches including the possibility to compare the MD networks by their self-similarity by comparing their fractal dimension, as reported elsewhere. [18] Here, we apply the same box-counting method to determine the fractal dimension of the networks. The initial box size is 1= ffiffiffi ffi A p and gets bisected consecutively resulting in a fourfold number of boxes in each subsequent iteration step. The MDL for each box size is determined by summing up all boxes which contain at least one dislocation. This total MDL as a function of box size is shown for the MD networks depicted in Figure 5a-c in Figure 5d-f, respectively. For the initial iterations and as long as the box size is larger than the TDS, one expects a quadratic increase of the MD density, which is the case. For further iterations, the boxes are small enough to capture the self-similarity of the MD network. For determining the fractal dimension of the datasets, the data points displaying quadratic increase were excluded from fitting. The extracted box-dimension for Figure 5d-f is then 1.24, 1.15, and 1.25. Hence, the fractal dimension of the MD network in Si/SiGe is somewhere between the fractal dimension of the coastline of Australia and the Koch curve. [19] The fractal dimension seems to capture the qualitative similarity between the experimental MD network and Figure 5c. This suggests a quantitative way of determining the nature of the MD formation by comparing their fractal dimension.

MD Formation and Cross-Hatch pattern
Finally, we seek to understand why the model of non-simultaneous MD formation fits the experimental MD network best. In principle, this suggests that there exist local strain fluctuations that cause some MDs to grow much earlier than others. Here, we propose that these could arise from the inherent cross-hatch pattern in SiGe heterostructures. For this purpose, we created an artificial www.advancedsciencenews.com www.pss-rapid.com 2D sinoidal undulation in the form of sinð x λ CH Þ þ sinð y λ CH Þ, where x and y are spatial coordinates and λ CH is the undulation period. This artificial cross-hatch pattern is depicted in Figure 6a,d, for λ CH > TDS and λ CH ¼ TDS, respectively.
The respective MD networks for simultaneous and nonsimultaneous MD formation are shown in Figure 6b,c,e,f, respectively. Here, simultaneous means all MDs start growing at the same time, but the speed is proportional to the local value of the cross-hatch modulation. Non-simultaneous in this model means that the TDs form MDs one after each other beginning with the one with the highest local value of the cross-hatch modulation to the last one with the lowest value, consecutively. As before, the speed is proportional to the cross-hatch modulation.
Interestingly, the cross-hatch modulation has basically no influence on the MD network arising from simultaneous MD formation. Furthermore, both MD networks are comparable to the MD network of Figure 3a. In contrast, for non-simultaneous MD formation, the cross-hatch modulation has a superior relevance. For λ CH > TDS, the cross-hatch modulation gets imprinted in the MD network, while for λ CH ¼ TDS this is not the case, and the MD formation seems random, such as in Figure 3b.
Hence, these results motivate the assumption that the driving force for MD formation depends on local strain variations caused by the cross-hatch pattern in SiGe. The effects are most pronounced for λ CH TDS but are still visible for λ CH ¼ TDS. For our Si/SiGe-sample, the average λ CH is 5 À 10 μm and the TDS is % 23.6μm (Figure 2b), consistent with our model.

Conclusion
In this work, the formation of MD networks, and respective accommodation of misfit relaxation in Si/SiGe heterostructures have been investigated by comparing experimental observations of MDs in a real Si/SiGe stack to MD networks calculated by Monte Carlo simulations. The Si-layer was grown coherently strained on strain-relaxed Si 0.7 Ge 0.3 above the critical thickness according to Matthews and Blakeslee. The arising MD network is observed by ECCI. It was possible to show that all MDs follow the [110] and ½110 directions in the (001) heterointerface. MDs propagate until they are blocked by other MDs. The observed spacing of MDs was 10.3 μm with a ratio of TDS/MDS of 2.3. The MD networks calculated by 2D-Monte Carlo simulations were derived by different MD formation rules, either all MDs form simultaneously or MDs form one after each other in a random order (non-simultaneously). The arsing MD networks look qualitatively different, which can be quantified in several key figures such as TDS/MDS ratio and fractal dimension. We showed that the ratio TDS/MDS is a constant with the exact value depending on the formation mechanism. Since the convergence of this ratio with the number of MDs is quite slow, a pure experimental analysis of MD networks would require several hundreds of dislocations. The ratio for simultaneous MD formation is 1.95, while it is 2.9 for non-simultaneous formation. Alternatively, we utilized the self-similar structure of MD networks to characterize them. The fractal dimension determined by the box-counting method is 1.15 and 1.25 for the simultaneous and non-simultaneous MD Figure 6. a,d) Artificial cross-hatch-like modulation with different periodicity utilized λ CH to evaluate the influence on the MD network. (b,e) arose from for simultaneous MD formation, while (c,f ) form non-simultaneous MD formation. The dislocations located in the region with a high cross-hatch value grow first.
www.advancedsciencenews.com www.pss-rapid.com networks, respectively. The experimental results show a fractal dimension of the MD network of 1.24, indicating non-simultaneous MD formation. This is also in agreement with the overall qualitative appearance of the experimental MD network and its TDS/MDS ratio of 2.3. The non-simultaneous MD formation may be caused by the cross-hatch pattern typically present in SiGe-heterostructures and its corresponding local strain modulation. The relaxation of the quantum well layer caused by MDs was quantified and identified as a potentially significant disturbance for electron-based quantum technologies.