The Interplay of Interstitial and Substitutional Copper in Zinc Oxide

Cu impurities are reported to have significant effects on the electrical and optical properties of bulk ZnO. In this work, we study the defect properties of Cu in ZnO using hybrid quantum mechanical/molecular mechanical (QM/MM)–embedded cluster calculations based on a multi-region approach that allows us to model defects at the true dilute limit, with polarization effects described in an accurate and consistent manner. We compute the electronic structure, energetics, and geometries of Cu impurities, including substitutional and interstitial configurations, and analyze their effects on the electronic structure. Under ambient conditions, CuZn is the dominant defect in the d9 state and remains electronically passive. We find that, however, as we approach typical vacuum conditions, the interstitial Cu defect becomes significant and can act as an electron trap.


INTRODUCTION
ZnO, with a wide bandgap of 3.44 eV (Reynolds et al., 1999), is one of the most widely studied transparent semiconductors, with applications in solar cells (Keis et al., 2002), light-emitting diodes (Tsukazaki et al., 2005a;Tsukazaki et al., 2005b), photocatalysts (Barick et al., 2010), and piezoelectric devices (Wang and Song, 2006). To optimize its performance in applications, it is important to understand and control the properties of impurities in ZnO. The Cu/ZnO system is a very important industrial methanol catalyst (Waugh, 1992;Baltes et al., 2008). Substitution of copper into ZnO (Cu Zn q ) (here and elsewhere in this article, we denote the effective charge of the defect, q, with respect to the lattice site explicitly with the superscript) has been reported to improve the photocatalytic activity (Mohan et al., 2012), ferromagnetism (Xing et al., 2008), and gas sensitivity (Gong et al., 2006) of ZnO. By admittance spectroscopy experiments, Cu Zn q is found to possess an (0/ −) acceptor level at 0.17 eV below the bottom of the conduction band (CBM) in ZnO (McCluskey et al., 2015). Doping with copper has been proposed as a route for producing stable p-type ZnO (Özgür et al., 2005), which, however, has not been successful to date and which partly inspired our current investigation. Moreover, although there are many investigations both experimental and computational on the effects of CuZn on the electrical and optical properties of ZnO, information about Cu interstitials (Cu i q ) in ZnO is limited. In this article, we report the properties of Cu in both substitutional and interstitial forms in ZnO using a hybrid quantum mechanical/molecular mechanical (QM/MM)-embedded cluster approach. For the description of point defects in crystals, the commonly used implementation of density functional theory (DFT) with periodic boundary conditions suffers from finite-size effects (Freysoldt et al., 2014). In contrast, the QM/MM method defines the vacuum reference level unambiguously and describes accurately the short-and long-range polarization effects of a charged defect in a host material. The approach allows us to compare the energetics of different defect configurations and charge states on an equal footing. We focus here on isolated defects, but copper impurities may also form complexes in a variety of ways, an investigation into which is underway and will be reported elsewhere.
To reduce the computational load, we have removed f functions from the oxygen basis set and some of the highly diffuse functions from the cation basis sets, which do not contribute to the bonding in these ionic solids.
For electron exchange and correlation, we have employed the BB1k functional (Zhao et al., 2004), which has been fitted to both thermochemical and kinetic data including 42% exact exchange, and the PBE0 functional (Adamo and Barone, 1999), which is frequently used in plane-wave basis calculations including 25% exact exchange for comparison. In order to embed the QM cluster within a polar environment, the MM region containing 10,460 atoms is treated with a previously derived interatomic potential (Catlow et al., 2008).
The hybrid QM/MM-embedded cluster approach used is implemented in the ChemShell (Sherwood et al., 2003) package. The QM/MM energy is obtained in an additive approach as a sum of QM and MM terms with the interaction energy between the two regions accounted for the QM term whose Hamiltonian includes the embedding potential. The GAMESS-UK (Guest et al., 2005) code is employed in the QM calculations, while the GULP package has been used to calculate the MM contributions.
The formation energies of the Cu substitution with charge q are calculated according to the following reactions: under Zn-rich/O-poor condition and under Zn-poor/O-rich condition.
The formation energies of the Cu interstitial are calculated according to the following reaction: under Zn-rich/O-poor condition and Under Zn-poor/O-rich condition. The chemical potentials of O 2 molecular and single Zn atoms are calculated using GAMESS-UK with the corresponding basis set and density functional; the standard state energy of ZnO is derived from the experimental heat of formation (Haynes, 2014).

RESULTS AND DISCUSSION
We first present the local geometries as well as the defect formation of the Cu interstitial Cu i q and Cu substitutional Cu Zn q ; we next calculate the self-consistent Fermi energies and charge carrier and defect concentrations of Cu-doped ZnO. Finally, we report on the equilibrium between Cu Zn q and Cu i q .

