Impact of metalorganic vapor phase epitaxy growth conditions on compressive strain relaxation in polar III-nitride heterostructures

A novel approach to estimating the critical thicknesses (CTs) of compressively strained III-nitride layers is suggested, based on a quasi-thermodynamic growth model and accounted for the effect of material decomposition during dislocation half-loop formation on the CT value. The approach provides good quantitative agreement with available data on CTs of MOVPE-grown InGaN/GaN and AlGaN/AlN epilayers. The extremely large CTs observed for high Al-content AlGaN alloys grown on bulk AlN substrates may be attributed, in particular, to the dominant AlGaN decomposition mechanism, producing group-III metallic liquid and gaseous nitrogen. The suggested approach may also be helpful for analysis of threading dislocation inclination in compressively strained layers and applicable to studying point defect formation in semiconductors and its dependence on growth conditions.


Introduction
III-nitride alloys are materials with high lattice constant mismatch between binary constituents, GaN, InN, and AlN. In the case of epitaxy on a binary substrate, the mismatch results in plastic stress relaxation via dislocation formation, limiting the crystal quality of the relaxed epitaxial layers. This limitation is very critical for III-nitride optoelectronic devices, leading to considerable reduction of their emission efficiency because of enhanced non-radiative electron and hole recombination at the generated dislocations. Therefore, knowledge of the critical thicknesses (CTs) of various epitaxial layers used in device structures, and understanding of the factors affecting the thicknesses is extremely important for fabrication technology of the devices.
In conventional zinc-blende III-V semiconductors, the principal relaxation mechanism is gliding of misfit dislocations (MDs). 1) The 60°dislocations glide in the (111) plane of the crystals normally grown on the (001) surfaces. Therefore, there always exists a non-zero component of the shear stress inducing the dislocation gliding. The CT of strained epitaxial layers, corresponding to the onset of plastic relaxation, can be accurately predicted for cubic III-V compounds by the Matthews-Blakeslee (MB) model. 2) That is not the case in polar III-nitride heterostructures, however, where dominant a-type basal-plane dislocations tend to glide in the basal (0001) plane. Here, the resolved shear stress vanishes but the length of the gliding path tends to infinity, resulting in a finite work of mismatch stress during MD formation. In spite of this fact, application of the MB model accounting for anisotropic elastic properties of III-nitrides 3) to polar InGaN/GaN and AlGaN/AlN layers provides CTs one to two orders of magnitude smaller than are obtained experimentally. Possible reasons for this are a rather low resolved shear stress, which may be insufficient for overcoming the Peierls barrier in the basal plane, and a too large gliding path, on which the MDs may be pinned by defects and impurities, suppressing the gliding. Therefore, other mechanisms of stress relaxation may become dominant in the polar structures.
Recently, an alternative mechanism of plastic stress relaxation in (0001) InGaN/GaN layers has been suggested invoking generation of V-shaped dislocation half-loops (DHLs) on the growth surface followed by their penetration to the lower interface of the strained epilayer. 4) The importance of the mechanism was supported by numerous observations of the DHLs-see, in particular, the review of experiments reported in Ref. 4. The mechanical model of CT extended to the case of DHLs predicted CTs much more accurately than the MB model. Nevertheless, the experimental CTs were still noticeably larger than the corresponding theoretical predictions (see Fig. 2

from Ref. 4).
A well-established mechanism of stress relaxation in (0001) AlGaN/AlN layers is inclination of threading dislocations (TDs). 5,6) This mechanism is dominant, however, at a high density of TDs. If AlGaN is grown on a bulk AlN with extremely low TD density (TDD), typically of 10 3 -10 4 cm −2 , inclination of existing TDs becomes insufficient for effective stress relaxation and DHLs are observed in such epilayers 7) similarly to as in the case of InGaN/GaN. However, our attempts to predict CTs for DHLs in AlGaN/AlN layers using the model of Ref. 4 failed to reproduce the available data: 8,9) the experimental CTs were much larger than the theoretical ones. This fact highlights the insufficient understanding of the relaxation processes occurring in AlGaN grown on bulk AlN.
In polar III-nitride structures, TDs forming "arms" of DHLs are inclined with respect to the hexagonal axis of the crystal by a certain angle α. Their specific feature is that the Burgers vector of the dislocations is normal to the inclination plane, i.e. is normal to both dislocation lines. 10) This means that such dislocations are formed by climbing rather than by gliding. Therefore, their formation in compressively strained layers should be accompanied by release of the material confined between the dislocation "arms". The release occurs by decomposition of the material (see Fig. 1). So, it is naturally expected that the CTs of such epilayers become dependent not only on the mechanical (elastic) properties of dislocations and mismatch strain but also on the thermodynamic properties of the solid alloys (InGaN or AlGaN) and the products of their decomposition during dislocation formation.
The goals of this study are: (i) to assess the importance of the thermodynamic aspect for stress relaxation in InGaN/ GaN and AlGaN/AlN layers grown by metalorganic vapor phase epitaxy (MOVPE) and (ii) to provide a more accurate estimation of CTs by accounting for this aspect.

