Identification of the Kirkendall effect as a mechanism responsible for thermal decomposition of the InGaN/GaN MQWs system

Our calculations were performed in the framework of Kohn–Sham DFT by means of SIESTA program [22, 23]. The exchange-correlation interaction was described by the general gradient approximation, and we used a version of the Perdew–Burke–Ernzerholf potential with parameters β, µ and κ fixed by the jellium surface (Js), jellium response (Jr), and Lieb–Oxford bound criteria, respectively [24, 25]. The electron–ion core interactions were represented by Troullier–Martins type pseudopotentials and the electron wavefunctions were expanded into the atomic orbital basis set using either single-ζ, triple-ζ, or triple-ζ polarized (N: $2s^\mathrm{TZ}$, $2p^\mathrm{TZP}$; Ga: $4s^\mathrm{TZ}$, $4p^\mathrm{TZ}$, $3d^\mathrm{SZ}$; In: $5s^\mathrm{TZ}$, $5p^\mathrm{TZ}$, $4d^\mathrm{SZ}$, $5d^\mathrm{SZ}$). The cutoff energy of 500 Ry was used to produce a real-space mesh, and Brillouin zone sampling was carried out using a $15\times15\times9$ Monkhorst–Pack k-point grid for a primitive unit cell of bulk GaN. The positions of all atoms in the GaN unit cell were relaxed with the convergence criteria imposed on the atomic forces being less than 0.001 eV Å⁻¹. This geometry optimization led to the following lattice constants for bulk GaN: $a = b = 3.215$ Å and c = 5.220 Å. They can be compared with the experimental data obtained by x-ray diffraction $a = b = 3.189$ Å and c = 5.185 Å [26].

An alloy $\mathrm{In_\textit xGa_{1-\textit{x}}N}$ (x = 0.11 and x = 0.22) was represented by means of a supercell model of size $3\times3\times3$. Inside each GaN metal layer, a certain number of In dopant atoms were randomly placed at the lattice sites to reflect the compositions most similar to the cases of $\mathrm{In_{0.125}Ga_{0.875}N}$ and $\mathrm{In_{0.25}Ga_{0.75}N}$. The geometric model of the supercell with 108 atoms is shown in figure 1. All atomic positions of this alloy were optimized again with the force restriction imposed to be less than 0.002 eV Å⁻¹. The prepared supercells were used to calculate the heights of the migration energy barriers of Ga and In vacancies along the hexagonal c axis by means of the nudged elastic band method [27–29].

Zoom In Zoom Out Reset image size

To investigate the thermal properties of the discussed system, such as the temperature dependence of the diffusion coefficients of In and Ga atoms, we performed phononic calculations using a direct method in harmonic approximation. This means that phonon frequencies are considered volume independent, and the thermal expansion coefficient is equal to zero at all temperatures [30]. For phononic calculations, a larger supercell was built by multiplying the structure presented in figure 1 by a factor $3\times3$, which led to a system $\mathrm{In_\textit xGa_{1-\textit{x}}N}$ of 963 atoms. This larger system was used to calculate force constant matrices by means of the finite-displacement method. Only the Γ point sampling was used at this stage of the calculations. A displacement of 0.03 Bohr (0.016 Å) was applied in each direction $(-x, +x, -y, +y, -z, +z)$ for every atom in a centrally located cell. In the harmonic approximation, the obtained force constant matrix is transformed into the dynamical matrix D, and the vibration frequencies of the atoms are determined by solving an eigenvalue problem of the matrix D [30–32].

Based on the vibrational spectra, several properties, such as the zero-point vibrational energy $E^\mathrm{ZP}$, as well as the vibrational energy $E^\mathrm{vib}$, the vibrational entropy $S^\mathrm{vib}$ and the free energy $F^\mathrm{vib}$ can be determined as functions of temperature T [33]:

$$\begin{equation} E^\mathrm{ZP} = \sum_j \frac{\hbar\omega_j}{2}, \end{equation} \tag{ 1 }$$

$$\begin{equation*} E^\mathrm{vib} (T) = E^\mathrm{ZP} + k_\mathrm B T \sum_j \frac{x_j}{\exp(x_j)-1}, \end{equation*}$$

