Raman Characterization of the In-Plane Stress Tensor of Gallium Nitride

Experimental characterization of the in-plane stress tensor is a basic requirement for the development of GaN strain engineering. In this work, a theoretical model of stress characterization for GaN using polarized micro-Raman spectroscopy was developed based on elasticity theory and lattice dynamics. Compared with other works, the presented model can give the quantitative relationship between all components of the in-plane stress tensor and the measured Raman shift. The model was verified by a calibration experiment under step-by-step uniaxial compression. By combining the stress characterization model with the expanding cavity model, the in-plane residual stress component field around Berkovich indentation on the (0001) plane GaN was achieved. The experimental results show that the distributions of the stress components, which significantly differed from the distribution of the Raman shift, were closely related to the GaN crystal structure and exhibited a gradient along each crystal direction.


Introduction
Gallium nitride (GaN) is regarded as one of the representatives of the third generation of semiconductor materials due to its fascinating intrinsic properties, such as a wide direct band gap (3.4 eV), strong binding energies and excellent thermal stability and conductivity [1,2], which are widely employed in light-emitting diodes (LEDs), high-electronmobility transistors and high-temperature, high-frequency and high-power semiconductor devices [3][4][5]. Semiconductor devices using GaN materials and related strain engineering require the introduction of state-and size-controllable stress in the GaN layer in the manufacturing process to regulate device functions and to improve the quality and reliability of the devices. Therefore, it is necessary to accurately measure and analyze the complex stress state of GaN materials.
Micro-Raman spectroscopy (MRS) is an effective method to achieve high-resolution [6,7] and nondestructive [8][9][10] detection of stress in materials. Early studies of light-scattering behavior in condensed matter physics found that the Raman spectra of a variety of crystal materials were affected by stress/strain. De Wolf et al. [11,12] gave a simplified quantitative relationship between the Raman shift and stress in single-crystal silicon (Si) and characterized the stress introduced by different processes in a variety of microelectronic structures based on MRS. Lei and Kang et al. [13,14] applied MRS to systematically analyze the characterization model and distribution law of residual stress in porous silicon. Xu et al. [15] investigated the interfacial mechanical parameters of graphene by MRS, and discovered for the first time that there were significant spatial/time scale effects on these interfacial mechanical parameters. Zhao et al. [16,17] carried out experimental studies on the crack extension and fracture toughness of graphene and the interfacial mechanical behavior of bilayer graphene using MRS. In recent years, with the development of GaN semiconductor devices, a series of advances have been made in research on the Raman stress characterization of GaN materials. For example, Briggs et al. [18] gave a relationship between the Raman shift and strain of each mode of GaN and a solution formula for the deformation potential constants of each mode based on linear deformation potential theory. Davydov et al. [19] obtained specific values of deformation potential constants of the partial mode for GaN based on the bulk Grueneisen parameter of the Raman shift of mode during hydrostatic compression under the condition of disregarding the effect of shear strain. Darakchieva et al. [20] explored the relationship between the anisotropic strain of GaN and the mode deformation potential by combining the elliptical polarization method, and calculated the specific value of the deformation potential constants of the mode associated with the tangential strain. Schustek et al. [21] constructed a linear set of equations for the Raman shift-strain relationship for different modes by linearly simplifying the Raman shift-strain equation. Amilusik et al. [22] gave a simplified model of the stress characterization for GaN using the property of the doubly degenerate E 1 and E 2 modes under anisotropic strain.
Based on the above basic research results, Amilusik et al. [22] characterized the residual stress field of ammonothermal grown GaN on substrates in lateral and vertical direction. Hossain et al. [23] analyzed the residual stress in GaN films on Si substrates of different sizes and determined that the stress state of GaN films could be modulated by varying the space size of the Si substrate. Wang et al. [24] investigated the effect of high temperature and mechanical stress on the breakdown voltage of GaN devices and found that both temperature and mechanical stress reduced the performance of the devices by decreasing the breakdown voltage. Park et al. [25] characterized the stress state of GaN films on different substrates, revealing the potential of lithium silicate as an ideal substrate. Chai et al. [26] analyzed the gradient of in-plane residual stress existing in gallium nitride/silicon carbide (GaN/SiC) microstructures of different sizes and found that the residual stress in the GaN film is proportional to the residual stress in the SiC buffer layer and that the stress state in the GaN film can be regulated by modulating the stress state in the SiC layer.
In the published work, the characterization model for Raman stress analysis of GaN could only give the analytical relationship between the sum of the principal stresses and the Raman shift. Moreover, the previous model did not consider the effect of each stress component and could not provide detailed information about the stress state. To address the above problems, this paper developed a model for the stress characterization of GaN based on polarized Raman spectroscopy and established an analytical relationship between the stress component and the measured polarized Raman shift. Our derived model was verified by a calibration experiment to realize the accurate characterization of the complex stress field on the GaN surface.