Criterion for dislocation formation
A dislocation may be generated in an epitaxial layer by gliding, if the work of the lattice-mismatch stress W exceeds the elastic energy of the dislocation E dis . For dislocations formed by climbing, which is accompanied by either material release (in compressively strained layers) or material consumption (in the layers under tensile strain), the criteria of dislocation formation should include additionally the chemical energy E ch related to the difference in chemical potentials of the released/consumed material in the final and initial states. The latter is valid if the kinetics of the dislocation formation is not the limiting stage of the process, i.e. if the dislocations are formed faster, for example, than growth of the epitaxial layer occurs and if the nucleation of DHLs on the growth surface does not limit their overall rate of formation. In this case, plastic relaxation of the mismatch stress is expected when , n and S are the normal to the plane of the DHL and the loop area, respectively, ŝ is the mismatch stress tensor, and b is the Burgers vector of the dislocation.
Here, we consider a-type edge dislocations with the Burgers vector = á ñ b 1120 1 3 that are formed in compressively strained polar epilayers. We will distinguish between two types of dislocation: (i) MDs with length L much longer than the thickness h of an epitaxial layer and (ii) DHLs with the dislocation "arms" inclined to the hexagonal axis of the wurtzite crystal by the angle α (at a   90 both types become equivalent to each other). In the former case, the elastic energy of dislocation can be approximated by the expression 2,11,12) q where K(θ) is the energy factor depending on the angle θ between the dislocation line and the Burgers vector 12) and R 0 is the dislocation core radius. In the case of DHLs, the energy E dis should be computed numerically by integrating the elastic fields of infinitesimal prismatic dislocation loops, known in analytical form, over the total DHL area. 4) As an example of the computations carried out within isotropicmedia approximations, Fig. 2 shows two-dimensional distributions of the mismatch stress (s yy ) around a DHL formed in an In 0.1 Ga 0.9 N/GaN layer. One can see that a DHL generates a tensile (positive) stress inside the area confined by the half-loop "arms" and a compressive (negative) stress outside the area. The stress field far from the DHL apex can be considered as a superposition of the fields produced by the individual edge dislocations forming the DHL "arms". Such a consideration fails, however, near the apex where the elastic fields of the "arms" merge substantially. It has also been checked that the value of E dis calculated for DHLs in the case of a   90 corresponds to that predicted by Eq. (1) for MDs.  We assume that material release during dislocation formation under compressive strain occurs via decomposition of the strained material followed by a fast removal of the decomposition products from the area confined by the dislocation lines. If Δμ is the difference in the chemical potentials of the decomposition products and the initial solid phase, corresponding to 1 mol of the material, then the chemical energy E ch can be calculated by the expression: where Ω s is the molar volume of the solid phase.
For dislocations considered here, the normal n and the Burgers vector b are parallel to each other and oppositely directed. Then the resolved mismatch stress s = where M is a specific elastic modulus, h = -/ a a 1 S is the lattice mismatch, and a and a S are the unstrained lattice constants of the epitaxial layer and of the template on which the epilayer is grown, respectively. Using the above expression and assuming that dislocations start to form at one can obtain a general equation For a p < /2 (DHLs), the function F can be found numerically by integrating the elastic energy density over the space. In Eqs. (2)-(3), we used the relationship a = / S L h sin , valid for both types of dislocation under consideration, and assumed that In the case of polar III-nitride structures, 2 33 with C ik being the stiffness constants of epitaxial materials, whereas the energy factor K can be calculated using the method 12) for θ = 90°. In isotropic elasticity theory, where G is the shear modulus and ν is the Poisson ratio.

