First-principles study of migration mechanisms and diffusion of carbon in GaN

Carbon related defects are readily incorporated in GaN due to its abundance during growth both with MBE and CVD techniques. Employing first-principles calculations we compute the migration barrier of neutral carbon interstitials in the wurtzite GaN crystal. The Minimum Energy Path (MEP) and the migration barriers of these defects are obtained using the Nudged Elastic Band (NEB) method with the climbing image modification (CI-NEB). In addition, the Dimer method is used to verify the results. The results yield a quantitative description of carbon diffusion in the crystal allowing for the determination of the most probable migration paths.

band method (CI-NEB) [21,22]. Using the CI-NEB method we ensure that the calculated energy corresponds to the migration barrier for a given path. In addition, we get information about the atomic configuration in the saddle point.

Background
It is empirically observed that most diffusion processes in solids exhibit a temperature dependence described by an Arrhenius law where D 0 is a temperature independent factor, Q is the activation energy for the atomic jump mechanism and k B T is the product of the Boltzmann constant with temperature. Deviations from this equation may be caused by the existence of short circuiting diffusion paths such as dislocations and grain boundaries, multiple diffusion mechanisms and impurity effects [23].
The jump rate is characterized by an attempt frequency and a thermodynamic factor that dictates the probability of an attempt resulting in a successful jump.
The attempt frequency, Γ 0 , can be obtained within harmonic transition state theory via the Vineyard equation [24], but often the Einstein or Debye frequency is used instead. ΔG represents the isothermal work associated with the motion of the particle from the initial to the final state.

Computational method
The calculations were performed with the vienna ab initio simulation package (VASP) [25] using the projector augmented wave (PAW) method [26,27]. The generalized gradient approximation (GGA) in the parametrization by Perdew, Burke and Ernzerhof (PBE) [28] was used. The investigation of the lowest energy paths was performed using the CI-NEB method [21,22] as well as the Dimer method [29] as implemented in VASP through the VTST-Tools by Henkelman, Jónsson and others [30]. The convergence of the size of the supercell was investigated using supercells of 64, 72, 96, 128 and 144 atoms. The 64-atom supercell is found to be inadequate to contain the deformations induced by the migration of carbon. The 72-atom supercell might be sufficient for certain migrations but not all. Thus, we chose the 96-atom supercell to be the safest and most efficient option at hand. The plane-wave energy cutoff was set at 450 eV and a Γ-centered 2 × 2 × 2 k-point mesh is used for the sampling of the Brillouin zone. The atomic configurations were relaxed until the maximum force per atom was less than 5 × 10 −3 eV/Å. In the case of the NEB calculations, the images were relaxed until the maximum force per atom was no more than 10 −2 eV/Å.
The crystallographic parameters of the w-GaN we obtained from our calculations are a = 3.211Å, c/a = 1.629 and u = 0.377, which are in good agreement with the experimental values a = 3.189Å, c/a = 1.626 and u = 0.377 and consistent with previous DFT calculations [31]. The bandgap at the Γ point is 1.757 eV which underestimates the experimental value of 3.4 eV as expected from GGA.
The main sources of error in DFT calculations of point defects are the electrostatic and elastic interactions of the defect between neighboring supercells and the underestimation of the bandgap. Lany and Zunger [32] describe a number of methods to overcome these shortcomings of DFT. However, migration barriers, the main focus of the present work, are calculated as  energy differences of electronically similar configurations. Thus, there is no need to apply any correction scheme. A typical systematic error in the formation energies calculations is considered to be approximately 0.1 eV [33]. Performing identical migrations at different sites we verified that the uncertainty of the migration barrier also lies within 0.1 eV.
In the NEB method [34], a set of "images" of the system is used to represent the migration path from the initial to the final configuration. The images (i.e., atomic configurations along the migration path) are connected with springs to resemble a string (or band). An optimization algorithm is then applied to relax the string down towards the minimum energy path. In the climbing image modification the highest energy image is driven to the saddle point by neutralizing the forces along the band. Since the image converges to the saddle point, the exact migration barrier can be calculated. In addition, the Dimer method was used to verify the saddle point location and the migration barrier acquired via the CI-NEB.

