Curvature and bow of bulk GaN substrates

We investigate the bow of free standing (0001) oriented hydride vapor phase epitaxy grown GaN substrates and demonstrate that their curvature is consistent with a compressive to tensile stress gradient (bottom to top) present in the substrates. The origin of the stress gradient and the curvature is attributed to the correlated inclination of edge threading dislocation (TD) lines away from the [0001] direction. A model is proposed and a relation is derived for bulk GaN substrate curvature dependence on the inclination angle and the density of TDs. The model is used to analyze the curvature for commercially available GaN substrates as determined by high resolution x-ray diffraction. The results show a close correlation between the experimentally determined parameters and those predicted from theoretical model. V C 2016 Author(s). All article content, except where otherwise noted, is licensed under a Creative Commons Attribution (CC BY) license (http://creativecommons.org/licenses/by/4.0/) . [http://dx.doi.org/10.1063/1.4959073]


I. INTRODUCTION
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. 1 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 8 -10 9 cm À2 . 2 Studies have shown that TDs in GaN are electrically active and detrimental to carrier transport in emission properties of semiconductor heterostructures. [3][4][5][6][7][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 3 -10 6 cm À2 ). [9][10][11][12][13][14] The ammonothermal method has produced very high quality, low TDD (<10 4 cm À2 ) GaN wafers; however, poor growth rates (<4 lm/h) and purity issues limited the commercial availability of ammonothermal substrates. Imanishi et al. have demonstrated GaN crystals with TDD $10 3 cm À2 via Na-flux growth; 15 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 4 cm À2 to 10 6 cm À2 . An advantage over other methods is its high growth rate (>100 lm/h), which is approximately two orders of magnitude greater than the ammonothermal and Na-flux method. 16 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.

II. BACKGROUND
A. Substrate curvature and bow 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.

B. Derivation of curvature in bulk multi-layered structure
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 e Ã i , also known as stress free (eigen) strain, 17 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 e T are shown as 18 @ 2 e T z ð Þ @z 2 ¼ 0: ( 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, e T , is maintained by the layer's elastic strain layers, e 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 zdependence; therefore, e i can be rewritten as The elastic bending shown in Fig. 2(d) accompanies nonuniform 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 r 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 j can be found using the relation j which reduces to j ¼ ÀB: The approach described above is consistent with the results presented in the literature. [19][20][21] C. Inclined threading dislocations and stress gradients in III-nitride structures This report attributes the compressive to tensile stress gradient in freestanding GaN to curvature caused by threading dislocations 22 23 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 24 explains the mechanism behind the relaxation of compresses stresses in Al x Ga 1Àx N. Further growth with "frozen-in" inclined dislocations is shown to generate a tensile strain gradient as demonstrated by Cantu et al., 25 which is responsible for cracking at a critical thickness. 26 Accord et al. and Manning et al. further demonstrated in situ growth wafer curvature measurements revealing compression to tension stress transition in compressed Al x Ga 1Àx N. 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 , where the total biaxial far-field plastic relaxation at the top layer surface resulting from the triangular MD grid e l pl is given by where b is the Burgers vector, q TD is the threading dislocation density, and L is the projected length directly related to the layer thickness h and the inclination angle a, such that L ¼ h tan(a). The inclination angle, a, 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 a as large as 20 . The TD inclines towards the h1 100i direction. From this, the average plastic relaxation for a layer is then

III. EXPERIMENTAL PROCEDURE AND ANALYSIS TECHNIQUE
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 cplane (0001). The substrate thicknesses were 300 lm, 630 lm, 420 lm, and 400 lm, respectively.
The radius of curvature was measured via XRD using a Rigaku Smartlab high-resolution x-ray diffractometer. The Rigaku diffractometer is a horizontal h-h 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 29 where x 1 and x 2 are respective positions along the xdirection and x 1 and x 2 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 FIG. 3. Schematic of TD inclination forming at the interface between two layers. A straight dislocation (1) bending at the interface by inclination angle, a, to form an angular dislocation (2).
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 Dx apart on the lab horizontal frame. The angle of the incident beam with the lab horizontal is x, the angle of the incident beam with the crystal vertical is v, the angle between the lab and crystal horizontals is U, and the angle between the crystal vertical and lab horizontal is d.
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, v 1 ¼ v 2 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, b, between the scan positions relative to the center of the R C The assumption is then made that the arc length, c ¼ Rb, 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 Dx, r r , is negligible. Mathematically, this condition would be satisfied if Dx ) r r or Dx ( 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 yaxis to maintain a positive R vector. This shift changes the sign on the Dx 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.

A. XRD curvature measurements
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 6 0.97 m. Substrate B and substrate C had similar radii of curvature, 5.42 6 0.23 m and 6.44 6 0.16 m, respectively. The smallest R C was recorded from substrate D with 2.75 6 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 x for one-dimensional translations where Dx > 0, hence a negative slope in Equation (13) is realized. For simplicity, the R C was represented as an absolute value.

B. Cathodoluminescence measurements
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 5 cm À2 . Substrate B shown in Fig. 7(b) had a TDD of 1.35 Â 10 6 cm À2 , 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 6 cm À2 ; 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 6 cm À2 .

V. THEORETICAL MODEL AND DISCUSSION
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 e Ã 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 e Ã 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. 31 They reported TDD of $5 Â 10 3 cm À2 . 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 10 cm À2 , where the TD was inclined by 10 -25 in the 1 100 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 9 cm À2 . 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 9 cm À2 .
At the TDD observed in this report for HVPE bulk GaN (10 5 -10 6 cm À2 ), it is difficult to measure TDI via crosssectional 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 4 cm À2 -10 6 cm À2 ), 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

VI. CONCLUSIONS
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.