First-principles study of vacancy-assisted impurity diffusion in ZnO

Group-III elements act as donors in ZnO when incorporated on the Zn site. Their incorporation and behavior upon annealing is governed by diffusion, which proceeds mainly through a vacancy-assisted process. We report first-principles calculations for the migration of Al, Ga, and In donors in ZnO, based on density functional theory using a hybrid functional. From the calculated migration barriers and formation energies, we determine diffusion activation energies and estimate annealing temperatures. Impurity-vacancy binding energies and migration barriers decrease from Al to In. Activation energies for vacancy-assisted diffusion are lowest for In and highest for Al.

Zinc oxide has a band gap of 3.4 eV.It can support high levels of n-type doping and still exhibit a high degree of optical transparency in the visible range, making it a technologically attractive transparent conductor. 1Group-III impurities substitute on the Zn site and act as shallow donors.Al, Ga, or In are sometimes unintentionally present, leading to background n-type doping. 2 These elements are also used for intentional doping up to degenerate electron concentrations. 3ompensating vacancies form at high dopant concentrations. 4Controlling dopant concentrations and profiles is essential for most applications, and this requires understanding of the diffusion process.
Here we address the stability and diffusion of Al, Ga, and In in ZnO by calculating formation energies and diffusion activation energies from first principles.The diffusion activation energy Q governs the thermal dependence of the diffusion constant D via the Arrhenius law D = D 0 exp (−Q/k B T), where k B is the Boltzmann constant.
Impurity diffusion is invariably mediated by point defects. 5In the case of n-type ZnO, interstitialassisted processes are unfavorable because of the high formation energy of self-interstitials in donordoped ZnO. 6,7 n addition, Zn interstitials have donor character and are positively charged, which means they would exhibit a repulsive interaction with donor impurities.In contrast, Zn vacancies have low formation energies and are thus likely to be present in n-ZnO. 6Recent experimental studies 8,9 indeed suggest that impurity diffusion is mediated by a defect in the −2 charge state, consistent with the charge state of V Zn in n-ZnO. 6Other experiments have shown that complexes between a Zn vacancy and the substitutional Ga atom, which we denote as V Zn Ga Zn complexes, act as acceptors and can be important for Ga migration and segregation during the growth of the material. 10,11 ompensating V Zn Al Zn (and larger) acceptor complexes have been found in Al-doped ZnO. 12 Previous theoretical studies 13 for the migration of group-III impurities in ZnO have been based on conventional density functional theory (DFT) using the generalized gradient approximation (GGA) in the PW91-parametrization (Perdew-Wang 91). 14Conventional DFT calculations in local or semilocal approximations are known to severely underestimate band gaps, and this shortcoming also affects the position of defect levels in the band gap and thus formation energies, binding energies, and migration barriers. 15Even when the GGA+U approach is used, where on-site Coulomb correlations are included via a Hubbard-like repulsion term, the band gap is only partially corrected.Calculations for the migration of the zinc vacancy using this technique have lead to different results. 6,16 n the present work we devote specific attention to overcoming this problem by performing calculations using a hybrid functional, which allows for accurate description of band gaps, and hence also of formation energies and energy barriers.After a discussion of the vacancy-assisted migration mechanism, activation barriers are calculated, and temperatures at which diffusion becomes possible are estimated.
We perform first-principles hybrid functional calculations within DFT, 17,18 using the projectoraugmented wave method [19][20][21] as implemented in the Vienna Ab initio Simulation Package. 22,23 e use the hybrid functional of Heyd, Scuseria, and Ernzerhof (HSE) 24,25 with a Hartree-Fock mixing parameter of 0.36, which has been shown 26 to produce accurate values for the band gap and structural parameters.The calculated band gap is 3.29 eV (experimental value 27 3.44 eV), and the equilibrium cell parameters are a = 3.25 Å, c/a = 1.60, and u = 0.38, in good agreement with experiment. 27,28 int-defect calculations are performed in a 96-atom supercell with two special offk-points in the Brillouin zone. 29efect formation energies are given by 30 where q is the charge state of the defect, E tot are the total energies obtained from the calculations, and F is the chemical potential of electrons in the solid, i.e., the Fermi level, referenced to the valence-band maximum (VBM).The Fermi level also determines the concentration of free carriers in the conduction band and is in the vicinity of the conduction-band minimum (CBM) under n-type conditions.Total energies of supercells for charged defects are corrected for artificial interaction between periodic images. 31,32 he atomic chemical potentials µ i ≤ 0 account for the atoms added or removed in the defect supercell with respect to the perfect bulk crystal.The µ i are given relative to the energies of the pure phases, µ 0 O = 1 2 E tot (O 2 ), and µ 0 Zn = E tot (Zn bulk), and are subject to the stability condition: H f is the enthalpy of formation, obtained as the total energy difference between the compound and its constituent elements (experimentally, 33 H f = 3.63 eV).The limits are set by µ O = 0 (oxygen-rich) and µ Zn = 0 (zinc-rich).Additionally, the chemical potentials of impurity atoms are limited by the enthalpies of formation of their oxides: Al 2 O 3 , Ga 2 O 3 , and In 2 O 3 .
Table I shows the bond lengths for isolated substitutional Al, Ga, or In donors on a Zn site, in the +1 charge state, as well as for the Zn vacancy in the −2 charge state (V −2 Zn ).These are the relevant charge states under n-type conditions.For Al, there is a −8% relaxation of the surrounding O atoms towards the impurity.For Ga, the relaxation is also inwards but less strong (−4%), while In shows an outward relaxation of +6%.This trend is consistent with increasing atomic size of the group-III impurity. 34he zinc vacancy leads to an anisotropic outward relaxation of the oxygen atoms (+11%, +15%), as shown in Fig. 1.Its formation energy is obtained according to Eq. (1) as This formation energy increases when going towards zinc-rich conditions, and also when lowering the Fermi level.Comparing the formation energy with previous studies (for F = 0 eV and µ O = 0 eV), Ref. 6 reported 4.93 eV for a scaled GGA+U, Ref. 35 shows 6.4 eV for HSE (at 37%), and Ref. 16 reported 3.34 eV for GGA+U and 6.95 eV for HSE (at 37.5%).We now consider the complexes [V Zn Al Zn ] − , [V Zn Ga Zn ] − , and [V Zn In Zn ] − , in which the substitutional impurity can be either in the same zinc plane as the vacancy ("in-plane") or in the neighboring zinc plane with an oxygen plane in between ("out-of-plane").In the latter case the two different stacking sequences of the planes containing the impurity and the vacancy are inequivalent in the   wurtzite structure, see Fig. 1.However, we find the energy difference between the two out-of-plane configurations to be small (less than 0.1 eV in GGA-PBE) (Perdew-Burke-Ernzerhof, Ref. 36).The binding energy of the complex is given by the difference of formation energies, 37 e.g., which is positive when actual binding occurs.The results for binding energies are shown in Table II; the binding energy is lowest for the In-containing complex and increases for Ga and Al.Our values are systematically larger than those reported in Ref. 13.We attribute the differences to the following factors: (1) The authors of Ref. 13 calculated the binding energy with respect to a reference state in which the two species (vacancy and impurity) are present within the same supercell.In contrast, our reference for the binding energy consists of the isolated species calculated in separate supercells [Eq.( 4)], thus eliminating any residual interaction.(2) Their use of the GGA functional, as opposed to the HSE hybrid functional used in our present work.We find that the GGA-PBE gives binding energies that are 0.5 eV smaller that those obtained using the HSE hybrid functional.(3) The neglect of charge-state corrections, which we treat according to Refs. 31 and 32.
Migration barriers E m are calculated as total-energy differences between the stable and the saddle-point configurations.The saddle-point configuration is obtained from a nudged elastic band 38 This article is copyrighted as indicated in the article.Reuse of AIP content is subject to the terms at: http://aplmaterials.  calculation for the exchange process, employing GGA in the PBE 36 parametrization.Afterwards, the geometry is scaled to the HSE equilibrium values, and the forces are minimized using a quasi-Newton method. 39The calculated migration barriers are given in Table II.These barriers correspond to an exchange process of the vacancy with a host cation or the impurity where the atom jumps into the vacancy.Comparing with previous studies for the exchange of the vacancy with the host cation, Ref. 6 obtained 1.4 eV with a small anisotropy of 0.1 eV, and Ref. 16 obtained 1.23 eV for the in-plane and 0.78 eV for the out-of-plane barriers.For the vacancy-impurity exchange processes, we find the same trend as for binding energies: the migration barrier increases in the sequence In < Ga < Al.The barrier for In is lowest because its starting point is already unfavorable (rather than the saddle point being much lower for In than for Ga or Al).The migration barriers obtained by using the HSE hybrid functional are about 0.2-0.4eV higher than those obtained using GGA-PBE (which are consistent with the GGA-PW91 values reported in Ref. 13).
To discuss what actually happens during diffusion, additional considerations are necessary.When a vacancy-impurity complex forms with a positive binding energy, three fundamental processes can occur related to vacancy-assisted migration. 40 II.This process is necessary (because it is the only one that changes the location of the impurity), but not sufficient on its own.For diffusion, exchange needs to be accompanied by either rotation or dissociation.(b) The rotation process: a host Zn atom next to the impurity jumps into the vacancy, leaving the complex bound but with a rotated orientation.(c) The dissociation of the complex due to a jump of a different host atom into the vacancy.For the rotation process, we assume that the migration barrier is the same as the barrier for an isolated Zn vacancy (Table II).For the dissociation process, it is a good approximation to assume the barrier to be the sum of the binding energy and the migration barrier for an isolated Zn vacancy.The barrier for dissociation is therefore always higher than the barrier for rotation.Thus, diffusion will proceed through the exchange+rotation mechanism.
The activation energy for diffusion Q has two contributions (see, e. g., Ref. 40).First, the formation energy of the vacancy next to the impurity; this is given by the formation energy of the vacancy reduced by the binding energy to the impurity (see Table II).Second, the larger of the two barriers for either the exchange or the rotation process; both are necessary for diffusion to take place, but since they are independent processes, the rate will be limited by the higher of the two barriers.Thus, Q depends on the atomic and electronic chemical potentials via the formation energy of the mediating defect and thus both on the sample and the chemical environment.(If anharmonic contributions to the vibrational entropies of the stable and saddle points are substantial, changes to the activation energy beyond harmonic transition state theory are conceivable.)To enable a comparison between the three different impurities, and also to experimental results, we make a specific choice.We assume µ O = −1.35eV, which corresponds to equilibrium with O 2 gas at T = 1173 K and p = 0.200 bar.
For the Fermi level, we assume that F is at the CBM.Under these assumptions, E f (V −2 Zn ) = 1.57eV.Note that under these conditions E f (V −2 Zn ) − E bind ([XV Zn ]) is always positive for all three impurities, i. e., E f ([V Zn X Zn ] − ) > E f (X + Zn ).To estimate a temperature above which diffusion effects become experimentally observable, we determine the temperature at which the transition rate reaches one jump per second and call this the defect annealing temperature T a .The transition rate can be estimated as = 0 exp (−Q/k B T), where 0 is the attempt frequency of the process.0 can be approximated by a typical phonon frequency (10 13 s −1 ), leading to T a /Q = 388 K/eV.Note that a different choice of 0 changes T a only logarithmically.
Table III shows the activation energies and annealing temperatures that we obtain from our calculations.The listed activation energies activate all necessary processes in all directions.(The barrier for diffusion within the c plane is lower by 0.17 eV for Al, and by 0.24 eV for Ga.Indium diffusion is easier along the c-axis by 0.21 eV.)In has the lowest activation energy and thus diffuses the easiest.The activation energy for Ga is considerably higher, by 0.68 eV, which gives a difference in annealing temperatures of 260 K.The activation energy for Al is higher again, by 0.21 eV.
The table also lists experimental results for the diffusion prefactors and activation energies.We note that there are considerable differences between the reported experimental activation energies, which makes it difficult to compare with our results.One reason for the variation could be that different chemical environments can produce different activation energies because of variations in the vacancy formation energy, as described above.However, there is also a large variation in prefactors (over as much as 10 orders of magnitude).An estimate for D 0 can be made from the  III).
uncorrelated random walk, where it is approximately equal to 0 d 2 with d, the jump distance, equal to 3.25 Å in ZnO.This leads to D 0 ≈ 10 −2 cm 2 /s.Experimental results for which D 0 is closer to this value are likely to be more reliable.There is potentially an increase in the prefactor by a factor exp (S/k B ), due to additional entropy differences along the migration path or formation entropy of the assisting vacancy.
The diffusivities for the experimentally obtained parameters are shown in Fig. 3.For the calculated values, the S is varied between 0 and 3 k B , the latter being a conservative estimate for an upper limit.This increases the diffusivity by a factor of about 20, which amounts to a parallel upwards shift in the Arrhenius plot.For Al and Ga, the experimental diffusivities fall into the region of values that is in agreement with our calculations.We note, however, that the diffusivities reported in Ref. 41 are extracted from a model that is not necessarily compatible with the main assumption in the present paper, namely, that Al diffuses via a vacancy-assisted process.Thus, a direct comparison with the values obtained by our calculations is not straightforward.For In, the experimental diffusivities are significantly smaller than the theoretical values.However, the calculations address the situation in which only the impurity and native point defects at their equilibrium concentration are present.Furthermore, additional defect sites in real samples could act as trap sites and significantly slow down diffusion.There are indications that this may indeed be the case for In: Reference 9 introduced In by implantation, a process that may well lead to more damage than in the case of the lighter Ga and Al atoms because of the larger mass of In.Reference 42 used indiffusion, but reported precipitation of a new phase and dislocation formation, which again may provide defect sites that act as traps for the diffusing impurity.
In summary, we have studied the activation energy for the diffusion of group-III impurities in ZnO, mediated by the Zn vacancy.We obtained first-principles values by performing hybrid functional calculations for formation energies and migration barriers.The activation energy itself depends on the electronic and atomic chemical potentials.The numbers indicate that In diffuses the easiest, while Ga and Al become mobile only at higher temperatures.
FIG.1.Ball-and-stick model43 of a vacancy (position at ×) in wurtzite ZnO, showing the neighboring Zn atoms (in grey) and oxygen (in red).The six Zn atoms in the middle plane surrounding the vacancy are all in equivalent positions and are denoted as in-plane neighbors.The Zn atoms in the upper and lower planes (all three equivalent atoms), denoted as out-of-plane neighbors, are inequivalent.
Fig. 2 illustrates them schematically.(a) The exchange process, for which we reported migration barriers in Table

FIG. 3 .
FIG. 3. Arrhenius plot of diffusivities based on the experimentally obtained parameters and on the first-principles calculations (TableIII).

TABLE I .
Calculated bond lengths to neighboring oxygen atoms.The entry for Zn gives the values for the perfect wurtzite crystal.The three planar bonds connect to the nearest oxygen plane, while the axial bond leads to the next-nearest oxygen plane.For the Zn vacancy, the values refer to distances measured from the nominal Zn site and are unequal due to symmetry breaking in the supercell.

TABLE II .
Calculated binding energies E bind of impurity-vacancy complexes, and migration barriers E m for the exchange process.The first line gives the migration barrier for the isolated vacancy.Binding energies and migration barriers are given both for the in-plane and the out-of-plane configurations.
FIG.2.Migration processes for a vacancy next to a substitutional impurity (grey): (a) exchange, (b) rotation, after which the impurity and vacancy are still bound on neighboring sites, (c) dissociation of the complex, where the impurity and vacancy are separated from each other.

TABLE III .
Calculated activation energies and annealing temperatures, and experimental results.