$$\begin{equation} \textrm{where} \hspace{0.6cm} x_j = \frac{\hbar \omega_j}{k_\mathrm B T} , \end{equation} \tag{ 2 }$$

$$\begin{equation} S^\mathrm{vib}(T)\! = \! k_\mathrm{B}\! \sum_{j} \!\left\lbrace \frac{x_j}{\exp(x_j)\! -\! 1}\! -\! \ln\left[1\! -\! \exp(-x_j)\right] \right\rbrace, \end{equation} \tag{ 3 }$$

$$\begin{equation} F^\mathrm{vib} (T) = E^\mathrm{ZP} + k_\mathrm BT \sum_j \ln[1-\exp(-x_j)], \end{equation} \tag{ 4 }$$

where $\hbar$, k_B , and ω_j are the reduced Planck constant, the Boltzmann constant, and the angular frequency of each jth mode of vibration, respectively. In order to apply equations (1)–(4) for crystals, it is necessary to include a frequency dependence on the wave vector, that is, $\omega_j (\vec{k})$. Thus, integration over the first Brillouin zone or summation over special k-points must be performed [34, 35]. In this study, the sampling of the k-space was tested and finally carried out on a grid of $4\times4\times4$, which is equivalent to 64 points.

Usually, the temperature dependence of the diffusion coefficient D is consistent with the Arrhenius formula [36]:

$$\begin{equation} D = D^{0} \exp \left(-\frac{\Delta H}{k_\mathrm B T}\right) , \end{equation} \tag{ 5 }$$

where $\Delta H$ is the activation enthalpy of diffusion and D⁰ denotes the pre-exponential factor also called the frequency factor. This is an empirical relationship according to which the diffusion coefficient is described at the macroscopic level by fitting to the experimental results.

At the atomic level description, the diffusion coefficient D is usually written in the following form:

$$\begin{equation} D = f\, l^2~\Gamma. \end{equation} \tag{ 6 }$$

This means that D is proportional to the jump rate constant Γ and the distance l that the atom moves to a neighboring site, multiplied by any correlation factor f that includes the number of possible jump sites in a given lattice type and diffusion dimension. For a hexagonal closed packed hcp lattice, the factor f equals 0.781 [37]. In the case of the vacancy mediated mechanism, the probability that a neighboring location could accept a particle should also be taken into account, so the vacancy concentration c_V should be included as the scaling factor of the coefficient $D^{*} = c_\mathrm V D$. Vacancy concentration is defined as the ratio between empty lattice sites and all lattice sites in a system. It is calculated according to the expression:

$$\begin{equation} c_\mathrm V = \exp \left(-\frac{\Delta G_f}{k_\mathrm B T}\right), \end{equation} \tag{ 7 }$$

where $\Delta G_f$ is the Gibbs energy required to create a vacancy. In general, the formation energy depends on the Fermi level position and the charge state of the defect. Because our goal is to compare the mobility of the gallium and indium atoms in the same material at any vacancy concentration, we will focus our efforts on determining the part related to the barriers to atomic hopping. Therefore, in the following considerations, we shall not discuss the vacancy formation energy and vacancy concentrations.

Individual jumps of atoms are promoted by thermal activation and can be described by the Eyring equation derived from the transition state theory [38]:

$$\begin{equation} \Gamma = \nu^{0} \exp\left({-\frac{\Delta G}{k_\mathrm B T}}\right) = \nu^{0} \exp\left({\frac{\Delta S}{k_\mathrm B}}\right) \exp\left({-\frac{\Delta H}{k_\mathrm B T}}\right), \end{equation} \tag{ 8 }$$

where the prefactor $\nu^{0} = \frac{k_\mathrm B T}{h}$ is a universal frequency and $\Delta G$ is the difference in the Gibbs free energy between the barrier saddle point and the equilibrium position. According to thermodynamic rules, it can be separated into enthalpy and entropy terms: $\Delta G = \Delta H - T\Delta S$. Equation (8) has been shown to be consistent with the simplified Vineyard formula [39, 40]:

$$\begin{equation} \Gamma = \frac{\prod \limits_{i}^{N} \nu_i}{\prod \limits_{j}^{N-1} \nu_j} \exp\left({-\frac{\Delta H}{k_\mathrm B T}}\right), \end{equation} \tag{ 9 }$$

where ν_i and ν_j are the vibration frequencies in the initial and transition states and N is the total number of frequencies. By comparing expressions (6) and (8) with (5), the prefactor D⁰ can be written as the part depending on the change in entropy during the migration of the atom in the system. The entropy of diffusion corresponds to the change in lattice vibrations associated with the displacement of the atom from its equilibrium to the saddle point configuration.

Using DFT calculations, it is possible to determine the changes in the total energy of the system ($\Delta E^\mathrm{DFT}$) between the transition state and the equilibrium point, that is, the energy barrier ($E_\mathrm b$) for atom migration. Additionally, phonon calculations make it possible to determine the frequencies of vibrations of atoms and, on their basis, to determine the vibrational free energy ($F^\mathrm{vib}$) as a function of temperature. To simplify the calculations, pressure-induced effects in the crystal are neglected, and only a constant-volume system is considered here. Thus, an approximation of the Gibbs energy changes by the Helmholtz free energy can be used:

$$\begin{equation} \Delta G(p,T) \approx \Delta F(V,T) = \Delta E^\mathrm{DFT} + \Delta F^\mathrm{vib}(T) \equiv E_\mathrm b(T). \end{equation} \tag{ 10 }$$

Thus, the jump rate equation based on the data from the DFT calculations will be as follows:

$$\begin{equation} \Gamma = \nu^{0} \exp\left({-\frac{E_\mathrm{b}(T)}{k_\mathrm B T}}\right) = \nu^{0} \exp\left({-\frac{\Delta E^\mathrm{DFT}}{k_\mathrm B T}}\right) \exp\left({-\frac{\Delta F^\mathrm{vib}(T)}{k_\mathrm B T}}\right). \end{equation} \tag{ 11 }$$

Dependence (11) is consistent with the Eyring and Vineyard equations.

In order to calculate the diffusion coefficients of Ga and In atoms in $\mathrm{In_\textit xGa_{1-\textit{x}}N}$ alloys, it is necessary to find the heights of the migration energy barriers of the corresponding point defects (see equations (6) and (11)). We assumed the diffusion of Ga and In atoms took place in the c direction of the hexagonal supercell, via the vacancy mediated mechanism, see figure 1.

The first assumption is based on the fact that the experimental papers discussed in the Introduction [2–5, 11–15] concern single- or multi-QW InGaN/GaN systems, grown just along the hexagonal c direction, and that the thermal decomposition of the InGaN/GaN QWs system occurred also in this direction. The thermal decomposition of the InGaN/GaN MQW system observed by us [13, 14] starts from the first QW, and later, with the increase in temperature and time, the next QW in the c direction is decomposed, up to the last one (see figure 12 [13]). Therefore, the direction c is crucial from the experimental point of view. The second assumption concerns the choice of the diffusion mechanism. It is known that the vacancy mediated diffusion mechanism is based on a usual jump between an empty lattice site and the site occupied by an atom, and that this mechanism is of particular importance in the situation when atoms migrate within the same sublattice, e.g. metal atoms having a similar ionic radius (In or Ga) in the metal sublattice. This mechanism is very important [41–43] in nitrides. It also follows from the experiment that in the case of thermal decomposition of the $\mathrm{In_\textit xGa_{1-\textit{x}}N/GaN}$ QWs, there is a high concentration of vacancies in the structure [2–5, 11–15]. Moreover, as was mentioned in the Introduction, the Kirkendall effect, which we propose as the mechanism responsible for the thermal decomposition of the $\mathrm{In_\textit{x}Ga_{1-\textit{x}}N/GaN}$ QWs, is based just on this mechanism. Therefore, we chose the vacancy mediated diffusion mechanism as the main one to be explored herein.