Positions of Cu
As noted, Cu impurities have been widely studied due to their possible influence on the optical properties of ZnO. Moreover, by using the emission channeling technique, 60-70% of the Cu atoms are found to occupy the substitutional Zn site with root-mean-square displacements from the site of 0.16-0.17 Å (Wahl et al., 2004). In our calculation, for the singly negatively charged state of Cu Zn − , the four nearest O neighbors of Cu are displaced outward by 0.11-0.14 Å (using the BB1k functional) as shown in Figure 1A. After an electron is removed from this negatively charged system, the neutrally charged state Cu Zn is formed. As shown by the spin density in Figure 1B, the resulting d 9 Cu impurity drives a Jahn-Teller distortion, with the axial neighbor O relaxing inward by 0.06 Å and the other three non-axial O ions relaxing outward by 0.03 Å. In the +1 charge state, four neighbor O ions all relax inward by 0.04-0.07 Å around the d 10 Cu ion ( Figure 1C).
The Cu impurities could also be present as interstitials in ZnO, in two possible positions: octahedral and tetrahedral sites. As discussed in previous studies by Janotti and Van de Walle (Janotti and Van de Walle, 2007) and Sokol et al. (Sokol et al., 2007), the Zn interstitial is expected to be more stable at the octahedral site than at the tetrahedral site. Hence, here, we only consider the interstitial at the octahedral site.
The calculated configurations are illustrated in Figure 2B. The Cu + interstitial has a lower coordination by electron-rich O 2ions and forms a trigonal pyramid with the closest Cu i -O separation distance of 1.98 Å and the two other distances of 2.01 and 2.05 Å (BB1k structures are shown in Figure 2A). On ionization, this nearly symmetric configuration is broken, with the Cu 2+ ion moving toward one of the lattice oxygens (1.97, 1.98, and 2.11 Å). The next nearest O ions move now toward the interstitial Cu (by 0.16, 0.50, and 0.38 Å) but do not approach close enough to coordinate to this ion directly by a dative bond.

Formation Energies
The calculated formation energies of both Cu Zn and Cu i are plotted in Figure 3. The (0/−) transition level of Cu Zn q is found to Frontiers in Chemistry | www.frontiersin.org December 2021 | Volume 9 | Article 780935 2 lie at 3.54 eV (PBE0) above the valence band maximum (VBM), which is in agreement with previous calculations by Lany and Zunger (Lany and Zunger, 2009), who reported 3.46 eV using generalized gradient approximation (GGA)+U with an additional hole-state correction for the Cu d state, and close to the calculations by Lyons et al. (Lyons et al., 2017), who reported 3.27 eV using the Heyd-Scuseria-Ernzerhof (HSE) hybrid functional. Our results contrast with Yan et al. (Yan et al., 2006), who reported 0.7 eV using local density approximation (LDA), and Gallino and Valentin (Gallino and Di Valentin, 2011), who reported 2.48 eV using B3LYP. The computed ε(0/−) using BB1k is, however, 4.40 eV. The (+/0) transition level of Cu Zn q is found to lie at 1.14 eV (BB1k) and 1.07 eV (PBE0) above the VBM, which yields a deep donor level, which is shallower than that of Lany and Zunger (Lany and Zunger, 2009), who reported 0.37 eV, and Lyons et al. (Lyons et al., 2017), who reported 0.46 eV.
In O-poor conditions, we observe that the Cu i q is more stable than Cu Zn q using BB1k until the Fermi level is very close to the CB from the calculated formation energies. In O-rich conditions, when the Fermi level is near the VBM, the interstitial Cu is still  Frontiers in Chemistry | www.frontiersin.org December 2021 | Volume 9 | Article 780935 3 more stable than substitutional Cu. While for PBE0 results, the substitutional Cu becomes more stable than the interstitial Cu for the Fermi level greater than 2.56 eV under the Zn-rich condition. Under O-rich conditions, the substitutional Cu is the most stable defect type in the d 9 state as a donor.

