Curvature and bow of bulk GaN substrates

Gallium Nitride (GaN) is a binary III-V direct band-gap semiconductor material with a Wurtzite crystal structure. Its unique physical properties have made it advantageous for practical applications in optoelectronic and high power/frequency electronic devices.¹ GaN-based device structures have traditionally been grown heteroepitaxially on sapphire or silicon carbide substrates. However, large lattice mismatch between GaN and these foreign substrates yields threading dislocation (TD) densities (TDDs) up to 10⁸–10⁹cm⁻².² Studies have shown that TDs in GaN are electrically active and detrimental to carrier transport in emission properties of semiconductor heterostructures.^3–8 Therefore, there is a strong need to get bulk GaN material for homoepitaxial growth of device structures with low TDD.

Reports in literature show that there have been successful methods to grow bulk GaN by the ammonothermal method, the sodium (Na)-flux method, and the hydride vapor phase epitaxy (HVPE), with reduced TDD (10³–10⁶cm⁻²).^9–14 The ammonothermal method has produced very high quality, low TDD (<10⁴cm⁻²) GaN wafers; however, poor growth rates (<4 μm/h) and purity issues limited the commercial availability of ammonothermal substrates. Imanishi et al. have demonstrated GaN crystals with TDD ∼10³cm⁻² via Na-flux growth;¹⁵ however, challenges with the Na-flux method are primarily the feasibility for upscaling and cost.

The majority of commercially available c-plane bulk GaN is grown by HVPE, with TDD on the order of 10⁴cm⁻² to 10⁶cm⁻². An advantage over other methods is its high growth rate (>100 μm/h), which is approximately two orders of magnitude greater than the ammonothermal and Na-flux method.¹⁶ This achievement has led to suitable GaN substrates for GaN laser diodes (LDs), light emitting diodes (LEDs), and vertical power electronic devices.

However, despite the high quality of HVPE materials, the freestanding substrates have significant curvature. The curvature present in freestanding GaN is caused by internal stresses from within the substrate. The stresses undergo a change in sign and change from the bottom to the top of the substrate, where they change from compressive to tensile stress.

In this report, the TDD and the radius of curvature (R_C) for commercially available bulk GaN substrates were recorded and compared alongside a derived theoretical model for dislocation inclination in freestanding c-plane GaN. From this, a relationship between the stress gradient induced curvature in c-plane GaN to the TDD and threading dislocation inclination (TDI) is demonstrated.

Figure 1(a) illustrates an example of a single crystal substrate with no bow or curvature. This is an example of stress/strain-free material. However, such substrates can also exhibit bowing due to the effects of substrate treatment, as shown in Fig. 1(b). Figure 1(c) is an example of a substrate with curvature due to the change in stress state from compressive to tensile in the material, as previously described. In this case, the curvature causes the substrate to bow. Such crystals can also be cut from the boule to be flat (without bow), but crystal planes remain curved from the stress gradient, as shown in Fig. 1(d).

Radius of curvature measurements via x-ray diffraction (XRD) is a true measure of the crystal strain as opposed to optical techniques (such as a multi-beam optical stress sensor) that use an array of parallel incident laser beams to reflect off the substrate surface and are collected by a detector. The latter method is heavily influenced by the substrate thickness and surface treatment, rather than actual crystal lattice elastic strain. For example, Figs. 1(c) and 1(d) have the same curvature, but only Fig. 1(c) is bowed. Therefore, to accurately measure the R_C of a substrate, regardless of which bowing state is present, XRD is the preferred method.

To accurately analyze the curvature of a bulk crystal substrate, a multi-layered structure model is used as shown in Fig. 2(a). The total thickness H is divided into N layers with arbitrary thickness h_i. The reference length, L, included a fixed number of lattice parameters (with lattice constant, a_ref) with no elastic strain, stresses, or bending. The layers are then separated as shown in Fig. 2(b), such that they are unconstrained and each layer undergoes a lattice transformation, expanding, or contracting with respect to its respective reference length. The number of lattice parameters within L remains the same, introducing transformation strain $εi*$⁠, also known as stress free (eigen) strain,¹⁷ expressed as

where a_i is the lattice parameter after undergoing the lattice transformation.

Coherent layers deform in a compatible manner, meaning a 1:1 correspondence between the lattices at the interface of contacting layers such that during the lattice transformation, the lattice maintains continuity across the interface. The lateral strains in the layers will have a z-dependence, with the coordinate origin at the bottom of the first layer. The compatible deformation of the multilayer structure with the respective coordinates is shown in Fig. 2(c), and the required conditions of the compatibility of total strain ε^T are shown as¹⁸ 

The general expression for total strain can then be expressed as