The calculated heights (by means of the nudged elastic band method) of the migration energy barriers $E_\mathrm b$ of vacancies $V_\mathrm{Ga}$ and $V_\mathrm{In}$ in $\mathrm{In_\textit xGa_{1-\textit{x}}N}$ alloys with various compositions ($x = 0, 0.11, 0.22$) are presented in table 1.

Table 1. Heights of the migration energy barriers of vacancies $V_\mathrm{Ga}$ and $V_\mathrm{In}$ in alloys $\mathrm{In_\textit xGa_{1-\textit{x}}N}$.

% of In in	Energy barrier $E_\mathrm b$, eV
	$V_\mathrm{Ga}$		$V_\mathrm{In}$
$\mathrm{In_\textit xGa_{1-\textit{x}}N}$	This paper	Others	This paper	Others
0	2.81	2.75, 2.8 [42, 44]	2.10	2.10 [44]
0.11	2.78		2.01
0.22	2.69		1.37

As can be seen in table 1, the height of the energy barrier of $V_\mathrm{Ga}$ in bulk GaN is in good agreement with the calculations presented in other papers [42, 44]. With the increase in the concentration of In atoms in alloys of $\mathrm{In_\textit xGa_{1-\textit{x}}N}$ up to 11% and 22%, the heights of the migration energy barriers of $V_\mathrm{Ga}$ decrease slightly. Figure 2 presents a comparison of the energy barriers profiles of the Ga and In atoms migrating in alloys $\mathrm{In_\textit xGa_{1-\textit{x}}N}$. When the concentration of In atoms is equal to 11%, the energy barriers related both to the migration of Ga and In atoms are lowered insignificantly. However, in the case of x = 0.22 an essential decrease in the height of the energy barrier of the migrating In atom is observed (see figure 2(b)). This indicates a possible greater diffusion of In atoms compared to Ga atoms in alloys $\mathrm{In_\textit xGa_{1-\textit{x}}N}$. It follows from the shape of the energy barriers presented in figures 2(a) and (b) that the local environment of the migrating Ga and In atoms has no significant influence on the initial and final state of the diffusion process, even at higher dopant concentrations. It can be proven by the fact that the position of the base and the bottom of barrier on the energy scale corresponds to the value of 0 eV, in most analyzed cases. It means that when the migrating Ga or In atom starts the diffusion from the lattice site (lattice node) corresponding to the initial state of the diffusion, it finishes it also in the lattice site (lattice node), corresponding to the final state of the diffusion. Only for the Ga atom migrating in the $\mathrm{In_{0.22}Ga_{0.78}N}$ alloy we can observe a small influence of the local environment on the final state of the diffusion of this atom. This migrating Ga atom after overcoming the energy barrier of 2.69 eV in the saddle point (transition state) returns not to the lattice site located to the lattice node, but in its neighborhood.

Zoom In Zoom Out Reset image size

The decrease in the energy barrier that is observed with the increase in indium concentration in the $\mathrm{In_\textit xGa_{1-\textit{x}}N}$ systems can be easily explained by the geometry relaxation of the system. As is known, the lattice constant of GaN $a \mathrm{[GaN]} = 3.190$ Å which is smaller than that of a[InN] = 3.545 Å. The $\mathrm{In_\textit xGa_{1-\textit{x}}N}$ structure is therefore strained. A manifestation of internal strain is the value of pressure inside the relaxed system. Our numerical calculations show that in the case of the relaxed $\mathrm{In_{0.11}Ga_{0.89}N}$ system, the internal pressure is 76.3 kbar, while in the case of $\mathrm{In_{0.22}Ga_{0.78}N}$ this pressure increases to the value of 163.9 kbar. The introduction of the single-point defect $V_\mathrm{Ga}$ or $V_\mathrm{In}$ into the $\mathrm{In_{0.11}Ga_{0.89}N}$ system and its relaxation leads to a decrease in pressure by about 17 kbar, while the introduction of an analogous single point defect into the $\mathrm{In_{0.22}Ga_{0.78}N}$ system and its relaxation results in a double decrease in pressure, i.e. by about 36 kbar. This situation occurs when defect $V_\mathrm{Ga}$ or $V_\mathrm{In}$ is introduced into the lattice node (initial or final states). When an In atom migrates in a defected $\mathrm{In_{0.11}Ga_{0.89}N}$ system and reaches the saddle point (transition state), then the decrease in pressure of the relaxed system reaches 21 kbar, but when an In atom migrates in the analogous $\mathrm{In_{0.22}Ga_{0.78}N}$ system and reaches the transition state, then the decrease in pressure is down to 44 kbar. This is, in our opinion, the main physical reason for the decrease in the migration energy barrier observed with the increase in the indium concentration in the $\mathrm{In_\textit xGa_{1-\textit{x}}N}$ systems.

