Role of hydrogen co-doping on opto-electronic behaviors of Na-H co-doped zinc oxide: a first principle study

In this work, the electronic structure and optical properties are investigated within the framework of the density functional theory (DFT) for different Na-H co-doping scenarios to find out the suitability of H co-doping technique for achieving p-type conductivity in ZnO. Very low formation energies were found for the H co-doped systems compared to others which suggests that they can suppress other n-type impurities and increase the effect of p-type NaZn defects in the lattice. From the electronic structure calculations, we have found that NaZn doped structures with 50% H co-doping produces the best p-type behavior indicating importance of controlling annealing time. Moreover, from the optical calculations, it has been found that NaZn creates impurity states 174 meV above the valence band and electron concentration in these states can be controlled by H co-doping concentration. H co-doping has not produced any substantial lattice strain as compared to other dopants and structures with Na-H co-doping is transparent in the visible light range.


Introduction
A wide range of useful properties displayed by ZnO have made this wide band gap semiconductor a research attraction in recent years [1,2]. ZnO shows a band gap of 3.37 eV at room temperature and a large exciton binding energy of 60 meV [3]. This enables applications in optoelectronics in the blue or UV region, including light-emitting diodes, laser diodes and photodetectors [4]. Due to a strong luminescence in the green-white region of the spectrum, ZnO is also a suitable material for phosphor applications. One of the most striking features of ZnO as a semiconductor is that large area single crystals are available, and epi-ready substrates are now commercialized.
Bulk crystals can be grown with a variety of techniques, including hydrothermal growth [5], vapor-phase transport [6] and pressurized melt growth [7]. Growth of thin films can be accomplished using chemical vapor deposition (MOCVD) [8], molecular-beam epitaxy [9], RF magnetron sputtering [10] and spray pyrolysis [11]. The epitaxial growth of ZnO on native substrates can potentially lead to high quality thin films with reduced concentrations of extended defects. Despite the many attractive applications of ZnO, obtaining p-type doping has proved to be a very difficult task [12].
The cause of this difficulty is the unintentional incorporation of impurities that act as shallow donors, such as hydrogen which is present in almost all growth and processing environments [13]. By means of densityfunctional calculations it has been shown that interstitial H forms a strong bond with O in ZnO and acts as a shallow donor, contrary to the amphoteric behavior of interstitial H in conventional semiconductors [14]. This shallow donor suppresses acceptor behavior generated by any p-type dopants specially by substitutionally doped Na in our case. This suppression of p-type behavior by H is termed as H compensation. Interestingly, hydrogen co-doped with Na into ZnO films can be the solution to the p-type doping difficulty as shown by S S Lin et al [15]. They showed Na-H co-doped sample shows reduced resistivity and p-type conductivity after annealing at 550°C. They proposed a model involving Na Zn , Na i , and H i to illustrate the mechanism of p-type behavior. In this qualitative model, the formation of compensating interstitials is severely suppressed, and the acceptor Any further distribution of this work must maintain attribution to the author(s) and the title of the work, journal citation and DOI. solubility is greatly enhanced by forming H-acceptor complexes during the growth process. The low Na Zn -H i dissociation energy means that it can be easily diffused out by annealing and thus achieving higher concentration of Na Zn resulting p-type conductivity. Many researchers have investigated behavior of native point defects [16], effect of group-I doping and hydrogenation [17], role of Na concentration on photoelectric properties of ZnO [18]. However, the effect of H co-doping with Na in ZnO on the electronic structure and optical behaviors is yet to investigated. Thus, in this work we have tried to investigate the effects of H co-doped with Na using first principle density functional method.
Our first-principles calculations are based on density functional theory (DFT) within the generalized gradient approximation (GGA) and the pseudopotential-plane-wave method. The impurity studies were carried out in supercells (3 × 3 × 2) containing 72 atoms for the Wurtzite structure as showed in figure 1(a). Eun-Cheol Lee et al have showed that in the regime of high dopant concentrations, the dominant defect is the dopant-H complex [19]. As we are modelling a system with high doping concentration the Na-H complex has been taken as the most stable complex in our models. U Wahl et al have shown that Na i prefers the octahedral interstitial sites rather than the tetrahedral interstitial sites [20]. Moreover, in the Na Zn -H i defect complex, the antibonding site of the O atom, which forms three Na Zn -O bonds off the c-axis, is the most stable position [17]. To investigate the overall effect of H co-doping it is compared with substitutional (Na Zn ) and interstitial (Na i ) doping (figures 2(a)-(c)). To understand the compensation process of H in Na Zn doped ZnO, the supercell was doped with two Na atoms in substitutional positions amounting about 5.56% doping. It is compared with 50% compensated (2Na Zn -H i ) and 100% compensated (2Na Zn -2H i ) models (figures 2(d)-(f)). Finally, to investigate the effect of p-type and n-type condition on H co-doping the 50% compensated (2Na Zn -H i ) model is compared with H compensated supercell also doped at the interstitial position (Na Zn -Na i -H i ) (figures 2(e), (h)). We have calculated total ground state energy for all possible configurations (except the symmetric ones) for multi-doped systems and have chosen the configuration with the lowest energy for further calculations. We expect that this work will provide a theoretical basis for understating the H co-doping technique for fabrication of p-type ZnO semiconductor thus broadening its current application space.