where A and B are unknown coefficients used to derive an expression for curvature. Given that each layer is subjected to a different magnitude of self-strain, the compatible total deformation mentioned above, ε^T, is maintained by the layer's elastic strain layers, ε_i, written as

The elastic strain in the ith layer will have a z-dependence as the compatible total deformation and self-strain have a z-dependence; therefore, ε_i can be rewritten as

The elastic bending shown in Fig. 2(d) accompanies non-uniform elastic strain. The inherent strain in the individual layers of the structure causes the multilayered structure's curvature, thus, to develop an expression for curvature, an analysis of the stress state is required. The z-dependent elastic strain is related to mechanical stress by assuming Hooke's Law and plane stress state in layers as shown

where σ_i(z) and M_i(z) are the z-dependent mechanical stress and the biaxial elastic modulus of the ith layer, respectively. The two unknown coefficients, A and B, can then be determined by using mechanical equilibrium conditions for force and bending moment to yield

where I_k and J_k are defined as

This is done in order to simplify the expressions for the coefficients. Assuming that the self-strain and biaxial moduli are known for the system, the curvature $κ$ can be found using the relation

which reduces to

The approach described above is consistent with the results presented in the literature.^19–21 

This report attributes the compressive to tensile stress gradient in freestanding GaN to curvature caused by threading dislocations²² and their inclination within the material. Cantu et al. reported relaxation of compressively strained Si-doped (0001) Al_0.49Ga_0.51N deposited on Al_0.62Ga_0.38N, where TEM discovered threading dislocations with initial [0001] lines directions in Al_0.62Ga_0.38N inclining toward the $⟨11¯00⟩$ directions within Si doped Al_0.49Ga_0.51N.²³ 

These threading dislocations have an effective misfit dislocation (MD) component that generates stress, which may be in the same sense or different senses to other stresses in the layer. Furthermore, a model describing the conditions for dislocation inclination and relaxation due to “effective climb” of edge dislocations by Romanov and Speck²⁴ explains the mechanism behind the relaxation of compresses stresses in Al_xGa_1−xN. Further growth with “frozen-in” inclined dislocations is shown to generate a tensile strain gradient as demonstrated by Cantu et al.,²⁵ which is responsible for cracking at a critical thickness.²⁶ Accord et al. and Manning et al. further demonstrated in situ growth wafer curvature measurements revealing compression to tension stress transition in compressed Al_xGa_1−xN.^27,28

The mechanism by which stresses relax/generate and form gradients are related to the inclination of the TD with respect to their original direction. The most common method of observed relaxation in lattice-mismatched semiconductor heterostructures is misfit dislocation (MD) formation at the interfaces, perpetuated by dislocation glide. This gliding process can proceed by bending the pre-existing threading dislocations or by generation of new dislocations. There is, however, a glide-free mechanism of elastic stress relaxation related to the inclination of the TD with respect to their original direction and their effective climb as shown by Romanov et al. (2003), where the total biaxial far-field plastic relaxation at the top layer surface resulting from the triangular MD grid ε^l_pl is given by

where b is the Burgers vector, ρ_TD is the threading dislocation density, and L is the projected length directly related to the layer thickness h and the inclination angle α, such that L = h tan(α). The inclination angle, α, is formed when the straight dislocation originating from the substrate bends at the substrate/layer interface and becomes an angular dislocation as shown in Fig. 3. The lines of these dislocations are inclined with respect to the [0001] growth direction, with α as large as 20°. The TD inclines towards the ⟨$11¯00$⟩ direction. From this, the average plastic relaxation for a layer is then

Research grade GaN crystals from the commercial vendors Nanowin, Furukawa, Mitsubishi Chemical Company (MCC), and Hitachi Cable, designated as samples A, B, C, and D, respectively, were acquired. All substrates are c-plane (0001). The substrate thicknesses were 300 μm, 630 μm, 420 μm, and 400 μm, respectively.

The radius of curvature was measured via XRD using a Rigaku Smartlab high-resolution x-ray diffractometer. The Rigaku diffractometer is a horizontal θ-θ diffractometer, which facilitates wafer placement on the stage to be free of constraints that otherwise may impact the curvature. An open detector in double axis configuration with a Ge (220) 4 bounce monochromator was used to create a parallel incident beam. A 3 × 5 mesh with 3.5 mm spacing was mapped to measure the R_C as shown in Fig. 4. Note that the direction of the scan is along the x direction which is coplanar with the scattering plane that is defined by the incident and scattered wave vector.

The R_C was determined by performing on-axis (0002) omega scans, then translating the sample and measuring the displacement of the omega peak positions. The difference in the incident beam angle from two different regions exposed to x-rays gave the R_C. It is expressed as²⁹ 

where x₁ and x₂ are respective positions along the x-direction and ω₁ and ω₂ are the omega peak positions at those positions.