Additional, essential information about the possibility of diffusion of gallium and indium atoms can be provided by the temperature dependence of $E_\mathrm b(T)$, which will be analyzed in the next section.

Next, we computed the phonon spectra of $\mathrm{In_\textit xGa_{1-\textit{x}}N}$ systems in the presence of migrating Ga or In atoms in the crystal lattice between the initial, transition (saddle point) and final sites, via the vacancy mediated mechanism. To verify the phonon calculations of the $\mathrm{In_\textit xGa_{1-\textit{x}}N}$ systems, we calculated the phonon spectrum of a bulk GaN crystal which is presented in figure 3. As can be seen, there is good agreement between the data obtained and other theoretical and experimental results [33, 45].

Zoom In Zoom Out Reset image size

Point defects such as vacancies or substitutional atoms present in the crystal lattice make a shift in the frequency of vibrational modes of the phonon spectrum in comparison to the spectrum of an ideal bulk crystal. This is related to a local change in the forces that act between the atoms. Figure 4 presents the phonon spectra of the system $\mathrm{In_{0.22}Ga_{0.78}N}$, where a migration of Ga and In atoms occurred between the initial and final lattice sites, throughout the saddle point (transition state). As can be found from figures 4(b) and (d) a single negative mode appears in the vibrational spectrum of the system only when the migrating atom is in the saddle point (transition state). The spectra presented in figure 4 corroborate both the appropriate convergence parameters and the supercell size chosen for the calculation [40, 46]. The frequencies of these negative modes of the systems studied were further eliminated in calculations of the jump rates Γ of the Ga and In atoms, according to equation (11) [40].

Zoom In Zoom Out Reset image size

Next, we tested various densities of k-point grids necessary to calculate the vibrational free energy and diffusion coefficient, using the standard procedure, presented in [40]. The tests were carried out on grids of points $1\times1\times1$, $2\times2\times2$, $4\times4\times4$, $6\times6\times6$ and $8\times8\times8$. Figure 5 presents a temperature dependence of the vibrational free energy difference $\Delta F^\mathrm{vib}$ of the Ga atom that migrates between the transition and initial states in GaN's supercell, taking into account the zero-point vibrational energy difference $\Delta E^\mathrm{ZP}$ (see equation (4)). As can be traced in figure 5, the choice of k point density does not influence only the zero-point vibrational energy difference; however, it is essential at higher temperatures [40], starting from 300 K. Note also that the differences in the values of vibrational free energy $\Delta F^\mathrm{vib}$ calculated for the $4\times4\times4$, $6\times6\times6$, and $8\times8\times8 \;k$-point grids almost coincide for temperatures higher than 300 K. The optimal density of the k-point grid to calculate the vibrational free energy of the $\mathrm{In_\textit xGa_{1-\textit{x}}N}$ system is therefore, in our opinion, $4\times4\times4$, which corresponds to 64 k points.

Zoom In Zoom Out Reset image size