Computational methods
Here, the calculations performed (i.e. scf, nscf and optical) were carried out using the Cambridge Serial Total Energy Package (CASTEP) code [21]. Plane wave basis set with pseudo-potential approach based on DFT was used as implemented in the code [22,23]. Accuracy of a DFT calculation mainly depends on choice of exchangecorrelation functionals used for a model. In our study we used Perdew-Burke-Ernzerhof (PBE) with generalized gradient approximation (GGA) parameterization scheme [24]. As GGA approach greatly underestimates the band gap and over or underestimates other properties [25], we used GGA+U method implemented via spin polarized, using formal spin as initial conditions with oxidation state of each element for charge neutralization. Hubbard (U) parameter of 10.5 eV and 7 eV were used for d-orbital of Zn and p-orbital of O, respectively. The interactions between ions and electrons are represented with Vanderbilt-type ultrasoft pseudopotential [26] for all atoms. All calculations including electronic structures (band structure and DOS) and optical were performed with a plane-wave cut off energy of 370 eV with a SCF tolerance of 1×10 −6 eV/atom. For the sampling of the Brillouin zone,´4 4 4k-point grids were generated according to the Monkhorst-Pack scheme [27]. Optical properties due to electronic transitions were calculated with polycrystalline geometry using a smearing value of  0.5 eV and drude damping of 0.05 eV. The standard Broyden-Fletcher-Goldfarb-Shanno algorithm (BFGS) was used for the variable cell optimization where the self-consistence convergence tolerance for the total energy is 1.4×10 −5 eV/atom, the maximum force of the atom is 0.03 eV/Å, the maximum stress is 0.05 GPa using code as implemented in Quantum espresso (pwSCF) [28]. These parameters are found to be sufficient for good convergence of total energy and internal forces.

Results and discussions
3.1. Structural properties Optimized lattice parameters of experimental and theoretical calculations and bond lengths of different ZnO models have been listed in table 1. Our first principle predictions of lattice parameters agree well with the experimental results, for example, for c lattice parameter it agrees within 0.04% [29], much better than the other theoretical calculations [30,31], indicating that our method of calculation is reliable. Lattice parameters a and c increase with increasing Na concentration, and this is due to the fact that the atomic radius of Na (0.95 Å) is larger than that of Zn (0.74 Å). However, Na i affects the parameters more and isotopically compared to Na Zn , as can be seen from the c/a values. These behaviors have been shown by others [18]. Position of Na in interstitial doped ZnO is at the octahedral site. However, four Na-O atomic distance along the positive c direction are different from the negative direction indicating shift from ideal octahedral position along the c axis which agrees with the shift determined by angular distribution of β − particles emission [20]. Co-doping of H with Na Zn increases c parameter more than both Na Zn and Na i models and the c/a is in between the two models. Change of the parameter c is important in case of hexagonal ZnO because it preferentially crystalizes along the c direction ([0001] direction) [33] thus controlling the in-plane stress of deposited thin films which can be modeled by the biaxial strain model [34] as, where c o is the corresponding value for bulk undoped ZnO (5.2047 Å) and c is the lattice parameter calculated from other models. The positive (negative) sign of s indicates that the films are in a state of tensile (compressive)