Materials and Experiments
The sample used in this paper was an undoped N-GaN (negative-gallium nitride) single crystal (Galactic Semiconductor Technology Co., Ltd., Anhui, China). The sample parameters are listed as follows: thickness of 0.374 mm, crystal plane pointing to the (0001) plane, main reference plane of the (1 − 100) plane, dislocation density < 2 × 10 6 cm −2 and etch pit density < 1.5 cm −2 . The 2-inch GaN sample was cut into several rectangular-shaped samples of 10 mm × 2 mm and 10 mm × 10 mm in length and width, respectively, using a scribing machine.
In this paper, uniaxial compression experiments were performed on GaN samples (10 mm × 2 mm × 0.374 mm) using a miniature in situ mechanics tester (CARE IBTC-300S, Care Measurement & control Co., Ltd., Tianjin, China), as shown in Figure 1a. The displacement resolution and load resolution of the mechanics tester were approximately 0.1 µm and 0.1 mN, respectively. The compression experiment was performed with a step loading of 5 N and loaded to 70 N. The Raman scattering signal from the sample surface was acquired three times after each loading using a laser confocal micro-Raman spectrometer (WITEC Alpha 300R, WITec Co., Ulm, Germany), as shown in Figure 1b. The spectral resolution and spectral stability of the Raman spectrometer were approximately 0.8 cm −1 and ±0.02 cm −1 , respectively. The spatial resolution of the displacement stage attached to the Raman spectrometer was approximately 50 nm. When acquiring the Raman scattering signal, the spectrometer parameters were set as follows: 50× objective (numerical aperture, N.A. = 0.55), laser wavelength of 532 nm, laser power of 20 mW, exposure time of 3 s, grating of 1800 L/mm and parallel polarization configuration (polarization directions of the incident laser parallel to the polarization direction of the scattered light). μm and 0.1 mN, respectively. The compression experiment was performed with a step loading of 5 N and loaded to 70 N. The Raman scattering signal from the sample surface was acquired three times after each loading using a laser confocal micro-Raman spectrometer (WITEC Alpha 300R, WITec Co., Ulm, Germany), as shown in Figure 1b. The spectral resolution and spectral stability of the Raman spectrometer were approximately 0.8 cm −1 and ±0.02 cm −1 , respectively. The spatial resolution of the displacement stage attached to the Raman spectrometer was approximately 50 nm. When acquiring the Raman scattering signal, the spectrometer parameters were set as follows: 50× objective (numerical aperture, N.A. = 0.55), laser wavelength of 532 nm, laser power of 20 mW, exposure time of 3 s, grating of 1800 L/mm and parallel polarization configuration (polarization directions of the incident laser parallel to the polarization direction of the scattered light). Nanoindentation experiments were performed on a GaN sample (10 mm × 10 mm × 0.374 mm) using a nanoindentation instrument (Agilent Nano Indenter G200, Agilent Technology Inc., Santa Clara, CA, USA), as shown in Figure 1c. The displacement resolution and load resolution of the nanoindentation instrument were approximately 0.01 nm and 50 nN, respectively. A Berkovich indenter with a radius of curvature of 50 nm was used to squeeze the GaN (0001) surface vertically at room temperature, with one side of the indenter parallel to the (11 00) plane of the sample. The nanoindentation experimental parameters are listed as follows: loading rate of 10 nm/s, allowable drift rate of 0.55 nm/s and maintenance of the maximum load state for 10 s before unloading. The final indentation depth of the nanoindentation experiment was 1300 nm. Then, a Raman mapping experiment was performed around the indentation. The Raman mapping experimental parameters are listed as follows: 100 × objective (N.A. = 0.90), laser wavelength 532 nm, output power of the laser device 20.0 mW (the actual laser power reaching the sample surface in the experiment was about 6.3 mW), exposure time 1 s, grating 1800 L/mm and parallel polarization configuration. The mapping area was a 50 μm × 50 μm square area centered on the indentation, with horizontal and vertical scan steps of 0.5 μm. The schematic is shown in Figure 2. Nanoindentation experiments were performed on a GaN sample (10 mm × 10 mm × 0.374 mm) using a nanoindentation instrument (Agilent Nano Indenter G200, Agilent Technology Inc., Santa Clara, CA, USA), as shown in Figure 1c. The displacement resolution and load resolution of the nanoindentation instrument were approximately 0.01 nm and 50 nN, respectively. A Berkovich indenter with a radius of curvature of 50 nm was used to squeeze the GaN (0001) surface vertically at room temperature, with one side of the indenter parallel to the (1 − 100) plane of the sample. The nanoindentation experimental parameters are listed as follows: loading rate of 10 nm/s, allowable drift rate of 0.55 nm/s and maintenance of the maximum load state for 10 s before unloading. The final indentation depth of the nanoindentation experiment was 1300 nm. Then, a Raman mapping experiment was performed around the indentation. The Raman mapping experimental parameters are listed as follows: 100 × objective (N.A. = 0.90), laser wavelength 532 nm, output power of the laser device 20.0 mW (the actual laser power reaching the sample surface in the experiment was about 6.3 mW), exposure time 1 s, grating 1800 L/mm and parallel polarization configuration. The mapping area was a 50 µm × 50 µm square area centered on the indentation, with horizontal and vertical scan steps of 0.5 µm. The schematic is shown in Figure 2.