As mentioned in section 3.1, essential information on the possibility of diffusion of Ga and In atoms in $\mathrm{In_\textit xGa_{1-\textit{x}}N}$ alloys can be provided by the shape of the temperature dependence of the migration energy barrier $E_\mathrm b(T)$ (see equation (10)). For the selected grid of $4\times4\times4 \;k$ points, we calculated $E_\mathrm b(T)$ of Ga and In atoms in $\mathrm{In_\textit xGa_{1-\textit{x}}N}$ alloys that are presented in figure 6. As can be seen, the behavior of both barriers is similar to that observed in T = 0 K, that is, $E_\mathrm b(\mathrm{In}) \lt E_\mathrm b(\mathrm{Ga})$, and their heights are essentially decreased with increasing concentration of In atoms (compare figure 2). Moreover, each barrier decreases with temperature; however, the functions $E_\mathrm b(T)$ exhibit different slopes for particular In atom concentrations in $\mathrm{In_\textit xGa_{1-\textit{x}}N}$. For the temperature T = 1200 K corresponding to the decomposition of QWs in alloys $\mathrm{In_\textit xGa_{1-\textit{x}}N}$ [13, 14], the $E_\mathrm b$ of In is equal to 0.87 eV in $\mathrm{In_{0.22}Ga_{0.78}N}$. As follows from table 1, the initial value of this barrier was 1.37 eV in T = 0 K. In the case of the Ga atom that migrates in this alloy, $E_\mathrm b(1200\, \mathrm{K}) = 2.26$ eV, which means that the initial energy barrier decreased only by 0.43 eV. Therefore, we can state that the increase in the temperature in the system affects the diffusion of In atoms in the alloys $\mathrm{In_\textit xGa_{1-\textit{x}}N}$.

Zoom In Zoom Out Reset image size

As mentioned in the Introduction, the phenomenon of thermal decomposition of $\mathrm{In_\textit xGa_{1-\textit{x}}N}$ QWs has been intensively investigated over the past years [2–5, 11–15]. The last experimental results of our group show that the diffusion of point defects is a dominant factor that contributes to this process [13, 14, 47].

In order to fully explain the mechanism responsible for the decomposition process, it is necessary to take into account both the heights of the migration energy barriers of particular point defects and the diffusion coefficients of atoms migrating in the crystal lattices of InGaN and GaN. Recent theoretical calculation results show that the heights of the migration energy barriers of In and Ga atoms in bulk GaN are almost twice smaller than those of N atoms; $E_\mathrm b(\mathrm{N}) \gt E_\mathrm b (\mathrm{Ga}) \gt E_\mathrm b(\mathrm{In})$, see table 1 and [42, 44]. Furthermore, Harafuji [48] et al using molecular dynamics simulations showed that the self-diffusion coefficient of Ga atoms in bulk GaN was $3\times10^{-9}$ and $2\times10^{-7}\mathrm{cm^2s^{-1}}$ at 1500 K and 2000 K, respectively, while the self-diffusion coefficient of N atoms was less than $1\times10^{-9}\mathrm{cm^2s^{-1}}$ at 2000 K. This means that a nitrogen vacancy is less mobile and plays the role of a base for Ga movement through the generation and disappearance of Ga–N–associated vacancies [48]. The last statement was also confirmed by Ambacher et al [49]. They measured the self-diffusion coefficient of nitrogen versus temperature in the temperature range 770 ${}^{\circ}$C–970 ${}^{\circ}$C in a hexagonal GaN crystal, based on the results of concentration depth profiles determined from SIMS signals. Their results show that the nitrogen self-diffusion coefficient at 800 ^∘C and 900 ${}^{\circ}$C was $9.3\times10^{-17}\mathrm{cm^2s^{-1}}$ and $7.3\times10^{-15}\mathrm{cm^2s^{-1}}$ [49]. The high stability of the nitrogen sublattice of GaN was also confirmed by Sadovyi et al by determining the low values of the oxygen diffusion coefficients $10^{-12}\mathrm{cm^2s^{-1}}$ and $10^{-13}\mathrm{cm^2s^{-1}}$ at 2775 K and 3375 K, respectively [50].