Formation energy
Assuming thermodynamic equilibrium and neglecting defect-defect interactions (i.e. in the dilute regime), the concentration of a native defect in a solid is determined by its formation energy E f through the relation [36], where N sites is the number of sites (including different configurations) per unit volume the defect can be incorporated on, k B is the Boltzmann constant and T the temperature in kelvin. Equation (2) shows that defects with high formation energies will occur in low concentrations. The energy appearing in equation (2) is, in principle, a free energy of formation. We have calculated the formation energy of the different neutral defect complex models using following equations: Where, E f (Na Zn ), E f (Na i ) and E f (Na Zn +H i ) indicate the formation energies of the substitutional and interstitial Na defects and substitutional Na and interstitial H co-doped defect complexes, respectively, E(ZnO+Na Zn ), E(ZnO+Na i ) and E(Na Zn +H i ) are the total energies of the supercell containing substitution, interstitial Na defects and co-doped Na Zn -H i complexes, E(ZnO undoped ) is the total energy of a perfect supercell of ZnO, E(Zn), E(Na) and E(H) are the total energies per atom of Zn, Na and H in their reference states, respectively. They have been calculated using pseudo atomic calculation scheme. Pseudo atomic calculation is a method of calculating total energy of isolated atom as their reference state using the energy of their outermost orbitals. In our case, 3d 10 -4s 2 , 2s 2 -2p 4 , 2s 2 -2p 6 -3s 1 and 1s 1 orbitals have been used to calculate pseudo atomic total energy for Zn, O, Na and H, respectively. Lower formation energy indicates a more stable system. All formation energy calculated from the above equations are used only for qualitative comparison. For more detailed quantitative approach readers can refer to these studies [14,19].
As shown in table 2 and figure 3(b), the formation energy of the interstitial model is negative and lower than the substitutional one which indicates that, the interstitial model is more stable than the substitutional one. The same behavior has been reported by other groups [19,37,38]. Higher concentration of Na Zn is unstable as has been shown by theoretical [18] and experimental results [35]. Presence of Na i defects increases stability of Na Zn defects thus, high solubility of Na atoms in ZnO lattice is only possible when both types of defect are present. H co-doping stabilizes both p-type (2Na Zn ) and n-type (Na Zn -Na i ) structures. Thus, H co-doping is a good approach for higher Na Zn concentration without increasing compensating Na i defects. Although, H also acts as a donor in ZnO, it is much easier to diffuse out than Na i defects [19,39].

Electronic structure
According to the optimized lattice structure, we have calculated the electronic band structures along with the high symmetry directions in the Brillouin zone for both undoped ZnO and doped ZnO models. Figure 4 shows the band structure and total density of states (TDOS) of pure ZnO. The band structures for different defect positions, different levels of H co-doping and type of doping conditions (i.e. p or n type) are shown in figure 5, 6 The Fermi level is set to zero. and 7, respectively. The partial density of states (PDOS) of undoped ZnO is shown in figure 8. According to figure 4, we can see that ZnO is a direct bandgap semiconductor with a bandgap of 3.394 eV at the G point, which is close to the experimental value 3.44 eV [40]. Our result is comparatively closer to the experimental measurement than those obtained by other methods such as LDA (0.76 eV) [41], GGA (0.731 eV) [42] and TB −mBJ (2.7 eV) [43]. The fermi level lies in the middle of the VBM (valence band maximum) and CBM (conduction band minimum) as expected from an intrinsic semiconductor. To clearly understand, which band corresponds to which orbital of which element, we have calculated the TDOS and PDOS for both undoped and doped models. From figures 4 and 8 we can see that, the top and bottom of the valence band and conduction band for undoped ZnO is mainly made out of O-2p and Zn-4s orbitals, respectively, which agrees with other results [30,44].
respectively. The Fermi level is set to zero.
In figure 5, we have calculated band structures of Na Zn , Na i and Na Zn -H i . All of them show direct band gap characteristic. For Na Zn the fermi level is near the VBM, indication of p-type conductivity. For both Na i and Na Zn −H i the fermi level is inside the conduction band, though shallower for Na Zn -H i . Thus, both show degenerate n-type conductivity and can act as compensating centers in p-type ZnO which have been shown by previous studies [12,14]. Band gap increases for Na Zn and decreases for both Na i and Na Zn -H i (table 2). For Na Zn two localize bands of Na-3s and Na-2p are present between-5 eV to 0 eV range thus contributing to upper valence bands and the conduction band is made of delocalized Na-2p orbital ( figure 8). In Na i model the conduction band has contribution from both Na-2s and Na-2p delocalized states at a deeper energy level of around 5 eV. However, the upper valence band has much less contribution from Na-3s and 2p than Na Zn . In case  of H co-doped model, the DOS is similar to Na Zn model with a shift towards lower energy. The upper conduction band has a strong contribution from H-1s orbital. In figure 6, effect of different levels of H compensation on band structure has been shown. All band structures indicate direct band gap while decreasing with compensation. Fermi level indicates that both 2Na Zn and 2Na Zn -H i models (corresponding to 0 and 50% compensation) are of p-type. However, for 50% compensation the fermi level is closer to the VBM (Δ 456 meV) indicating better p-type conductivity than no H co-doping. 100% compensation produces degenerate n-type conductivity. From figure 9, we can see that H-1s interacts with Na-2p and 3s channels and shifts the bands slightly towards higher energy which can explain the increased p-type behavior. In all the cases the valence band has a strong contribution from localized Na-3s and 2p orbitals and conduction band from delocalized Na-2p and H-1s (in case of compensation) orbitals. In figure 7, effect of p-type and n-type conditions on band structure is shown. Both Na i -Na Zn and Na i -Na Zn -H i show n-type degeneracy. However, in Na i -Na Zn    and conduction band respectively are influenced more by the Na zn −3s and 2p channels than Nai defect complexes. Delocalization of the 3 s channel in the Na i -Na Zn -H i model specially for Na Zn defect complex is due to interaction with H-1s channel. The interaction is more pronounced in the n-type model (Na i -Na Zn -H i ) than in the p-type model (2Na Zn -H i ). Contribution from H-1s orbital is much more pronounced in the n-type model than p-type in the lower conduction band region.