Chemical potentials
In order to calculate the CT h c from Eq. (2), one needs to evaluate the difference in the chemical potentials of the decomposition products and the initial solid material, which is either InGaN or AlGaN in our study. The chemical potential μ S of the initial solid phase can be straightforwardly calculated. In particular, for In x Ga 1-x N/GaN alloys  (4) is derived for InGaN alloys; in the case of other ternary compounds, e. g. AlGaN, the chemical potential of the solid phase is estimated in the same manner. Evaluation of the chemical potentials of the crystaldecomposition products is a more complicated task. It is closely related to the "thermodynamic reservoir" problem normally introduced for calculating the formation energies of dopants and point defects. 14,15) Here, thermodynamic properties of the reaction products are assumed to be the same as in the equilibrated reservoir with which atoms/molecules are being exchanged. Normally, thermodynamic potentials of species in the reservoir are considered in two limiting cases: the so-called metal-rich and nitrogen-rich conditions. In the former case, the chemical potentials of group-III metals are assumed to be equal to that of the solid metal; in the latter case, the chemical potential of nitrogen corresponds to that of molecular nitrogen (N 2 ) in the gas phase. Both cases are weakly relevant to the real growth conditions used for growth of III-nitride heterostructures by different techniques. In addition, growth of nitride compounds occurs under strong influence of N 2 adsorption/desorption kinetics originating from a high kinetic barrier for N 2 dissociation preceding the adsorption stage. 16) The kinetic effect modifies considerably the properties of N 2 outgoing from the growth surfaces of nitride compounds, which should be accounted for in the choice of an appropriate "thermodynamic reservoir". In any case, such a choice should correlate with a particular growth technique and growth conditions used for fabrication of nitride heterostructures.
Earlier, we developed a quasi-thermodynamic (QT) model of III-nitride compound growth by MOVPE with N 2 adsorption/desorption kinetics taken into consideration. 17) Accounting for atomic conservation on the growth surface, the model relates the flux J k of the k species at a certain point on the wafer with the partial pressure p k of the species at this point via the expression: A k 1 2 are the sticking coefficient and the Hertz-Knudsen factor of the k species, respectively, and N A is Avogadro's number. In this expression, p k 0 is the so-called thermodynamic pressure of the k species. The thermodynamic pressures are related to the chemical potentials of the species via mass-action law equations written for all independent chemical reactions involved in growth. At known p k , the mass-action law equations together with the atomic conservation equations represent a full system for finding all p , k 0 the composition of the growing alloy, and its growth rate. 17) In turn, p k can be found from the species' partial pressures p k in at the reactor inlet, maintained according to the growth recipe by either numerical simulations 18) or an approximate expression derived for a horizontal reactor: 19) = - where B k and D k are the diffusion conductivity and diffusion coefficient of the k species in the gas phase, respectively, H is the reactor height, V F is the mean gas flow velocity over the wafer, and z is the local lateral coordinate counted downstream from the beginning of the growth zone.
Excluding p k from the Hertz-Knudsen representation of the species fluxes and Eq. (5), one can obtain If the sticking coefficient α k is close to unity, then normally a b  B k k k and the expression (6) for the species flux is reduced to @ - The QT-model of MOVPE growth is derived under the assumption that atoms in the adsorption layer are nearly in equilibrium with the growing crystal. Since the a b p k k k 0 product represents the flux of the k species desorbing from the growth surface to the gas phase, it is reasonable to choose p k 0 as the pressure characterizing the thermodynamic properties of the k species in the "thermodynamic reservoir". In this case, the chemical potential of the gaseous species is: 0 where m k 0 is the standard chemical potential of the gaseous k species.