Model
GaN crystals have a wurtzite structure (as shown in Figure 3) and belong to the hexagonal system. Its lattice group is C6v 4 , and the basic unit cell is composed of four atoms.

Model
GaN crystals have a wurtzite structure (as shown in Figure 3) and belong to the hexagonal system. Its lattice group is C 6v 4 , and the basic unit cell is composed of four atoms. According to the factorial group analysis, the following relationship exists for the optical phonon modes of wurtzite GaN at the Γ point in the Brillouin zone [27]: where the A 1 , E 1 and E 2 (E 2 low and E 2 high ) modes are Raman active and the B mode is silent. Due to the macroscopic electric field associated with the relative atomic displacement of the longitudinal modes, the A 1 and E 1 modes are split into longitudinal optical (LO) components and transverse optical (TO) components [28].

Model
GaN crystals have a wurtzite structure (as shown in Figure 3) and belong to the hexagonal system. Its lattice group is C6v 4 , and the basic unit cell is composed of four atoms. According to the factorial group analysis, the following relationship exists for the optical phonon modes of wurtzite GaN at the Γ point in the Brillouin zone [27]: where the A1, E1 and E2 (E2 low and E2 high ) modes are Raman active and the B mode is silent. Due to the macroscopic electric field associated with the relative atomic displacement of the longitudinal modes, the A1 and E1 modes are split into longitudinal optical (LO) components and transverse optical (TO) components [28]. Amilusk's study showed that the E1(TO) and E2 high modes are sensitive to biaxial stress in the plane perpendicular to the c-axis, the A1(TO) mode is sensitive to uniaxial stress in the direction parallel to the c-axis and the E2 low mode is rather insensitive to biaxial stress and hydrostatic compression [22]. Figure 4 shows a typical Raman spectrum of the (0001) plane GaN we measured in the vertical backscattering configuration. The characteristic peak corresponding to the A1(TO), A1(LO), E2 low and E2 high modes is observed. The Amilusk's study showed that the E 1 (TO) and E 2 high modes are sensitive to biaxial stress in the plane perpendicular to the c-axis, the A 1 (TO) mode is sensitive to uniaxial stress in the direction parallel to the c-axis and the E 2 low mode is rather insensitive to biaxial stress and hydrostatic compression [22]. Figure 4 shows a typical Raman spectrum of the (0001) plane GaN we measured in the vertical backscattering configuration. The characteristic peak corresponding to the A 1 (TO), A 1 (LO), E 2 low and E 2 high modes is observed. The characteristic peak corresponding to the E 1 mode is not observed, as the E 1 mode is forbidden due to the Raman selection rule in the vertical backscattering configuration of the (0001) plane. Since the signal intensity of the E 2 high mode is higher than that of the E 2 low mode, we chose the E 2 high mode to establish an analytical relationship between the Raman shift and the stress component. characteristic peak corresponding to the E1 mode is not observed, as the E1 mode is forbidden due to the Raman selection rule in the vertical backscattering configuration of the (0001) plane. Since the signal intensity of the E2 high mode is higher than that of the E2 low mode, we chose the E2 high mode to establish an analytical relationship between the Raman shift and the stress component. The general process of establishing the relationship between the Raman shift and stress is shown in Figure 5 [29]. The general process of establishing the relationship between the Raman shift and stress is shown in Figure 5 [29]. The general process of establishing the relationship between the Raman shift and stress is shown in Figure 5 [29]. Expanding the lattice dynamics equation along the Cartesian coordinate system and setting the determinant of its coefficients equal to zero, the characteristic equation of lattice dynamics is obtained: 11 where εuv is the component of the strain tensor ε and Kuvij (i, j = 1, 2, 3) is a component of the mode variable state tensor K, which is used to describe the change in the elastic constant caused by strain. According to the generalized Hooke's law, there is a linear relationship between stress and strain in the case of small deformation: = ε Sσ (3) where σ is the stress tensor and S is the elasticity tensor. Introducing Equation (3) to Equation (2), the eigenvalues λk (k = 1, 2, 3) of the characteristic equation of lattice dynamics can be solved. For the eigenvalues λk, there exists the relationship λ k  ω k 2  ω 0 2 , where ω 0 and ω k are the Raman shifts of the Raman vibration mode in the stress-free state and stressed state, respectively. In the case of small deformation, the Raman shift increment ∆ω k is much smaller than ω 0 , so the following relationship exists: In the actual measurement process, the intensity of each characteristic peak depends on the Raman selection rule, that is, the Raman intensity I of the Raman vibrational mode Expanding the lattice dynamics equation along the Cartesian coordinate system and setting the determinant of its coefficients equal to zero, the characteristic equation of lattice dynamics is obtained: where ε uv is the component of the strain tensor ε and K uvij (i, j = 1, 2, 3) is a component of the mode variable state tensor K, which is used to describe the change in the elastic constant caused by strain. According to the generalized Hooke's law, there is a linear relationship between stress and strain in the case of small deformation: where σ is the stress tensor and S is the elasticity tensor. Introducing Equation (3) to Equation (2), the eigenvalues λ k (k = 1, 2, 3) of the characteristic equation of lattice dynamics can be solved. For the eigenvalues λ k , there exists the relationship λ k = ω 2 k − ω 2 0 , where ω 0 and ω k are the Raman shifts of the Raman vibration mode in the stress-free state and stressed state, respectively. In the case of small deformation, the Raman shift increment ∆ω k is much smaller than ω 0 , so the following relationship exists: In the actual measurement process, the intensity of each characteristic peak depends on the Raman selection rule, that is, the Raman intensity I of the Raman vibrational mode is related to the Raman tensor and the polarization vector of the incident laser and scattered light, as shown in Equation (5): where Q is a constant determined by the scattering law of the sample, the optical properties and the experimental equipment; e I and e S are the polarization vectors of the incident laser and scattered light, respectively; and R k is the Raman tensor of the k mode. The actual measured Raman shift ∆ω obs after the deformation is a linear combination of the mode Raman shift ∆ω k within the same characteristic peak weighted by its respective contributions to the total scattered intensity:  [22], the characteristic equation of lattice dynamics is expressed as: where a E 2 , b E 2 and c E 2 are the phonon deformation potential constants of the E 2 mode [19,20]. Combining Equations (4) and (7), the following relationship is obtained: The Raman tensor corresponding to the E 2 mode is [30]: where d is a constant that depends on the Raman polarization of the sample. The measured GaN surface is in the state of in-plane stress: where σ x and σ y are the normal stress components of the stress tensor σ in the sample coordinate system along the x direction and y direction, respectively, and τ is the shear stress component in the x-y plane. Combining Equations (3) and (10), the stress-strain equation for GaN in the in-plane stress state is obtained: By combining Equations (3), (8), (9) and (11), the relationships between the stress component and the Raman shift of (0001) plane GaN in the polarization configuration is obtained: where HH and HV are the polarization directions of the incident laser, parallel and perpendicular, respectively, to the polarization direction of the scattered light, and α is the angle of polarization of the incident laser.