Based on equations (6) and (11) we calculated the diffusion coefficients of the Ga and In atoms, $D_\mathrm{Ga}$ and $D_\mathrm{In}$, in alloys $\mathrm{In_\textit xGa_{1-\textit{x}}N}$ with various compositions. The temperature dependencies of $D_\mathrm{Ga}$ and $D_\mathrm{In}$ are presented in figure 7. As can be observed, an In atom is more mobile in the $\mathrm{In_{0.22}Ga_{0.78}N}$ alloy than in bulk GaN, while the Ga atom is slightly more mobile in bulk GaN than in alloys $\mathrm{In_\textit xGa_{1-\textit{x}}N}$ ($x = 0.11, 0.22$). Moreover, indium atoms diffuse faster than gallium atoms both in $\mathrm{In_\textit xGa_{1-\textit{x}}N}$ alloys and in GaN. It follows from our calculations that the diffusion coefficient $D_\mathrm{In}$ of indium atoms in $\mathrm{In_{0.22}Ga_{0.78}N}$ is by four orders higher than that of Ga atoms in GaN. For example, at 1200 K, $D_\mathrm{In}(\mathrm{In_{0.22}Ga_{0.78}N}) = 9.29\times10^{-7}\mathrm{cm^2s^{-1}}$, $D_\mathrm{In}(\mathrm{GaN}) = 4.85\times10^{-10}\mathrm{cm^2s^{-1}}$, $D_\mathrm{Ga}(\mathrm{In_{0.22}Ga_{0.78}N}) = 7.76\times10^{-12}\mathrm{cm^2s^{-1}}$, and $D_\mathrm{Ga}(\mathrm{GaN})=$ $3.18\times10^{-11}\mathrm{cm^2s^{-1}}$. This means that at the annealing temperature $940\,^{\circ}$C of $\mathrm{In_{0.18}Ga_{0.82}N/GaN}$ MQWs system [13] approximated value of the indium atom diffusion coefficient, $D_\mathrm{In}$, is of the order of ${\sim}1.5\times10^{-8}\mathrm{cm^2s^{-1}}$, and that of the Ga atom, $D_\mathrm{Ga}$, in the GaN barrier area is $1.4\times10^{-11}\mathrm{cm^2s^{-1}}$. The results obtained above agree well with the data reported by Harafuji et al [48] for GaN.

Zoom In Zoom Out Reset image size

Based on the results presented in section 3.3, the corresponding diffusion fluxes of metal atoms between the GaN barrier area and the InGaN QW can be presented schematically in figure 8(upper panel). During the annealing of two stacked materials, GaN and $\mathrm{In_\textit xGa_{1-\textit{x}}N}$, at a temperature high enough to thermally activate the diffusion of the atoms, the migration of the atoms can occur at the interface of these materials, where the atoms will diffuse from $\mathrm{In_\textit xGa_{1-\textit{x}}N}$ to GaN and vice versa. Since the In atom diffuses into InGaN much faster than the Ga atom in GaN, the flux density of the In atoms that migrate from $\mathrm{In_\textit xGa_{1-\textit{x}}N}$ to GaN ($J_\mathrm{In}$) will be much higher than the flux density of the Ga atoms that diffuse in the opposite direction ($J_\mathrm{Ga}$). As a result of such an annealing process, a layer of $\mathrm{GaN/In_\textit xGa_{1-\textit{x}}N}$ alloy is formed with a higher In content, compared to the initial In concentration in the $\mathrm{In_\textit xGa_{1-\textit{x}}N}$ QW, located between the two sides of the phase boundary. The final thickness of this $\mathrm{GaN/In_\textit xGa_{1-\textit{x}}N}$ layer depends on both the temperature and the annealing time. It is important that as a consequence of the unbalanced diffusion rates of atoms between the $\mathrm{In_\textit xGa_{1-\textit{x}}N}$ and GaN materials, vacancies will be introduced into the interfacial area within $\mathrm{In_\textit xGa_{1-\textit{x}}N}$, whose flux density in figure 8 is denoted $J_\mathrm V$. Therefore, the alloy region with a high In content will be more extended within the GaN barrier region, and vacancies will be introduced in the interface region in the $\mathrm{In_\textit xGa_{1-\textit{x}}N}$ QW.

Zoom In Zoom Out Reset image size