Decomposition routes
Dislocations considered in this study are assumed to nucleate on the growth surface by fluctuations. Therefore, there may be various routes of their formation depending on the particular way of the solid phase decomposition. The most natural way of decomposition is product release to the gas phase by the reverse reactions of those producing crystal growth. In MOVPE of InGaN, these decomposition reactions may be It is an intrinsic feature of the QT-model that the difference in the chemical potentials of the initial crystal and its decomposition product does not depend on the particular choice of decomposition reactions. This fact can be checked by using the above expressions for the chemical potentials of solid and gaseous species along with the massaction law equations relating the thermodynamic pressures p . can be calculated from tabulated chemical potentials of the relevant solid, liquid, and gaseous species. In contrast to the decomposition product release into the gas phase, Δμ is not equal to zero here, as the liquid phase is not in equilibrium with the bulk of growing crystal. A similar consideration can also be applied to AlGaN alloys. In order to estimate the value of Δμ, it is necessary to evaluate p .  At the known partial pressures of species at the reactor inlet, numerical solution of Eq. (9) provides the thermodynamic pressures necessary for calculating Δμ.

Dislocations in InGaN/GaN heterostructures
First of all, the modified model of CT has been applied to MOVPE-grown InGaN/GaN heterostructures. Normally, the InGaN composition is controlled in MOVPE by variation of InGaN growth temperature. In order to find a relationship between the InGaN composition and its growth temperature, we used data from both experimental and theoretical studies on InGaN/GaN growth in the horizontal AIX 200/4 RF-S reactor (H = 3 cm). 20) The growth was carried out at a pressure of 200 mbar in a N 2 atmosphere and an NH 3 molar fraction at the reactor inlet of 0.352 ( = p 7040 NH in 3 Pa). Growth temperature variation from 700°С to 940°С provided a decrease of InN molar fraction in InGaN from x = 0.36 to x = 0.02; the detailed dependence of the InGaN composition on growth temperature was reported earlier in Ref. 20. The gas flow velocity in the reactor V F was about 0.5 m s −1 , weakly depending on the growth temperature.
The diffusion conductivities of the species were calculated analytically using Eq. (5) for the center of a two inch wafer (z = 2.54 cm), using diffusion coefficients estimated from conventional molecular-kinetics theory. 21) Thermodynamic properties of all the species were borrowed from Ref. 13 except for the enthalpy of InN formation, for which the value of 53 kJ mol -1 was chosen, providing good agreement between numerical simulations and experiment. 20) The interaction energy W and its dependence on the InGaN composition was taken from Ref. 22. The sticking coefficient a N 2 of N 2 on the InGaN growth surface was interpolated between those of pure GaN and InN as follows: The sticking coefficients of N 2 at the surfaces of the binary compounds, a N ] 1103 directions. 23) As there is a lack of reliable data on dislocation inclination at other InGaN compositions, we have chosen this value for the whole range of the alloy compositions. Figure 3 shows the CT calculated when accounting for the above Δμ (solid line) and with Δμ = 0 (dash-dotted line). One can see that non-zero Δμ increases the CT, especially at low InGaN compositions where Δμ is high. Moreover, in the case of InGaN decomposition accompanied by liquid phase formation (Δμ ≠ 0), the CT tends to infinity at x ≈ 0.05, i.e. there is no stress relaxation via DHLs at smaller x. Symbols in Fig. 3 show the experimental data on CT collected earlier [24][25][26][27][28][29][30][31][32] and borrowed additionally from Ref. 33. One can see that the theoretical curves calculated for Δμ = 0 and Δμ ≠ 0 represent nearly the lower and upper boundaries for the array of experimental points, respectively. This may be explained by competition between the InGaN decomposition mechanisms, with the product release entirely to the gas phase (Δμ = 0) and that involving liquid phase formation (Δμ ≠ 0). It should be noted, however, that there exists an alternative interpretation for the observed large scatter of the data on CT, which is especially large at low InGaN compositions. In the case of a relatively high density of TDs, one more relaxation mechanism may come into play, reducing the CT: inclination of TDs in the strained epilayer. 5,6) Indirect evidence for this is provided by the data 33) demonstrating dependence of InGaN/GaN CT on the TDD controlled by using either sapphire or bulk GaN substrate for MOVPE growth (the data are shown in Fig. 3 by half-filled pentagons and stars, respectively). Unfortunately, this study did not report detailed investigation into particular relaxation mechanisms of stress relaxation in high-TDD and low-TDD InGaN/GaN epilayers. So, further efforts are needed to identify the dominant routes of stress relaxation in such structures in a wide range of InGaN compositions. Besides, InGaN layers with similar compositions can be grown under significantly different growth conditions, which may also contribute to scatter of experimental CT values.