Charge Carrier and Defect Concentrations
From the computed formation energies, the self-consistent Fermi energy and equilibrium defect and carrier concentrations can be determined. Here, we use the code "SC-FERMI" (Buckeridge, 2019). We focus on the results obtained using the BB1k functional, which reproduces the localization of holes on the oxygen sublattice more accurately than that on other functionals we have tested.
The concentration of each defect X in each charge state q is given by: where N X is the density of sites in which the defect may form, g X q is the degeneracy of the charge state, E F is the self-consistent Fermi energy, and k is the Boltzmann constant.
The electron (n 0 ) and hole (p 0 ) carrier concentrations can be determined by integrating the density of states weighed by the appropriate Fermi-Dirac function: where f e (E) [exp ((E F -E)/kT) +1] −1 is the Fermi-Dirac distribution function and f h (E) 1-f e (E).
The computed self-consistent E F and equilibrium carrier and equilibrium concentrations of Cu impurities with native defects in ZnO as a function of T are shown in Figure 4. The range of temperatures is 0-1500 K, which encompasses common synthesis In O-rich conditions, the E F remains deep in the bandgap, between 1.9 and 2.2 eV above the VBM as T is increased, as shown in the inset of Figure 4A. The carrier concentrations remain below 10 16 cm −3 for T ≤ 1500 K, with the Cu interstitial concentration [Cu i ] three-order of magnitude below. The [Cu Zn ] in the neutrally charged state is the dominant defect in this range of E F, with the concentration above 10 18 cm −3 for T > 600 K, which is close to the experimental result of ∼10 18 cm −3 at room temperature (Kanai, 1991).  In O-poor conditions, due to the lower formation energies, E F moves closer to the CB and even above the CBM as shown in the inset of Figure 4B. From our analysis, ZnO is found to be n-type with electron concentrations n 0 of 10 16 cm −3 for T > 453 K (Hou, 2021) (details of the properties of native defects in ZnO will be published in the future).
We next investigated the equilibrium defect concentrations of Cu impurities in ZnO with fixed E F as a function of T (Figures 5, 6). We mainly considered four conditions of E F at: (A) 0.1 eV above the VBM, (B) 0.1 eV below the CBM, (C) 0.1 eV above the CBM, and (D) 1 eV above the CBM. We note here that the concentrations are not available when the formation energies of the defects are negative based on Eq. 7. In O-rich conditions, the [Cu Zn ] in the neutrally charged state is the dominant defect for all four conditions. In O-poor conditions, when the E F is near to the CBM, the [Cu i + ] is close to [Cu Zn 0 ], but both are in relatively low concentrations.

O Partial Pressure Variation
To compare our results to the experiment, it is important to relate the theoretically defined O-rich and O-poor condition to the oxygen chemical potential under different temperature and partial pressure conditions. The chemical potential of oxygen gas at varying oxygen partial pressures at a given temperature is expressed by: By setting the zero state of μ O (T, p) to be the total energy of oxygen at T 0 K, which is μ O (0, p 0 ) 1/2E total O2 0, the temperature dependence of the oxygen chemical potential at a constant oxygen pressure p 0 is defined as: where H is the enthalpy, and S is the entropy. Based on the data from thermochemical tables (Stull and Prophet, 1971), μ O (T, p 0 )  Frontiers in Chemistry | www.frontiersin.org December 2021 | Volume 9 | Article 780935 6 at p 0 1 atm was calculated by Taylor et al. (Taylor et al., 2016) and Reuter and Scheffler (Reuter and Scheffler, 2001) Our procedure then follows the method by Reuter and Scheffler (Reuter and Scheffler, 2001). The formation energies of Cu impurities in ZnO as a function of the O partial pressure from 10 -22 to 1 atm at 300 and 1000 K are shown in Figure 7. At 300 K, the Cu Zn in the neutrally charged state remains as the dominant defect for O partial pressures from 10 -22 to 1 atm, which is consistent with the defect concentration results. At high temperatures of 1000 K, the Cu i + becomes the dominant defect under very low O partial pressures, with the crossover from Cu Zn 0 to Cu i + at 6.5*10 -9 atm.

Balance Between Cu Substitutional and Interstitial
Bulk ZnO contains significant concentrations of intrinsic defects, to which are attributed the intrinsic n-type conductivity in ZnO (Harrison, 1954;Thomas, 1957;Hagemark, 1976). In the presence of the Cu dopants, there can be interchange of electrons or holes between Cu substitutional with extrinsic defects and Cu interstitial with intrinsic defect Zn vacancy or interstitial, for e.g., The corresponding processes and their reaction energies ΔE f (in eV) are listed in Table 1, with the electron in the CB and the hole in the VB. In general, Cu Zn 0 remains energetically preferable. Cu will not migrate directly from the substitutional to the interstitial site under ambient conditions, owing to the high reaction energy at 6.21 and 6.15 eV. However, under O-poor/Zn-rich conditions, in the presence of the native defect Zn interstitial, Cu will spontaneously migrate to the interstitial site with a reaction energy of −0.45 eV; the Cu i can trap electrons from CB.

Summary and Conclusion
We have investigated the copper dopants in both substitutional and interstitial forms in ZnO from embedded cluster calculations. By computing defect formation energies, we find that the Cu substitutional in the neutrally charged state is the dominant defect under O-rich conditions, which acts as a deep donor, while under Zn-rich conditions Cu interstitial becomes more stable than the Cu substitutional, which is consistent with the results of the equilibrium carrier and effect concentrations as a function of temperature. The Cu will not migrate directly from the substitutional to the interstitial site under ambient conditions, but under Zn-rich conditions Cu will spontaneously migrate to the interstitial site and trap electrons in the presence of Zn interstitial.

DATA AVAILABILITY STATEMENT
The original contributions presented in the study are included in the article/Supplementary Material; further inquiries can be directed to the corresponding authors.

AUTHOR CONTRIBUTIONS
All authors contributed to the components of the science and techniques used in the work and to the development and revision of the manuscript.