The threading dislocation density was determined using cathodoluminescence (CL) imaging. CL was performed using an FEI Scanning electron microscope (SEM) with a Gatan MonoCL4 system to obtain panchromatic cathodoluminescence images as well as CL spectra. Images were taken with beam excitation energies in the range of 5–10 kV, corresponding to electron penetration depths in GaN of approximately 150–500 nm as calculated using the method described in Ref. 30.

Equation (13) can be derived through geometric analysis of the general experimental setup and implementation of basic diffractive principles (Fig. 5). Modeling the substrate surface as the circumference of a circle with radius, R, a geometric relationship can be established between the two scan points, 1 and 2, which are Δx apart on the lab horizontal frame. The angle of the incident beam with the lab horizontal is ω, the angle of the incident beam with the crystal vertical is χ, the angle between the lab and crystal horizontals is Φ, and the angle between the crystal vertical and lab horizontal is δ.

Since the same peak in reciprocal space is being observed, it can be assumed that the Bragg's condition is identically satisfied at both x positions, that is, χ₁ = χ₂ in Fig. 5 for all positions on the surface when observing a single peak in reciprocal space. The relation between Bragg peaks can be expressed as the angle, β, between the scan positions relative to the center of the R_C

The assumption is then made that the arc length, γ = Rβ, is equal to the horizontal shift of the sample stage. This assumption holds true as long as the deviation of the arc from the chord defined by Δx, σ_r, is negligible. Mathematically, this condition would be satisfied if Δx ≫ σ_r or Δx ≪ R. Both are reasonable assumptions considering that radii of curvature are generally on the order of 10 m–1000 m, while translations are on the order of 5 mm. This makes the following:

Our initial analysis of this system for a convex system was done with the lab frame using the same Cartesian unit vectors as the R_C frame. However, upon shifting to a concave derivation, the R_C frame had to be rotated about the y-axis to maintain a positive R vector. This shift changes the sign on the Δx term as the x-axis unit vector rotates about the y-axis. This can now be simplified to

as previously described by Fewster and thus derived regardless of the radial frame of reference used.

The radii of curvature for the commercial bulk GaN crystals are displayed in Table I.

Substrate A had the largest R_C with the largest standard error, 8.08 ± 0.97 m. Substrate B and substrate C had similar radii of curvature, 5.42 ± 0.23 m and 6.44 ± 0.16 m, respectively. The smallest R_C was recorded from substrate D with 2.75 ± 0.17 m.

As mentioned previously, R_C was measured through a series of omega scans displaced along the x direction, parallel to the source and detector of the diffractometer. The radii of curvature of the bulk GaN crystals further indicate that the substrates were bowed concave up as shown in Fig. 6. With the substrate bowed concave up, the incident beam would reach the Bragg condition at reduced angles as it traversed the curved surface. This would, in turn, result in a Bragg peak at a smaller ω for one-dimensional translations where Δx > 0, hence a negative slope in Equation (13) is realized. For simplicity, the R_C was represented as an absolute value.

Panchromatic CL images for the commercial bulk GaN crystals along with the corresponding evaluated TD density are displayed in Fig. 7 and Table II, respectively. The CL images show predominantly bright areas corresponding to band gap luminescence from defect free GaN as well as dark spots corresponding to non-radiative regions around defects that have been shown to be associated with threading dislocations, as in Rosner et al. The CL images show a fairly uniform concentration of TDs along the region of the imaged surface, independent of length scale, while the corresponding SEM images of the bulk crystals all exhibited planar smooth surfaces (not shown).

Substrate A, shown in Fig. 7(a), had the lowest number of TD per scanned area, with a TDD of 4.16 × 10⁵cm⁻². Substrate B shown in Fig. 7(b) had a TDD of 1.35 × 10⁶cm⁻², with a low standard error compared with the rest of the sample set. The TDs in substrate B are distributed in small groups or clusters. Substrate C shown in Fig. 7(c) had a TDD of 1.26 × 10⁶cm⁻²; however, this value is associated with a larger standard error, due to large variations in the distribution of dislocations at larger scales (not shown). Substrate D shown in Fig. 7(d) had the highest density of threading dislocations, at 2.92 × 10⁶cm⁻².

To model the dislocation inclination in freestanding GaN, a relationship between curvature, dislocation density, and threading dislocation inclination angle must be established. The curvature of a bulk multi-layered structure was derived to be equal to the unknown coefficient B in Eq. (10), which takes into account the thickness dependence of the self strain $εi*$ (z), defined as the lattice transformation that each layer undergoes as it expands or contracts with respect to its reference length. The self strain subjected to each layer of the structure is equivalent to the average plastic relaxation resulting from the glide free mechanism of stress relaxation related to the inclination of TD from Romanov et al., as shown below