Dislocations in AlGaN/AlN heterostructures
MOVPE of the layers is normally carried out using H 2 carrier gas. Here, we consider AlGaN/AlN layers grown by MOVPE at 1100°С, a reactor pressure of 100 mbar, and an NH 3 molar fraction at the reactor inlet of 0.01 ( = p 100 NH in 3 Pa). Growth was assumed to be performed in the same AIX 200/4 RF-S reactor with Al x Ga 1-x N composition controlled by adjusting flow rates of metalorganic species.
The diffusion conductivities of the gaseous species were calculated in a manner similar to that discussed in Sect  Little is known about the inclination angle α of DHLs in AlGaN/AlN layers. To our knowledge, there is only one paper 7) reporting on DHL formation in an AlGaN-based structure grown on bulk AlN. In that study, the inclination angle of dislocation lines forming half-loops was found to vary from 20°to 60°. Unfortunately, stress relaxation was initiated in that study by inserting a 10-period Al 0.5 Ga 0.5 N/AlN superlattice between the main Al 0.5 Ga 0.5 N layer and the AlN. Therefore, the data on the inclination angle cannot be regarded as highly reliable. On the other hand, the CT of GaN layers grown on AlN by NH 3 MBE at 900°С was reported to be nearly that of three monolayers, i.e. 0.78 nm, 34) which is close to the value predicted by the MB model. Since the MB model deals with conventional MDs, the above correspondence implies α = 90°for GaN/ AlN layers.
Given the lack of reliable information on the inclination angle of DHLs, we have decided to use instead the data on composition-dependent inclination of TDs. Figure 4 summarizes available data [5][6][7]35) on the inclination angle observed in AlGaN/AlN layers as a function of AlGaN composition. Despite the data scatter, there is a distinct trend of inclination angle decrease with the AlN molar fraction in AlGaN. Approximation of the data with the smooth curve shown in Fig. 4 (see appendix) was then regarded in this study as the composition dependence of inclination angle for new DHLs. Figure 5 compares the CTs of AlGaN/AlN epilayers estimated for the cases of Δμ = 0 and Δμ ≠ 0. In the latter case, calculations were carried out for two growth temperatures: 1100°C (dotted line in Fig. 5) and 1250°С (solid line in Fig. 5). One can see that CT rises dramatically if the DHL formation occurs with group-III atoms releasing to the liquid phase. Moreover, at the same flow rates of metalorganic species at the reactor inlet, the CT depends noticeably on the growth temperature: a lower temperature provides a larger CT. In the case of Δμ ≠ 0, no plastic relaxation via DHL formation is predicted for x > 0.6 at 1250°С and for x > 0.55 at 1100°С. In contrast, the CT predicted for Δμ = 0 varies from 1 nm at x = 0 to 84 nm at x = 0.80 and 1012 nm at There is a lack of reported data on the CTs of AlGaN/AlN layers to compare with our theory. Two points are obtained for x = 0.6 and x = 0.7 by extrapolating the data 8) on relaxation degrees of various AlGaN layers grown on bulk AlN substrates to zero relaxation level (filled diamonds in Fig. 5). These points are close to the CT predicted for Δμ ≠ 0, exceeding that corresponding to Δμ = 0 by about two orders of magnitude. This may be evidence for the dominant route of stress relaxation via DHL formation involving liquid phase for group-III atomic release. The quantitative discrepancy between the experimental points and theoretical curve for Δμ ≠ 0 may be explained by differences in the growth conditions used in Ref. 8 and those assumed in our simulations. There are also data 9) on the growth of relaxation-free AlGaN on bulk AlN, which are shown in Fig. 5 by open squares. Two of the three points lie above the theoretical curve corresponding to Δμ = 0, supporting the above conclusion on the possible mechanism of stress relaxation in AlGaN/AlN.
To our knowledge, there are no data in the literature on the CT of MOVPE-grown Al x Ga 1-x N/AlN layers with small x, including GaN/AlN. In order to understand how stress relaxation occurs in such epilayers, we include the data on the CT of GaN on AlN grown at 900°C by NH 3 MBE 34) (filled pentagon in Fig. 5). According to Ref. 34, CT at this temperature is not affected by the type of substrate used for growth, either AlN/sapphire template or AlN bulk crystal. One can see from Fig. 5 that all the theoretical curves tend to merge at x → 0, becoming weakly dependent on the particular growth technique and being close to the above experimental point.