Calibration Experiment
In this paper, uniaxial compression experiments were conducted using GaN samples to verify the above model. In the stress-free state, the Raman wavenumber of the E 2 high mode we measured is 568.2 cm −1 , which is close to the value measured by others [31]. The Raman shift of the E 2 high mode is the difference between the Raman wavenumber when the E 2 high Materials 2023, 16, 2255 7 of 13 mode is stressed and that when it is stress-free. For the experimental data analysis, we used the three Raman spectra acquired at each load value to calculate the Raman shift separately and averaged them as average Raman shift values. The stress values obtained from the step-loading measurement and the average Raman shift values of the E 2 high mode at this stress value were plotted as the horizontal and vertical coordinates of the experimental data graph, and then the slope (Raman shift-stress coefficient) was obtained by linear fitting to the experimental measurement points, as shown in Figure 6. From the fitted results, the Raman shift-stress coefficient of the (0001) plane GaN in the uniaxial stress state is −2.69 cm −1 /GPa. Based on the relationship between the in-plane stress and the Raman shift of (0001) plane GaN in the HH configuration in Equation (12), the Raman shift-stress coefficient of the E 2 high mode in the same stress state is −2.79 cm −1 /GPa. The comparison shows that the relative error between the experimental results and the theoretical results is less than 3.6%, so the experimental results better verify the correctness of our derived theoretical model. angle of polarization of the incident laser.