Equation (17) can then be used to solve for the unknown coefficient B, so that the curvature can be determined. This result derives an expression where the curvature can be proportional to the TDI and TDD. This statement assumes that the biaxial modulus no longer has z-dependence as $εi*$ is assumed to be the z-dependent plastic relaxation. After appropriate calculations are made, the curvature is

The thickness dependence is eliminated, and the curvature developed by the average plastic relaxation is directly related to the TDD and the average inclination angle of the TDs. Figure 8 shows a three dimensional schematic for inclined TD in wafers grown in [0001] direction.

Dwilinski et al. reported on the growth of exceptionally high quality bulk GaN substrates using the ammonothermal method.³¹ They reported TDD of ∼5 × 10³cm⁻². Regardless of the crystal size, the R_C was reported to be greater than or equal to 1000 m. Figure 9(a) plots the derived model for varying inclination angles alongside the ammonothermal bulk GaN crystal grown by Dwilinski et al. The substrate supported the model, with an estimated inclination angle of ∼15°.

Figure 9(b) shows the dependence of R_C with respect to TDD for the GaN substrates investigated. The experimental R_C and TDD values are plotted alongside the derived theoretical model over varying inclination angles. All of the commercial bulk GaN crystals investigated correlated well with the theoretical relations derived, supporting the proposed model relating the curvature from stress gradients to TDD and TDI.

According to the derived theoretical relation, substrate A had a TDI angle of 20°, with minimal deviation in both the R_C and TDD. Substrate B had an inclination angle of 10°. The standard error was negligible for the R_C and TDD; therefore, this is the most precise data point of the batch. Substrate C crystal was speculated to have a ∼10° inclination angle, but the standard error is significant along the TDD; therefore, it could have an inclination ranging from 10° to 17°. Substrate D is speculated to have an inclination angle between 8° and 10°. From here, we can conclude that the theoretical model holds for bulk GaN crystals of large and small TDD and R_C.

The ranges of inclination angles (8°–20°) chosen for this study are based of the observed and measured inclination angles from the previous studies.^23,32,33 TEM studies done by Cantu et al. reported a TDD of ∼3 × 10¹⁰cm⁻², where the TD was inclined by 10°–25° in the $11¯00$ direction and had a misfit component, thus relieving misfit strain. Observations of Follstaedt et al. for compressively strained AlGaN indicated that the TD inclination occurred prior to the introduction of Si, where the average inclination angle was ∼19°, and the TDD is ∼6.9 × 10⁹cm⁻². Follstaedt et al. further examined AlGaN films with increasing Ga content, reporting TD inclination angles ranging from 6.7° to 17.8° and TDD of 6–7 × 10⁹ cm⁻².

At the TDD observed in this report for HVPE bulk GaN (10⁵–10⁶cm⁻²), it is difficult to measure TDI via cross-sectional TEM, as the TDD is too low. Another method to measure the TDI is via x-ray topography; however, the TDD observed is too high for this technique. TD can be observed at this density when using Electron Channeling Contrast Imaging as performed by Picard et al. and Carnevale et al.^34,35

The average radius of curvature for the commercial HVPE grown GaN substrates is approximately ∼5 m, leading to a large disparity in vicinality from edge to edge over the substrate (e.g., ∼0.6° for R_C = 5 m for a 50 mm substrate). For comparison, the vicinality of GaN on sapphire is optimized for ∼0.2° miscut in the m-direction of GaN which corresponds to 0.2° towards a-direction sapphire.

Therefore, even for GaN substrates suitable for LD and LEDs (TDD ∼ 10⁴cm⁻²–10⁶cm⁻²), there still exists a large non-negligible curvature. From this, one can see the importance of lowering the TDD in the future so as to achieve the flattest, least bowed substrates possible and avoid the detrimental effects to carrier transport and emission properties of TD have in semiconductor heterostructures as demonstrated by Rosner et al.

The radius of curvature and threading dislocation density for commercially available freestanding GaN substrates were measured. Substrates with higher threading dislocation densities have smaller radii of curvature. Furthermore, a model was derived to investigate the dislocation inclination in GaN substrates and its dependence to TDD and R_C. The model was developed by assuming that the curvature is caused by internal stresses within the substrate that undergoes a change from compressive to tensile stress from the bottom to the top of the substrate. The experimentally determined R_C and TDD for substrates with uniform distribution of TDD confirm and support the derived theoretical model.

This work was supported by the King Abduallah Center for Science and Technology and King Abdullah University of Science and Technology (KACST/KAUST) as well as the Materials Research Laboratory and California Nanosystems Institute at UC Santa Barbara for providing access and training to their laboratories. A.E.R. got the support for the conducting theoretical part of this study from Russian Science Foundation Project No. 14-29-00086.