Determination of diffusion paths
The two most stable configurations for the neutral carbon interstitial are the type-1 and type-2 C-N split complexes shown in Fig. 1. They differ by less than 0.1 eV [35]. Hence, both cases are investigated as potential initial and final configurations for the diffusion of carbon. This small energy difference and the similarity of the two configurations allows for migrations accompanied by transformation from one type to the other. In this case, the migration barrier is defined as the energy difference between the saddle point and the lowest energy configuration. As a result, the barrier in this case could be smaller by ∼ 0.1 eV following the opposite direction of the reaction.
The migration paths investigated in the present work have components both parallel and perpendicular to the c-axis of the wurtzite crystal. The endpoints of the migrations are either a type-1 or type-2 split interstitial. Fig. 2 presents the different migration paths considering jumps among the first and second nearest neighbors. The different combinations of paths and endpoint configurations result in nine different calculations.
Process A, shown in Fig. 2 refers to motion between two out of plane first nearest neighbors. A single gallium atom forms four bonds with nitrogens in the ideal crystal. Three of these nitrogens lie on a plane which is perpendicular to the [0001] direction. Process A describes the migration of a carbon atom forming a split with one of these three nitrogens to the fourth nitrogen which is out of plane.
In the case of process B, migration occurs between two adjacent nitrogens lying on a plane perpendicular to the [0001] direction. This process is a first nearest neighbor migration process and together with process A is expected to exhibit the lowest migration barriers.
As an additional process, we consider process C in which migration occurs between the second nearest neighbors which lie in different planes with respect to the [0001] direction. Even though this process utilizes the hexagonal channel of the wurtzite structure, it is expected to exhibit the highest barriers since it is a second nearest neighbor migration.

Results
The difference of the formation energy between the two types of interstitials is 0.13 eV with the type-2 being the lowest in energy. The bond length between the carbon and the nitrogen is 1.29Å and 1.31Å for the type-1 and type-2 respectively.
The results obtained using the CI-NEB method are summarized in Table 1. The lowest migration barrier is observed in the case of process A from a type-2 interstitial to a type-1. This barrier becomes even smaller by approximately 0.1 eV if the reaction path is reversed. The diffusion potential energy along the reaction path is shown in Fig. 3.
The Dimer method verifies the results acquired by the CI-NEB for all the process A migrations. Both methods result in the same atomic configuration for the saddle point and consequently the same migration barrier. In the case of the lowest barrier migration, the saddle point lies closer to the type-1 interstitial. With respect to the straight line between the initial and final configurations, the saddle point forms an angle of 29 • and 16 • with the type-1 and type-2 ends respectively.
First nearest neighbor in plane migrations (process B) exhibit barriers close to 3 eV making them unfavorable compared to their out of plane counterparts. As expected, second nearest neighbor out of plane migrations require a much greater activation energy.  In order to estimate the temperature where carbon interstitials are activated we assume a jump rate of 1 Hz and a typical Debye frequency of 10 THz. Then, using eq. (2) and the calculated migration barriers, one can estimate the temperature where migration occurs. Hence, using Γ = 1 Hz, Γ 0 = 10 13 Hz and ΔG = 2.3 eV we estimate that carbon interstitials are activated at temperatures close to 890 K. This temperature is relevant during growth of GaN both with MBE and CVD techniques.

Conclusions
We employed the Climbing Image Nudged Elastic Band technique to explore migration barriers of carbon diffusion in GaN. The Dimer method was also used in order to verify the results obtained by the CI-NEB in the case of the three lowest energy barrier migrations. The two most stable interstitial types of carbon are investigated with respect to diffusion via first and second nearest neighbor mechanisms. The first out of plane nearest neighbor mechanism exhibits the lowest migration barriers indicating a favorable direction of migration along the [0001] axis of the crystal. However, for growth temperatures, where carbon migration is relevant, the in plane migration is also competing. Migration via the second nearest neighbor mechanism exhibits high energy barriers minimizing the probability for this mechanism to contribute to the carbon diffusion.