Calibration Experiment
In this paper, uniaxial compression experiments were conducted using GaN samples to verify the above model. In the stress-free state, the Raman wavenumber of the E2 high mode we measured is 568.2 cm −1 , which is close to the value measured by others [31]. The Raman shift of the E2 high mode is the difference between the Raman wavenumber when the E2 high mode is stressed and that when it is stress-free. For the experimental data analysis, we used the three Raman spectra acquired at each load value to calculate the Raman shift separately and averaged them as average Raman shift values. The stress values obtained from the step-loading measurement and the average Raman shift values of the E2 high mode at this stress value were plotted as the horizontal and vertical coordinates of the experimental data graph, and then the slope (Raman shift-stress coefficient) was obtained by linear fitting to the experimental measurement points, as shown in Figure 6. From the fitted results, the Raman shift-stress coefficient of the (0001) plane GaN in the uniaxial stress state is −2.69 cm −1 /GPa. Based on the relationship between the in-plane stress and the Raman shift of (0001) plane GaN in the HH configuration in Equation (12), the Raman shift-stress coefficient of the E2 high mode in the same stress state is −2.79 cm −1 /GPa. The comparison shows that the relative error between the experimental results and the theoretical results is less than 3.6%, so the experimental results better verify the correctness of our derived theoretical model. Other scholars had also measured the value of the Raman shift-stress coefficient by experiments. Demangeot et al. [32] gave a value of −2.9 cm −1 /GPa for the Raman shiftstress coefficient, which is closer to our experimental result. This is because the stress state in this work is consistent with the stress state of this paper (uniaxial compressive stress state). The value of the Raman shift-stress coefficient measured by Kisielowski et al. [33]  Other scholars had also measured the value of the Raman shift-stress coefficient by experiments. Demangeot et al. [32] gave a value of −2.9 cm −1 /GPa for the Raman shiftstress coefficient, which is closer to our experimental result. This is because the stress state in this work is consistent with the stress state of this paper (uniaxial compressive stress state). The value of the Raman shift-stress coefficient measured by Kisielowski et al. [33] is −4.2 cm −1 /GPa, which differs from our experimental result. The reason is that the stress state of the measured sample, as well as some material constants used in Raman stress model, are different from those in this paper. Kozawa et al. [34] obtained a value of −6.2 cm −1 /Gpa for the Raman shift-stress coefficient, which is much larger than that obtained in this work. The main reason is that the Olsen and Ettenberg model used in Kozawa's work establishes the relationship between the film stress and the bending deformation based on the numbers of material and structural idealizations. It is not applicable to multilayer problems with unknown or variable elastic parameters or thickness or stress state. The sample in this work had a buffer layer neglected by the authors. Kisielowski et al. [33] showed that the buffer layer affects the residual stress state in the GaN layer.