Optical properties
In the linear response range, the solid macroscopic optical response function can usually be described by the frequency dependent dielectric function In figure 12, absorption spectrum for undoped and doped models have been calculated for a certain range of photon energy. The absorption edge in the higher energy area (3-4 eV/300-400 nm) in undoped ZnO spectra can be attributed to the transition between VB and CB, corresponds to absorption in the ultra-violate range. The absorption edge corresponds to the band gap, and the absorption spectrum is induced by the inter-band transition of electrons. The absorption edge moves to the higher energy area (3.5-4.5 eV) for Na Zn , indicating that the band gap increase after doping, shows a blue shift compared to un-doped ZnO. Na Zn shows a second absorption edge at a very low energy level range (92-790 meV) corresponding to electron transition between VB and impurity states above it. Absorption spectrum of Na Zn -H i is similar to Na Zn except that the second absorption edge is missing, indicating the impurity states have been filled by H-1s electrons. All our models show very little absorption in the visible range (see insets in the 1.8-3.1 eV range) which is appropriate for transparent conductive oxide (TCO) applications. The calculated refractive index ( ) w n and extinction coefficient ( ) w k of undoped, Na Zn , Na i single-doped, and Na Zn -H i co-doped ZnO are shown in figure 13. The static refractive indices ( ) n 0 of undoped, Na Zn , Na i single-doped, and Na Zn -H i co-doped ZnO are about 1.58, 3.92, 1.573 and 1.573, respectively. In the energy range of 0-16 eV, the ( ) w k of undoped, Na i single-doped, and Na Zn -H i co-doped ZnO increases. Though, Na Zn shows similar behavior, at 500 meV the value of k is non-zero (1.96), similar to the e i behavior.
The effect of H co-doping on absorption and dielectric function have been shown in figure 14. Other models show similar behavior to that of Na Zn , the only difference is that 100% compensation (2Na Zn -H i ) shows a red shift for absorption spectrum. Second absorption edge and non-zero value at near zero energy level for 0 and 50% compensation is the evidence of empty impurity states. 100% compensation fills the empty states. Decreased absorption near the second absorption edge of 2Na Zn -H i model compared to 2Na Zn , is evidence of partially filled impurity states because of H co-doping. In both compensated and non-compensated cases, for n-type condition a blue shift near the first absorption edge can be seen in their absorption spectrum ( figure 15) with respect to p-type condition. Na Zn -Na i does not has a second absorption edge like 2Na Zn which is evidence of Na i -3s electrons filling empty states created by Na Zn in O-2p orbitals.

Conclusion
In this study, the structural, electronic and optical properties of undoped, Na single-doped, and Na-H co-doped ZnO were evaluated by first-principles calculations based on DFT. Na i sits in the slightly shifted octahedral position and increases both lattice parameters (a and c) compared to Na Zn for the same doping concentration. Co-doping of H increases lattice parameter c much more than both substitutional and interstitial doping. However, compensation can help stabilize the structure by reducing strain as it greatly reduces Na-O bond distance parallel to c axis, however, it has to be done at a high concentration level. Strong O-H bond is formed in all our H co-doped models and they are not affected by both compensation level and doping condition. From all our models, H co-doped ZnO are the most stable and substitutional Na doped are the least stable. Additionally, the stability increases with H co-doping and decreases with Na doping in ZnO respectively for higher concentrations. Thus, p-type conductivity cannot be achieved in ZnO by mono-doping only, rather co-doping with donor atoms is necessary. From our calculations, H is a better candidate as it retains the host structure, forms stable complexes, and can be diffused out at a low temperature. From the electronic structure calculations, we have seen that 50% H compensation produces the best p-type conductivity as H-1s interacts with Na-2p and 3 s channels and shifts the bands slightly towards higher energy. Thus, during annealing process retention of some percentage of H can be beneficial for p-type conductivity. In summary, annealing temperature is not only the important factor here, annealing duration must be controlled too to get good p-type behavior. From the optical calculations, we have found that all our models are transparent in the visible range. A shallow acceptor level can be derived from the fact that a new and strong peak at about 174 meV is present for Na Zn which is due to transition of electrons from Na-3s to the acceptor levels created by the unoccupied p states of O. This acceptor level can be filled at different levels depending on the H co-doping level. From the results above, Na-H codoping is a good strategy for achieving good p-type conductivity with the optimized level of H compensation by controlling annealing temperature and time.