I. INTRODUCTION
For the last few decades, zinc oxide has been under intensive investigation because of its potential applications in electronic, optoelectronic, semiconductor, and spintronic devices.Reference Dietl and Ohno1, Reference Kamilla and Basu2 Particularly, the wurtzite structure composed of alternating planes, which are occupied by tetrahedrally coordinated O2− and Zn2+ ions, largely enhances the piezoelectricity and pyroelectricity properties of ZnO. Furthermore, point defects and impurities during the crystal growing process significantly influence the structural and electronic properties of materials. Therefore, knowledge of migration behaviors of point defects in ZnO material is essential to control the formation of native point defects and explore novel application fields. A large amount of literature in this area has already been published.Reference Oba, Choi, Togo and Tanaka3−Reference Look, Claflin, Alivov and Park6 Oxygen vacancy, zinc interstitial, and zinc Frenkel pairs in ZnO have been studied by electron paramagnetic resonance.Reference Vlasenko and Watkins7, Reference Vlasenko8 A detailed study of the electronic and structural properties and stability of all native point defects in ZnO was performed by first-principles methods.Reference Janotti and Van de Walle9 Erhart and AlbeReference Erhart and Albe10, Reference Erhart and Albe11 determined the migration paths of intrinsic defects in ZnO by using density functional calculations combined with the climbing image nudged elastic band (NEB) method. Defect clusters in ZnO could give valuable insights into the dynamics of microstructural change in materials and point defect incorporation and dissociation, whereas isolated defects show a high tendency to form defect complexes by combining with other defects under radiation damage conditions. However, experimental information about defect cluster migration in ZnO is ambiguously interpretable on the atomic scale. Since the influence of defect clusters on structural properties of ZnO is still controversial partly because of a lack of theoretical studies, computer simulations may provide important insights into the migration behaviors of point defects and defect clusters in ZnO.
The NEBReference Henkelman and Jónsson12 and the dimerReference Henkelman and Jónsson13 methods have been widely used to search for minimum energy paths (MEPs) of the migrations of defects and defect clusters in materials. For a given initial and final state of a transition, the MEP is easily found using the NEB method. However, it is very challenging to find the MEP for the transition, if the final state or the mechanism of transition is unknown. The dimer method is an independent implementation of hybrid eigenvector following approachReference Munro and Wales14 and was developed by Henkelman and Jónsson. Compared to the NEB method, the dimer technique is an especially useful tool for a cluster when the lowest energy final state of its transition is not easily predictable. This method has been used to study possible transition states of defects in SiC and α-Fe,Reference Gao, Henkelman, Weber, Corrales and Jónsson15 the migration behaviors of tungsten clusters on nanoparticle,Reference Chen, Hu, Gao, Deng and Sun16 and Cr-vacancy clusters in α-Fe.Reference Chen, Gao, Hu, Hu, Terentyev, Heinisch, Henager and Khaleel17
In this paper, we have used the dimer method to probe possible transition states, MEPs, and migration mechanisms of vacancies, interstitials, and an interstitial-vacancy pair in zinc oxide. The migrations of zinc and oxygen interstitials occur by a kick-out mechanism and are energetically preferred over the vacancies. Finally, we devote more attention to an interstitial–vacancy pair as a complex consisting of an octahedral oxygen interstitial and a zinc vacancy and explore how a vacancy affects the migration behavior of an IO during the migration processes.
II. DETAILS OF THE CALCULATION
In each configuration, the force calculations were carried out by using version 2.19 of the DL_POLY package.Reference Smith and Forester18 The computational cell of wurtzite structure consisted of 7 × 7 × 7 unit cells containing 1372 atoms. Three-dimensional periodic boundary conditions were applied and calculations were carried out at constant volume. We use the Buckingham potential for ZnO to describe the atomic interaction with parameters previously published by Kulkarni et al.Reference Kulkarni, Zhou and Ke19 This potential has been widely used to simulate defect and monovalent ion incorporation, to produce the wurtzite crystalline lattice parameters, and also to predict elastic and mechanical properties of ZnO nanostructures.Reference Kulkarni, Zhou and Ke19–Reference Dai, Cheong, Sow, Lim and Tan21 The long-range Coulomb force was considered by using the smoothed particle mesh Ewald method with a cutoff distance of 10 Å. The van der Waals cutoff distance chosen in the simulations was 10 Å.
The possible transition states and migration paths of defects are determined by the dimer method. In the dimer method, the saddle point search is completed through moving the dimer on the potential energy surface (PES), which is the 3N-dimensional surface of all possible configurations of the N atoms in the system, from the vicinity of the potential energy minimum of the initial state up to a saddle point. Basins of the PES are stable states, while saddle points correspond to the lowest energy transition points between states, from which the minimum energy barrier for a transition to take place can be determined. Thus, during the search, each time the dimer is displaced; it is also rotated to find the minimum energy configuration. Once the saddle points for the transitions involved in defect migration are determined, it is possible to trace out the MEPs. For each initial configuration, we carried out 100 dimer searches with different initial random vectors that have nonzero components only for a 20–30-atom subset of the system near the defect, which allows the system to be displaced from a distribution of nearly equal potential energy minima in the initial basin. The dimer separation was set at 1 × 10−4 Å, and the values of the finite difference steps for rotation and translation were 10−4 and 10−3 Å, respectively. The maximum move distance was set at 0.15 Å and the search was stopped when the maximum force of atoms became less than 0.01 eV/Å.
III. RESULTS AND DISCUSSION
A. Migration of zinc vacancy
First-principles calculationsReference Janotti and Van de Walle22, Reference Wardle, Goss and Briddon23 have revealed that the migration of zinc vacancy occurs when the nearest neighbor atom directly jumps into the vacancy site. Here, we show two typical migration processes of a zinc vacancy that may occur in the wurtzite-structured ZnO.
The MEPs of a single zinc vacancy (VZn) using the dimer method and the corresponding processes are shown in Fig. 1. As expected, the common migration mechanism is the nearest neighbor Zn atom jumping into the vacancy position with the energy barrier of 1.77 eV (a→b) and 1.35 eV (a→b′), respectively. For path 1, the nearest neighbor zinc atom in the same basal plane as the VZn jumps into the vacancy position, which corresponds to the position exchange of the Zn atom and vacancy (see a→b in Fig. 1). It is obvious that the energy curve of path 1 is symmetric and the middle point of the migration path is a saddle point. Path 2 occurs when a nearest neighbor zinc atom in the adjacent basal plane hops into the vacancy site (see a→b′ in Fig. 1). The energy barrier for this path is 1.35 eV and is about 0.43 eV lower than that for path 1. Taking into account the large deformation due to the Zn atom movement following path 1, it is expected that path 1 has a relatively high barrier. Both migration energies calculated for the VZn are reasonably consistent with first-principles calculations reported by Erhart and Albe.Reference Erhart and Albe10 These two migration processes demonstrate that the diffusion of a zinc vacancy in ZnO is anisotropic.
Migration paths and calculated energy paths that correspond to this migration process of a zinc vacancy. (a) Starting configuration, (b) and (b′) two possible final configurations, and (c) calculated energy barriers corresponding to path 1 (a→b) and path 2 (a→b′). The red and light purple spheres represent the O atoms and Zn atoms, respectively. Also the dotted light purple circle indicates the zinc vacancy. Upper panels show top views, while lower panels show side views. Side views are perpendicular to the 
direction and top views are along the 
direction. In the first frame, the respective paths are labeled. The arrow indicates the vacancy movement at the next step.
B. Migration of oxygen vacancy
Similarly, the migration of an oxygen vacancy (VO) occurs by a nearest neighbor O atom directly jumping to the vacancy site. The migration paths determined by the dimer search and the corresponding MEPs are shown in Fig. 2. The nearest neighbor oxygen atom in the same basal plane of VO can directly move into the vacancy position (see a→b in Fig. 2), which leads to an energy barrier of 1.78 eV, similar to that of the VZn discussed above. Path 2 represents a transition of a nearest neighbor oxygen atom in the adjacent basal plane jumping to the vacancy site (see a→b′ in Fig. 2). This transition leads to a lower energy barrier of 1.36 eV, which is in agreement with the calculated value reported by Erhart and Albe.Reference Erhart and Albe11 In spite of the local geometry around the oxygen vacancy with tetrahedral symmetry, the migration of VO also shows an anisotropic behavior based on our calculations, which is similar to the anisotropic migration behavior of the VZn, as discussed above. It should be pointed out that vacancy migration is unlikely to be a dominant migration process in ZnO, considering its higher energy barrier, and hence, a complex defect needs to be studied.
Migration paths and calculated energy paths of an oxygen vacancy. (a) Starting configuration, (b) and (b′) two possible final configurations, and (c) calculated energy barriers corresponding to path 1 (a→b) and path 2 (a→b′). The dotted red circle indicates the oxygen vacancy.
C. Migration of zinc interstitial
A zinc interstitial (IZn) is initially placed in an octahedral interstitial position in the wurtzite lattice, and the possible migration process is described in Fig. 3. The migration path shows that the IZn pushes out the nearest neighbor zinc atom, which then transfers to a new interstitial atom (see a→b→c in Fig. 3). In the case of a kick-out (or interstitialcy) mechanism, the dimer method results in an energy barrier of 0.76 eV as shown in Fig. 3(d). The midway along the migration path between initial and final states is the saddle point. At the saddle point, the two Zn atoms share a lattice site forming a rather symmetric configuration as shown in Fig. 3(b). Clearly, the Zn interstitial is quite mobile via overcoming a small energy barrier, which is 0.59 eV lower than that for the VZn migration (1.35 eV). Taking into account the low energy barrier of zinc interstitial, zinc interstitials may be trapped by other defects or impurities and diffuse with these defects together, which was proposed in Ref. Reference Janotti and Van de Walle22.
The kick-out migration mechanism and the MEP of a zinc interstitial. (a) Starting configuration, (b) saddle point, (c) final configuration where a new zinc interstitial is formed in an adjacent basal plane, and (d) the calculated energy barrier. To easily follow the interstitial movements, we use the blue sphere to represent the initial zinc interstitial.
D. Migration of oxygen interstitial
Of the 100 dimer searches for an oxygen interstitial (IO) at the octahedral interstitial site, 12 searches converge to a saddle point with an energy of 0.13 eV, 13 searches converge to other saddle points with the energies of 0.58–0.96 eV, 2 searches fail to find saddle points, some of 73 searches could not converge to the original initial state, and the others lead to physically improbable energy values. As shown in Fig. 4, three typical migration processes occur, one involving a slight displacement for the transform between two nonsymmetric configurations (see a→c in Fig. 4) and two being the kick-out mechanisms (see a→b′→c′ and a→b″→c″ in Fig. 4), respectively. In Fig. 4(d), we show the MEPs corresponding to these migration processes. After the lattice relaxation of the starting configuration, the oxygen atom at the first nearest neighbor lattice site moves away from its original lattice position about 0.95 Å along the 
direction. As illustrated in Fig. 4(a), the oxygen atom is noticeably pushed upwards due to the existence of the IO. The similar nonsymmetric configuration of an oxygen interstitial was also found in the previous study.Reference Huang, Wang and Wang24 The IO slightly moves along the 
direction and then pushed another oxygen atom at the lattice site upwards (see a→c in Fig. 4), which is easily understood from the side views in Figs. 4(a) and 4(c). The calculated migration barrier for path 1 (a→c) is 0.13 eV. However, the path 1 of the IO only involves configuration transformation and is unlikely to contribute to the long distance migration behavior.
Migration paths and the corresponding energy paths of an oxygen interstitial. (a) Starting configuration, (b′) and (b′′) saddle points, (c), (c′), and (c′′) final configurations, and (d) the calculated energy barriers along path 1 (a→c), path 2 (a→b′→c′), and path 3(a→b′′→c′′). The transition to the nonsymmetric configuration requires the lowest energy barrier. The green sphere represents the initial oxygen interstitial.
Both path 2 (a→b′→c′) and path 3 (a→b′′→c′′) display a similar kick-out mechanism, where the IO moves close to a nearest neighbor oxygen atom, occupies the lattice site, and then pushes it toward the adjacent basal plane to form a new octahedral interstitial. The migration energy barrier for path 2 is determined to be 0.58 eV, which is comparable to that of the IZn migration in ZnO. Furthermore, the final state [see Fig. 4(c′)] is also a nonsymmetric configuration in which the original IO is pushed upwards around 0.25 Å along the 
direction. In the path 3, the IO migrates to an adjacent octahedral site with an energy barrier of 0.96 eV. The energy difference between these paths is about 0.38 eV. A main reason for this is due to the different geometries of the saddle states, leading to quite different energies for migration.
E. Migration of interstitial–vacancy pair
As mentioned above, the low energy barrier of an oxygen interstitial allows it to diffuse at a low temperature. Here, we investigate the migration of an interstitial–vacancy pair (IO–VZn) consisting of an octahedral oxygen interstitial and a single zinc vacancy, to explore how a vacancy affects the migration behavior of an IO. Five possible migration paths and saddle points are determined by running 100 dimer searches. Figure 5 shows the first step of four migration paths and the MEPs for the IO–VZn pair. As demonstrated in Fig. 5, the migration of the IO–VZn pair can be considered as either IO migration or zinc atoms hopping to and from the vacancy. A nearest neighbor zinc atom in the adjacent basal plane moves into the vacancy site (see a→b→c in Fig. 5), which needs to overcome an energy barrier of 0.65 eV. Compared to the out of plane motion of the VZn migration discussed in Sec. III. A, the energy path of path 1 (a→b→c) is distorted, which is affected by the nearby IO. The energy barrier largely decreases from 1.35 to 0.65 eV, and the energy difference between the initial and final states in Fig. 5(d) is about of 0.33 eV. Also the nearest neighbor zinc atom in the same basal plane of VZn could jump into the vacancy position (see a→b′′→c′′ in Fig. 5). The calculated migration energy for path 3 (a→b′′→c′′) is only 1.05 eV, which is 0.72 eV lower than that for a single VZn transition. It indicates that the mobility of VZn could be enhanced in the presence of the IO.
The first step of four possible migration processes and the MEPs for an IO–VZn pair consisting of an O interstitial (green atom) and a Zn vacancy (the dotted light purple circle). Path 1 (a→b→c) and path 3 (a→b′′→c′′) show the migrations of the VZn in the IO–VZn pair. Path 2 (a→b′→c′) and path 4 (a→b′′′→c′′′) illustrate the migrations of the IO in the IO–VZn pair. (d) The calculated energy paths along path 1, path 2, path 3, path 4, and path 5. Also, path 5 in (d) corresponds to the energy barrier for the first step of the multistep migration process in Fig. 6. Upper panels show top views, while lower panels show side views. In the first frame, the respective paths are labeled.
Compared to the mobility of oxygen interstitial discussed in Sec. III. D, it is unexpected for the IO in the IO–VZn pair to migrate rather slowly. The kick-out mechanism is still a dominant mechanism for the migration of the IO in the IO–VZn pair. The IO replaces the nearest neighbor oxygen atom and then pushes it to form a new octahedral interstitial (see a→b′→c′ and a→b′′′→c′′′ in Fig. 5). Due to the anisotropic diffusion nearby the VZn, the energy barriers for path 2 (a→b′→c′) and path 4 (a→b′′′→c′′′) are 1.38 and 1.14 eV, respectively. Compared to a single oxygen interstitial, the existence of zinc vacancy restricts the IO migration in the IO–VZn pair.
Figure 6 provides an example to illustrate a possible multistep migration process. The final state of the minimum energy path is used as the initial atomic configuration for subsequent dimer searches and this process is repeated until a migration process of the cluster completes. There are two steps needed to complete a migration process of the IO–VZn pair and the corresponding energy barrier is illustrated in Fig. 6. In this migration process, the nearest neighbor oxygen atom migrates into an adjacent octahedral interstitial via the kick-out mechanism (see a→b→c in Fig. 6) with an energy barrier of 0.77 eV [see a→b→c in Fig. 6(f)], and then, the out-of-plane Zn atom hops into the vacancy position (see c→d→e in Fig. 6) overcoming an energy barrier of 1.07 eV [see c→d→e in Fig. 6(f)]. Briefly, the IO–VZn pair needs to overcome an energy barrier of 1.07 eV following this path. The set of movements separates the IO–VZn pair a distance of almost 5 Å.
The multistep migration process of an IO–VZn pair. Lower right panel shows the migration energy path of the IO–VZn pair, where the migration process from c→e leads to the highest energy barrier. The black arrows in (a) and (c) indicate the defect movements in the next frame.
It is of interest to further consider the migration behavior of zinc vacancy–interstitial Frenkel pair consisting of a VZn and IZn with different separation distances, which will be discussed in details below. After the atomic relaxation of the starting configuration, we find that the pair is recombined even for a long distance of 3.8 Å between VZn and IZn. The long-range Coulombic interaction is responsible for the strong attraction between the vacancy and the interstitial. We also performed a simulation when the distance between the vacancy and the interstitial is initially at least 3.2 Å, finding the migration of an oxygen Frenkel pair gives the result similar to the case of zinc Frenkel pair discussed above.
IV. CONCLUSIONS
We have performed a detailed investigation of the migration paths and determined the corresponding MEPs of vacancies, interstitials, and an interstitial–vacancy pair in zinc oxide using the dimer method. A difference between the in-plane and out-of-plane migrations for zinc and oxygen vacancies is as large as 0.5 eV. Based on the calculated migration barriers for zinc and oxygen vacancies, we predict that the migration of vacancies in ZnO is anisotropic. Zinc and oxygen interstitials migrate via the kick out mechanism with rather low energy barriers, which is in agreement with experimental observations and first-principles study. In addition, an interstitial–vacancy pair containing an octahedral oxygen interstitial and a single zinc vacancy exhibits more complex migration behaviors. In the case of an interstitial–vacancy pair, the migration barrier of the zinc vacancy is 0.72 eV lower than that for a single VZn transition, whereas an oxygen interstitial diffuses more difficult, as compared to a single IO transition. Obviously, the influence of other defects on the migration behaviors of defect cluster is notable. Finally, a possible migration path of the IO–VZn pair indicates that the IO or VZn successive jumps can lead to the defect clusters migration, which can provide a guide to probe the migration behaviors of defect clusters under nonequilibrium conditions.
ACKNOWLEDGMENT
This research was supported by the National Natural Science Foundation under Contract No. 11147165. F. Gao was supported from the Division of Materials Sciences and Engineering, Office of Basic Energy Sciences, US Department of Energy under Contract DE-AC05-76RL01830.