Raman Shift-Stress Relationship around the Indentation
The study by Yoffe et al. showed that the stress around the indentation is consistent with the expanding cavity model constructed based on the cylindrical coordinate system [35][36][37][38]. Figure 7a gives a microphotograph of a nanoindentation on the (0001) plane GaN sample, and Figure 7b shows the sample coordinate system, crystal coordinate system and cylindrical coordinate system established at the surface center of the nanoin-Materials 2023, 16, 2255 8 of 13 dentation sample. In the cylindrical coordinate system, the z direction coincides with the axial direction of the cylindrical coordinate system in the (0001) plane, r denotes the radial direction and its anticlockwise angle with the x-axis is the tangential direction ϕ.

Raman Shift-Stress Relationship around the Indentation
The study by Yoffe et al. showed that the stress around the indentation is consistent with the expanding cavity model constructed based on the cylindrical coordinate system [35][36][37][38]. Figure 7a gives a microphotograph of a nanoindentation on the (0001) plane GaN sample, and Figure 7b shows the sample coordinate system, crystal coordinate system and cylindrical coordinate system established at the surface center of the nanoindentation sample. In the cylindrical coordinate system, the z direction coincides with the axial direction of the cylindrical coordinate system in the (0001) plane, r denotes the radial direction and its anticlockwise angle with the x-axis is the tangential direction φ.
where r, φ and z denote the radial direction, tangential direction and axial direction, respectively, under the cylindrical coordinate system; σrr is the radial stress component, σφφ is the tangential stress component and σrφ is the shear stress component; νrφ and νφr are Poisson's ratios and Grz and Gzφ are the shear moduli, which can be calculated using the elasticity tensor S ' in the cylindrical coordinate system; Brr and Bφφ are functions characterizing the degree of material bulging caused by plastic deformation near the indentation; fφφ is a function characterizing the nonconservative change in volume after plastic deformation of the material near the indentation; P is the peak load applied to the indenter; and Based on the expanding cavity model, each component of the in-plane stress state at any point around the indentation is expressed as follows [37]: where r, ϕ and z denote the radial direction, tangential direction and axial direction, respectively, under the cylindrical coordinate system; σ rr is the radial stress component, σ ϕϕ is the tangential stress component and σ rϕ is the shear stress component; ν rϕ and ν ϕr are Poisson's ratios and G rz and G zϕ are the shear moduli, which can be calculated using the elasticity tensor S in the cylindrical coordinate system; B rr and B ϕϕ are functions characterizing the degree of material bulging caused by plastic deformation near the indentation; f ϕϕ is a function characterizing the nonconservative change in volume after plastic deformation of the material near the indentation; P is the peak load applied to the indenter; and H is the Vickers hardness of the crystal. Since the radial shear stress component is 0, the stress at any point around the indentation given by the expanding cavity model is actually in a state of non-equal biaxial stress, and the directions of each principal stress are known. The ratio of the radial stress component to the tangential stress component is: From the definition of Poisson's ratio and the shear modulus, it follows that: (15) where S 11 , S 12 , S 21 , S 22 , S 55 and S 66 are the components of the elasticity tensor S in the cylindrical coordinate system. Combining Equations (14) and (15), and substituting the material parameters of the crystal plane, the ratio of the radial to the tangential stress component is obtained as follows: Combined with the rotation axis equation of the stress component and Equation (12) in this paper, the relationship between the measured Raman shift of the E 2 mode of GaN in the HH configuration and σ rr and σ ϕϕ in the cylindrical coordinate system is obtained: Combining Equation (16) with Equation (17), the relationship between each stress component and the Raman shift at any position in the region around the GaN indentation is obtained as follows:

Distribution of Residual Stress Component around the Indentation
In this paper, a Raman mapping experiment was performed on the region around the indentation of the GaN sample. The Raman peak of the E 2 high mode was fitted using the Lorenz function, and the experimental data were analyzed by calculating the Raman shift of the E 2 high mode. The distribution of the Raman shift of the E 2 high mode around the (0001) plane indentation is shown in Figure 8. The distribution diagram of the Raman shift shows that there is a region of increasing shift near the indentation with a certain directionality. Specifically, the Raman shift significantly increases at the center of the indentation along the [2   By substituting the measured Raman shift into Equation (18), the distribution of the residual stress components in the region around the indentation of the (0001) plane GaN are obtained, as shown in Figure 9. The residual stress components have obvious distribution characteristics along the [21 1 0], [1 21 0], [11 00] and [112 0] directions, indicating that the distribution trend is related to the lattice structure of the crystalline material. Relatively, the gradient of stress variation along the [11 00] and [112 0] directions is more obvious. The distribution of stress values indicates that σrr is negative (indicating compressive stresses) and that σφφ is positive (indicating tensile stresses). The extreme values of σrr and σφφ appear on the right side of the vertical side of the indentation near the center, which are approximately −4.20 GPa and 0.90 GPa, respectively. At the same point of the sample, the value of σrr is several times larger than that of σφφ, further confirming the dominance of residual compressive stresses on the sample surface. Figure 9. Distributions of (a) radial stress σrr and (b) tangential stress σφφ. Figure 9. Distributions of (a) radial stress σ rr and (b) tangential stress σ ϕϕ .
According to the distributions of σ rr and σ ϕϕ in Figure 9, the distributions of normal stresses σ x and σ y and shear stress τ on the (0001) plane GaN in the rectangular coordinate system can be obtained by using the rotation axis equation, as shown in Figure 10  As we previously mentioned, most scholars have applied the relationship between the Raman shift and the sum of the principal stresses for stress characterization of GaN. A comparison of Figure 8 with Figure 9 shows that there is a large difference between the distribution of the Raman shift and the distribution of the stress components, which illustrates the limitations of the previous analytical ideas of stress characterization. If we also use the previous method, then we will obtain a stress distribution consistent with the distribution of the Raman shift in Figure 8. This approach will not only give the incorrect trend of stress distribution but also increase the numerical level of the overall stress value with the calculation error up to 174%. In addition, that method cannot give the directions of true stress at any point of the sample. Using our derived model for stress characterization in combination with the expanding cavity model, it is possible to give the stress values and stress directions at any point in the region around the indentation on the (0001) plane GaN. Therefore, the advantage of the stress characterization model around the indentation obtained by considering the stress state on the sample surface is that it can give the true stress information and then realize the fine characterization of the stress component field.
The above discussion proved that it should be meaningful to consider the influence of stress state when establishing the Raman stress model of GaN. If the stress state of the sample was not considered, that is, using the relationship between the Raman shift and the sum of the principal stresses to analyze the stress distribution trend, it would not only Figure 10. Distributions of (a) normal stress σ x , (b) normal stress σ y and (c) shear stress τ.
As we previously mentioned, most scholars have applied the relationship between the Raman shift and the sum of the principal stresses for stress characterization of GaN. A comparison of Figure 8 with Figure 9 shows that there is a large difference between the distribution of the Raman shift and the distribution of the stress components, which illustrates the limitations of the previous analytical ideas of stress characterization. If we also use the previous method, then we will obtain a stress distribution consistent with the distribution of the Raman shift in Figure 8. This approach will not only give the incorrect trend of stress distribution but also increase the numerical level of the overall stress value with the calculation error up to 174%. In addition, that method cannot give the directions of true stress at any point of the sample. Using our derived model for stress characterization in combination with the expanding cavity model, it is possible to give the stress values and stress directions at any point in the region around the indentation on the (0001) plane GaN. Therefore, the advantage of the stress characterization model around the indentation obtained by considering the stress state on the sample surface is that it can give the true stress information and then realize the fine characterization of the stress component field.
The above discussion proved that it should be meaningful to consider the influence of stress state when establishing the Raman stress model of GaN. If the stress state of the sample was not considered, that is, using the relationship between the Raman shift and the sum of the principal stresses to analyze the stress distribution trend, it would not only give an incorrect conclusion, but also achieve an error of stress value up to 174%. Meanwhile, based on the model presented in this work, the decoupling of the stress components of the (0001) plane GaN can be realized under the condition that the principal stress direction is known, which will be beneficial for further studies to realize the complete decoupling of the stress components of GaN in any arbitrary stress state. At the same time, the GaN stress characterization model proposed in this paper can solve most of the engineering problems. This is because in practical engineering problems, there usually exist positions with known stress states (such as free boundaries) which can be used as boundary conditions. Then the continuity relationship given by elastic mechanics should be used to derive the principal stress value and its direction at the analyzed location, thus realizing stress decoupling. Of course, such an analysis is required to combine Raman measurement with the mathematical or simulation modeling of elastic mechanics, hence it is not a direct measurement. Moreover, this paper can provide an analytical idea for the establishment of a Raman stress model of other material. In addition, the GaN samples we used in our experiments were not intentionally doped, so we did not consider the effect of doping in the analysis of stress, nor in the theoretical model.

Conclusions
This work, by establishing the GaN stress characterization model under polarized Raman, presented the relationship between the Raman shift and each stress component. Compared with the traditional method using the relationship between the Raman shift and the sum of the principal stresses, the proposed model considered the influence of stress state. The Raman shift-stress coefficient was achieved by the calibration experiment. By combining the stress characterization model with the expanding cavity model, the residual stress component field around the indentation was achieved based on Raman experiments in the HH configuration. The distributions of the residual stress components showed a close correlation with the crystal structure of the sample material, furtherly verified the importance of considering the stress state. This work would be beneficial to the strain engineering and stress characterization of GaN materials and their related devices.

Data Availability Statement:
The data that support the findings of this study are available from the corresponding author upon reasonable request.

Conflicts of Interest:
The authors declare no conflict of interest.