Over time, the coalescence of excess vacancies will lead to the formation of small voids spread throughout the entire $\mathrm{In_\textit xGa_{1-\textit{x}}N}$ well area. This process is presented schematically in figure 8(lower panel). It should also be noted that in the system discussed, diffusion of nitrogen vacancies may also occur, which is observed experimentally at the InGaN/GaN interfaces [51–53]. However, due to the experimental values mentioned above of the diffusion coefficient of nitrogen atoms ($D_\mathrm N$), it should be assumed that the flux of diffusing nitrogen vacancies in the analyzed system will be a minority one. Although the experimental value of $D_\mathrm {N}$ in $\mathrm{In_\textit xGa_{1-\textit{x}}N}$ ($x = 0.11, 0.22$) is not known, it should not be significantly different from the known value for GaN, due to the relatively small concentration of In atoms in the alloy $\mathrm{In_\textit xGa_{1-\textit{x}}N}$ considered. Therefore, in our opinion, the described and illustrated process can be adequately explained by the well-known Kirkendall mechanism, which is already used to describe the processes taking place at the nanoscale. The values of the diffusion coefficients $D_\mathrm{Ga}$ and $D_\mathrm{In}$ calculated by us indicate that the main reason for the thermal decomposition of the $\mathrm{In_\textit xGa_{1-\textit{x}}N/ GaN}$ QW system is the diffusion of In atoms. As we showed in our paper [14], the temperature-induced decomposition of the $\mathrm{In_{0.2}Ga_{0.8}N/GaN}$ QW system is a process of several stages, initiated by the formation of 'initial voids', which eventually take the form of 'final voids' with areas where metallic indium precipitates and InGaN amorphous material can be identified (compare figures 3–5 in [14]). The In precipitate areas originate from the diffusion of In atoms that occurs through the dominant flux $J_\mathrm{In}$. Therefore, in our opinion, the described Kirkendall mechanism is the process responsible for the 'initial voids' stage of formation. Based on the calculated temperature dependencies of the diffusion coefficients $D_\mathrm{In}$ and $D_\mathrm{Ga}$, one can estimate the diffusion lengths L of both metal atoms in the alloy $\mathrm{In_{0.22}Ga_{0.78}N}$ at two experimental annealing temperatures applied in [14], i.e. $900\,^{\circ}$C and $930\,^{\circ}$C. Using the relation $L = \sqrt{2Dt}$ where t is the annealing time equal to 30 min. [14], we calculated the following diffusion length ratio $L_\mathrm{In}/L_\mathrm{Ga} = 14$ already for the lowest $900\,^{\circ}$C annealing temperature ($L_\mathrm{In}/L_\mathrm{Ga} = 13$ at $930\,^{\circ}$C). This means that indium atoms can diffuse in the $\mathrm{In_{0.2}Ga_{0.8}N/GaN}$ QWs system distances, which are by at least one order of magnitude longer than those of gallium atoms.

Moreover, the preferred diffusion direction of In atoms should be the direction toward the GaN substrate. This is due to two reasons: first, the GaN substrate is a much bigger reservoir of metal vacancies than the area of the GaN barrier and second, there is only one interface towards the GaN substrate, as opposed to their whole set, corresponding to the five InGaN/GaN QWs. As was shown in [44] there is an embedded internal electric field present at the InN/GaN interface area that influences the motion of the vacancies; moreover, the preferred migration direction of the In atoms, conditioned by the presence of a lower migration energy barrier of ${\sim}1.70$ eV, is just the direction toward the GaN substrate. Based on these two facts, it can be concluded that the thermal decomposition of the $\mathrm{In_{0.2}Ga_{0.8}N/GaN}$ QWs system, induced by the Kirkendall mechanism, should start from the lowest QW, as is observed experimentally [13, 14]. Our proposed explanation of the foundations of thermal decomposition of the $\mathrm{In_\textit xGa_{1-\textit{x}}N/GaN}$ MQWs system is also in good agreement with other experimental observations discussed in [2–7, 11–14].

It should be noted here that the presented mechanism of the Kirkendall effect was disscussed only in the simplest case, that is, we analyzed uncharged defects and one migration direction [0001]. In order to obtain a complete picture of this effect in $\mathrm{In_\textit xGa_{1-\textit{x}}N/GaN}$ systems, one should analyze the influence of various defects' charge states on the values of diffusion coefficients of In and Ga atoms. This issue requires additional time-consuming DFT calculations, which will be the subject of the next paper.