Growth parameters affecting CT
Our simulations discussed in the previous sections show that growth temperature is the primary factor affecting the CT of grown epilayers in the case of atomic release with group-III liquid formation. In all the cases considered in this study, а  typical Δμ variation with temperature is about 20-30 kJ mol −1 .
The factor of secondary importance is the reactor pressure. Its impact can be illustrated by the case of strained GaN on a relaxed Al 0.21 Ga 0.79 N buffer layer grown by MOVPE in a H 2 atmosphere. Figure 6 shows the CT dependence on the growth temperature calculated for conventional reactor pressure of 200 mbar (the NH 3 molar fraction at the reactor inlet is 0.2) and those for pressures an order of magnitude higher or lower. The typical Δμ variation with the reactor pressure is 5-10 kJ mol −1 in the practically important range of growth temperatures. Therefore, the CT variation with pressure is remarkably smaller than that with temperature. It is interesting that the CT of the GaN layer can be considerably enlarged by lowering the growth temperature from commonly used 1050°C-1100°C to 900°C as typically used for GaN self-doping with carbon. If TDD in the wafer is insufficiently high to provide efficient stress relaxation by TD inclination, the use of low growth temperature enables suppression of new dislocation formation.
The type of carrier gas (H 2 versus N 2 ) is found to influence weakly the CTs of the strained epilayers, since the Δμ variation with the gas composition is less than 2 kJ mol −1 .

Conclusions
As the Burgers vector of DHLs formed in strained InGaN/ GaN or AlGaN/AlN epilayers is normal to the inclination plane, formation of such dislocations should be accompanied by atomic release from the area confined by the dislocation "arms". During the release, solid alloy, InGaN or AlGaN, decomposes, so that the difference of chemical potentials between the decomposition products and the initial alloy contribute to the energy balance, controlling the CT of the epilayers. Since the chemical potential difference depends generally on the growth conditions, CTs of the strained layers become dependent on the growth conditions as well. In order to estimate the impact of the chemical factor on the CT value, we have suggested using the QT-model of MOVPE growth of III-nitride materials, accounting for the kinetic effects of N 2 adsorption/desorption. Considering the effect of the growth conditions on the energetics of dislocation formation, such an approach may also be helpful for analysis of TD inclination in compressively strained epilayers, and applicable to studying point defect generation in semiconductors.
The chemical factor is found to be negligible, if the products of InGaN or AlGaN decomposition are released to the gas phase. In contrast, if the alloys decompose into group-III metallic liquid and gaseous nitrogen, the chemical factor becomes important, generally increasing the CT. In the case of InGaN/GaN layers, the above two decomposition routes provide nearly the lower and upper boundaries for experimental data on the CTs of the alloys, which is explained by competition of the decomposition mechanisms. An alternative explanation invokes inclination of TDs in the InGaN layer as a competitive mechanism of stress relaxation. For AlGaN/AlN, the large CTs predicted for AlGaN decomposition with group-III liquid phase formation and nitrogen release into gas phase agree well with the few experimental data points available for these materials. The dominance of such a decomposition route may explain the fact that AlGaN/ AlN alloys with AlN fraction exceeding 60% can be grown without plastic relaxation up to thicknesses as large as 0.5-1.0 μm.
Our approach to the estimation of CTs of the compressively strained polar epilayers is based on the assumption that the kinetics of the atomic release and nucleation of DHLs at the growth surface are not limiting stages of dislocation formation. Accounting for the kinetic effects in future studies may improve further our understanding of stress relaxation in III-nitride heterostructures.