Diffusion of indium in single crystal zinc oxide: a comparison between group III donors

Doping of zinc oxide (ZnO) with donor atoms, such as the group III elements aluminum (Al) and gallium (Ga) is common practice in realizing transparent conductive oxides (TCO) [1–4]. The main driver for using ZnO-based TCO's has been, and still is, the higher cost of the more commonly used indium-tin-oxide (ITO). ZnO is more abundant, and with reported resistivity as low as 8–9 × 10⁻⁵ Ω cm for Al- and Ga -doped ZnO [5, 6], i.e. comparable to that of ITO (7.2 × 10⁻⁵ Ω cm [4, 7]), makes ZnO desirable for use as a transparent electrode. However, for thin films (<100–200 nm) ZnO-based TCO's are still inferior to ITO in terms of resistivity stability at elevated temperatures in air ambience [3, 8]. In addition to Al and Ga, In has also been demonstrated to produce highly conductive ZnO layers, where a resistivity of about 8–9 × 10⁻⁴ Ω cm was reported for In-doped thin films prepared by spray pyrolysis [9] and RF magnetron sputtering [1]. This should be compared with a resistivity of 1.9 × 10⁻⁴ Ω cm and 5.1 × 10⁻⁴ Ω cm for Al- and Ga -doped ZnO, respectively, using similar magnetron sputtering deposition techniques [1]. Furthermore, it has previously been shown that self-compensation occurs in Al and Ga doped ZnO [10, 11], resulting in a conductivity limit and reduced applicability as TCO. The self-compensation was explained by the formation of doubly negatively charged zinc vacancies (${{V}}_{\mathrm{Zn}}^{2-}$) and a complex between ${{V}}_{\mathrm{Zn}}^{2-}$ and substitutional Al at zinc site (${\mathrm{Al}}_{\mathrm{Zn}}^{+}$) for Al-doped ZnO [11, 12]. Thus, given the similar valence character of In and its observed donor behaviour in ZnO, ${{V}}_{\mathrm{Zn}}^{2-}$ will presumably play an important role also in the case of In-doped ZnO-based TCO's.

Studies of donor dopant diffusion in ZnO enable information about the donor-defect interactions involved. Accordingly, such studies are a viable path to elucidate the mechanisms limiting the net carrier concentration. Previous experimental investigations of In diffusion in ZnO report activation energies ranging from 1.17 to 3.16 eV, with diffusion constants spanning over 9 orders of magnitude [13–15]. Clearly, this large inconsistency in literature data, including the dissimilar diffusion energies reported for two identical experiments, [13] and [14], urges the need for further studies of the In diffusion in ZnO.

In more recent years, calculations using density functional theory (DFT) have emerged as a method to estimate defect formation energies and also favourable migration paths for self- and impurity diffusion. Hence, by presuming a prevailing diffusion mechanism, DFT-calculations enable predictions of dopant diffusion activation energies. In the case of In diffusion in ZnO, a migration barrier for the prevalent (${\mathrm{In}}_{\mathrm{Zn}}{{V}}_{\mathrm{Zn}}{)}^{-}$ pair has previously been predicted to be 1.22 eV (along the c-orientation), with an estimated overall diffusion activation energy of 1.77 eV under O-rich and n-type conditions [16].

In this work, we use secondary ion mass spectrometry (SIMS) to investigate the diffusion of In in single crystalline bulk ZnO, utilizing an effectively inexhaustible surface layer source. Similar to our previous reports on Al and Ga diffusion in ZnO [11, 17, 18], the experimental results are analyzed using a reaction–diffusion (RD) model that combines with DFT results from [16]. From this, we show that the formation- and migration energetics reported in [16] are self-consistent with our exeprimental results on In diffusion.

2.1. Experimental

2.2. Computational methods

2.2.1. Reaction diffusion model

To describe the experimental In diffusion depth profiles, a RD type model is employed that is based on a complete set of reaction diffusion differential equations (see e.g. [11, 17–22] for a similar and general treatment), and considers the formation, migration and dissociation of the (${\mathrm{In}}_{\mathrm{Zn}}{{V}}_{\mathrm{Zn}}{)}^{-}$ pair through the following reaction

$$\begin{eqnarray}&&{{V}}_{\mathrm{Zn}}^{2-}+{{\rm{In}}}_{\mathrm{Zn}}^{+}\rightleftharpoons {\left({{\rm{In}}}_{\mathrm{Zn}}{{V}}_{\mathrm{Zn}}\right)}^{-}.\end{eqnarray} \tag{ 1 }$$

Here, ${{V}}_{\mathrm{Zn}}^{2-}$ is highly mobile while ${\mathrm{In}}_{\mathrm{Zn}}^{+}$ is regarded as immobile at the studied temperatures. We can write the following system of RD equations for the migrating pairs

$$\begin{eqnarray}\begin{array}{rcl}\displaystyle \frac{\partial {C}_{s}}{\partial t} & = & {k}_{s}{C}_{p}-{k}_{p}{C}_{v}{C}_{s}\\ \displaystyle \frac{\partial {C}_{p}}{\partial t} & = & {D}_{p}\displaystyle \frac{{\partial }^{2}{C}_{p}}{\partial {x}^{2}}-\displaystyle \frac{\partial {C}_{s}}{\partial t},\end{array}\end{eqnarray} \tag{ 2 }$$

with C_s, C_p and C_v representing the concentration of substitutional In${}_{\mathrm{Zn}}^{+},{({{\rm{In}}}_{\mathrm{Zn}}{{V}}_{\mathrm{Zn}})}^{-}$ pairs and ${{V}}_{\mathrm{Zn}}^{2-}$, respectively. k_s is the dissociation rate of the pairs whilst ${k}_{p}{C}_{v}$ is the corresponding formation rate. It was shown previously for Ga diffusion in ZnO [17] that when the transport capacity of ${{V}}_{\mathrm{Zn}}^{2-}$ is much higher than that of the pair (i.e. ${C}_{v}{D}_{v}\gg {C}_{p}{D}_{p}$), k_s and k_p are reduced to depend only on the binding energy of the pair and the effective radius (1 nm) for capturing ${{V}}_{\mathrm{Zn}}^{2-}$ by the substitutional dopant, respectively. In this regard, a pair binding energy of ${E}_{b}{({{\rm{In}}}_{\mathrm{Zn}}{{V}}_{\mathrm{Zn}})}^{-}=1.19\,\mathrm{eV}$ has been predicted by Steiauf et al [16], and is used as a fixed parameter in our simulations.

In order to solve the above RD equations (equations (2)), C_v(x, t) needs to be known. Profiles of C_v(x, t) can be estimated from DFT predictions of the ${{V}}_{\mathrm{Zn}}^{2-}$ formation energy (${E}_{f}({V}_{\mathrm{Zn}}^{2-})$) and also accounting for the effect of band-gap narrowing at high temperatures. That is, the distribution of ${{V}}_{\mathrm{Zn}}^{2-}$ can be expressed as [11, 18]

$$\begin{eqnarray}&&{C}_{v}(x,t)={N}_{s}{e}^{-({E}_{f}({V}_{\mathrm{Zn}}^{2-})/{k}_{B}T)}{\left(\displaystyle \frac{n(x,t)}{{N}_{c}(T)}\right)}^{2},\end{eqnarray} \tag{ 3 }$$

where N_s is the number of substitutional zinc lattice sites, n equals $({C}_{s}-{C}_{p}-2{C}_{v})$ and accounts for the charge neutrality of the system, and N_c is the effective density of states in the conduction band. The high transport capacity of ${{V}}_{\mathrm{Zn}}^{2-}$ implies an almost instantaneous establishment of an equilibrium C_v that is governed by the local Fermi-level position. The vacancy formation energy can in turn be expressed by ${E}_{f}({V}_{\mathrm{Zn}}^{2-})\,={E}_{f,0}({V}_{\mathrm{Zn}}^{2-})-2{\epsilon }_{{\rm{F}}}$, where ${E}_{f,0}({V}_{\mathrm{Zn}}^{2-})$ is the formation energy at the valence band edge, set to 7.0 eV in our simulations as guided by previous DFT reports [16, 23–25], and _F is the Fermi-level position referenced to the valence band edge. At the simulated temperatures (900 °C–1150 °C), the band gap of ZnO (3.3 eV at 295 K) varies between 2.77 and 2.64 eV, assuming a similar temperature dependence for the band gap narrowing (−0.52 meV °C^–1) as that reported up to 500 °C [26]. We further consider a fixed valence band edge, that is, the absolute value for the conduction band edge equates the band gap at all temperatures. For a more detailed discussion of the RD model used in this work, see [11, 17, 18].

2.2.2. Density functional calculations

First-principles calculations were performed by using the Heyd–Scuseria–Ernzerhof [27] hybrid functional and the projector augmented wave method [28–30], as implemented in the VASP code [31, 32]. The fraction of screened Hartree–Fock exchange was set to α = 37.5% [33], which yields a band gap (3.42 eV) and lattice parameters (a = 3.244 Å  and c = 5.194 Å) that are in excellent agreement with experimental values [34, 35]. All defect calculations were performed using a plane-wave energy cutoff of 500 eV, a special k-point at $k=\left(\tfrac{1}{4},\tfrac{1}{4},\tfrac{1}{4}\right)$, and a 96 atom sized wurtzite supercell [25]. Defect formation energies were calculated by following the well established formalism outlined in [36, 37]. For instance, the formation energy of In_Zn in charge-state q is given by

$$\begin{eqnarray}&&{E}_{{f}}({\mathrm{In}}_{\mathrm{Zn}}^{q})={E}_{\mathrm{tot}}({\mathrm{In}}_{\mathrm{Zn}}^{q})-{E}_{\mathrm{tot}}^{\mathrm{bulk}}+{\mu }_{\mathrm{Zn}}-{\mu }_{\mathrm{In}}+q{\epsilon }_{{\rm{F}}},\end{eqnarray} \tag{ 4 }$$

where ${E}_{\mathrm{tot}}({\mathrm{In}}_{\mathrm{Zn}}^{q})$ and ${E}_{\mathrm{tot}}^{\mathrm{bulk}}$ denote the total energy of the defect-containing and pristine supercells, and μ_Zn and μ_In are the chemical potential of the removed Zn- and added In-atom, respectively. For charged defects, we applied the anisotropic [38] Freysoldt–Neugebauer–Van de Walle finite-size correction [39, 40]. Oxygen rich conditions are considered, where μ_Zn corresponds to the total energy per bulk metallic Zn atom plus the formation enthalpy of ZnO, i.e. ${\mu }_{\mathrm{Zn}}={E}_{\mathrm{tot}}(\mathrm{Zn})+{\rm{\Delta }}{H}_{{f}}(\mathrm{ZnO})$. The solubility of In is limited by the formation of In₂O₃, and under oxygen rich conditions ${\mu }_{\mathrm{In}}={E}_{\mathrm{tot}}(\mathrm{In})+\tfrac{1}{2}{\rm{\Delta }}{H}_{{f}}({\mathrm{In}}_{2}{{\rm{O}}}_{3})$.

Figure 1 shows the In concentration versus depth profiles for the isochronally heat treated sample, as obtained from SIMS measurements. A migration of In from the In-doped ZnO film into the ZnO bulk is observable above 800 °C, while heat treatments conducted at lower temperatures (600 °C–800 °C) revealed no change in the In distribution. For treatments below 1150 °C, the In concentration in the 0.3 μm thick deposited film is over 2 orders of magnitude higher than the resulting apparent solid solubility of In in ZnO (i.e. the concentration of In at the bulk surface). At 1150 °C, In start to deplete from the film, but still holds a concentration that is one order of magnitude higher than in the bulk. The box shape of the In diffusion profiles are of similar character as that previously reported for Al [11] and Ga [17] in ZnO. However, as shown in figure 2, the apparent solid solubility of In is about one order of magnitude lower than that observed for Al and Ga at similar temperatures [11, 17]. In addition, the diffusion length of In is significantly longer than that of Al and Ga at comparable temperatures (see inset in figure 2), with the deep end of the In profile extending ∼11 μm after the 1050 °C treatment, as compared to ∼6 μm and ∼0.5 μm for Ga [17] and Al [11], respectively.

Zoom In Zoom Out Reset image size

Zoom In Zoom Out Reset image size

Figure 3 shows the calculated formation energies of ${{V}}_{\mathrm{Zn}}$, and substitutional ${\mathrm{In}}_{\mathrm{Zn}}$, ${\mathrm{Al}}_{\mathrm{Zn}}$ and ${\mathrm{Ga}}_{\mathrm{Zn}}$, versus the Fermi-level position. As seen, the formation energy of ${\mathrm{In}}_{\mathrm{Zn}}$ is higher than those of ${\mathrm{Al}}_{\mathrm{Zn}}$ and ${\mathrm{Ga}}_{\mathrm{Zn}}$. The difference in the apparent solid solubilities between Al and Ga in figure 2 is in close agreement with that expected from the DFT results in figure 3; the calculated energy required to form ${\mathrm{Al}}_{\mathrm{Zn}}$ is 0.23 eV higher than that of ${\mathrm{Ga}}_{\mathrm{Zn}}$, which is consistent with the experimental values of (1.06 ± 0.07) eV and (0.70 ± 0.13) eV, respectively. Moreover, the pre-exponential factors for Al and Ga are close to the number of lattice sites in ZnO ($4\times {10}^{22}\ {\mathrm{cm}}^{-3}$ Zn-sites), i.e. the extracted apparent solid solubility corresponds to the actual solubility. The absolute DFT values for E_f of ${\mathrm{Al}}_{\mathrm{Zn}}$ and ${\mathrm{Ga}}_{\mathrm{Zn}}$ in the studied temperature range (see figure 3) are also close to the experimental ones in figure 2, corroborating that the experimental ambient conditions (air) are close to that of the O-rich conditions in DFT. Further, from the DFT-data in figure 3, E_f of ${\mathrm{In}}_{\mathrm{Zn}}$ is expected to be 0.5 eV higher than that of ${\mathrm{Al}}_{\mathrm{Zn}}$. However, the apparent solid solubility of In (figure 2) does not corroborate this trend, but rather exhibit a lower value (by 0.6 eV) as compared to Al. Here, it should be noted that the pre-exponential factor obtained for In (figure 2) is significantly lower than that found for Al and Ga, and may indicate that the incorporation of In is limited by a different process. For instance, one possibility is that the rate of transport of In from the polycrystalline surface film and into the bulk crystal is the limiting process. Indeed, previous diffusion experiments [13, 14] using In ions implanted into ZnO samples reveal an In diffusion plateau with a concentration of $\sim 1\times {10}^{19}\ {\mathrm{cm}}^{-3}$ already after annealing at 850 °C in air, which is notably higher than observed in our indiffused sample. On the other hand, even though the incorporation of In is transport limited, steady state surface conditions still apply, and the assumption of an inexhaustible diffusion source remains valid. ${\epsilon }_{{\rm{F}}}$ may, however, be pinned slightly below the conduction band edge causing an increase of E_f(${{V}}_{\mathrm{Zn}}^{2-}$) and thus also of the overall apparent diffusion activation energy.

Zoom In Zoom Out Reset image size

To investigate the apparent transport limit in more detail, high-angle annular dark-field (ADF) imaging (HAADF) STEM analysis was conducted on the deposited In-doped film (after the 1150 °C treatment) and, for comparison, ADF imaging was conducted on the deposited Ga-doped film studied in [17] (after the 1050 °C treatment). Figure 4 shows the results obtained from the STEM analysis and reveal vertical lines of increased brightness for both the In- and Ga doped films. These lines may indicate the presence of inversion domain boundaries (see [41–44]) decorated with In and Ga, respectively. The observed contrast comes from the Z contrast that may be observed in ADF and HAADF. The appearance of the vertical lines are more frequent in the In-doped film, which could point to the fact that a larger fraction of In is trapped at inversion domain boundaries. If this is indeed the case, it may partly explain the limited flux of In from the deposited film as compared to the Ga-doped film. It can be noted that inversion domain boundaries have previously been observed for several elements in ZnO, including In [41, 42], Fe [41, 42], and Li [44] containing ZnO, and suggest that this is a rather common phenomenon in impurity doped ZnO. However, to fully understand this effect a dedicated study of the thin film-properties is needed, which is beyond the scope of this paper. It should also be noted that even though the lines in figure 4 are more frequent in the In-doped film, each line in the Ga-doped film seem to consist of several (3–5) atomic layers of Ga as compared to the single atomic layers observed for In.

Zoom In Zoom Out Reset image size

The experimental diffusion profiles have been modelled using the RD model introduced in section 2.2.1, with the resulting best fits shown as solid lines in figure 1. The RD model with the migration barrier of (${\mathrm{In}}_{\mathrm{Zn}}{{V}}_{\mathrm{Zn}}{)}^{-}$ as the controlling/limiting process for diffusion, closely reproduce the experimental profiles. In the simulations, only C_p and D_p are used as fitting parameters, where the former is treated as an influx boundary of already formed (${\mathrm{In}}_{\mathrm{Zn}}{{V}}_{\mathrm{Zn}}{)}^{-}$ pairs (formed in the film and 'injected' into the bulk crystal). Under equilibrium conditions and given that In mainly dissolve at substitutional zinc site, and that C_s ≫ C_p, the effective diffusivity can be expressed as (see [20, 21, 45]):

$$\begin{eqnarray}&&{D}_{{\rm{eff}}}=\displaystyle \frac{{C}_{p}^{{\rm{eq}}}{D}_{p}}{{C}_{s}^{{\rm{eq}}}},\end{eqnarray} \tag{ 5 }$$

where ${C}_{s}^{{\rm{eq}}}$ is taken as the apparent solubility ${C}_{{\rm{tot}}}^{{\rm{In}}}$ and ${C}_{p}^{{\rm{eq}}}$ is the concentration of pairs at the bulk surface. Figure 5 shows the temperature dependence of D_eff, revealing two different activated processes. At temperatures above 950 °C an apparent activation energy of ${E}_{a}({\rm{In}})=2.2\,\mathrm{eV}$ with a diffusion constant of D₀ = 4 × 10⁻² cm² s⁻¹ follows for the In diffusion. However, at lower temperatures the extracted diffusivities reveal a very different activated process. Since the diffusion prefactor is merely a structural entity [46], ${D}_{0}\approx {10}^{-2}\ {{\rm{cm}}}^{2}\,{{\rm{s}}}^{-1}$ for ZnO (see e.g. [17]), the large prefactor of 2 × 10⁶ cm² s⁻¹ obtained at low temperatures suggest that it does not represent a different diffusion mechanism, but may instead imply a transient diffusion behavior at these initial annealing stages. Indeed, as can be seen from the measured In concentration in figure 2, at the initial temperatures ${C}_{{\rm{tot}}}^{{\rm{In}}}$ deviates slightly from the trend at higher temperatures, indicative of a transient regime. Therefore, only diffusion profiles above 950 °C is considered representative in the extraction of the In diffusion activation energy. Table 1 lists the present results together with those reported previously in the literature on the diffusion of In in ZnO.

Zoom In Zoom Out Reset image size

Table 1.  Extracted diffusion parameter values for In in ZnO, together with previous theoretical and experimental results from the literature. With D₀ being the pre-exponential factor, E_m is the migration barrier for the exchange process of ${\mathrm{In}}_{\mathrm{Zn}}{{V}}_{\mathrm{Zn}}$ and E_a is the apparent activation energy for diffusion.

	Type	D0	Ea(In)	Em(${\mathrm{In}}_{\mathrm{Zn}}{{\rm{V}}}_{\mathrm{Zn}}$)−
		(cm2 s–1)	(eV)	(eV)
Present work	Exp	4 × 10−2	2.2 ± 0.2	—
Thomas [15]	Exp	2 × 102	3.16	—
Sakaguchi et al [14]	Exp	1.1	2.68	—
Nakagawa et al [13]	Exp	2.9 × 10−7	1.17	—
Steiauf et al [16]	Theo	—	1.77	1.22

In our simulations, only the product ${C}_{p}^{{\rm{eq}}}{D}_{p}$ is a unique quantity, and from experiments alone it is only possible to separate the two quantities in cases where (i) the pair is not the controlling/limiting process for the diffusion, or (ii) by isoconcentration experiments. Assuming that equation (1) is the prevailing diffusion mechanism, both (i) and (ii) would result in a shape of the diffusion profiles being very different from the present ones in figure 1, resembling that of erfc(x) functions [18, 20, 47]. However, as we have recently demonstrated for Al [18] and Ga [17] diffusion, it may be possible to separate C_p and D_p also at conditions different from (i) and (ii) by introducing DFT results to the RD-simulations (see section 2.2.1). In contrast to Al and Ga, applying the same procedure (as outlined in [17]) to determine D_p fails in the present case of In. This reveals an inadequacy of our RD model at these conditions, and may be due to the much lower dopant concentration in the present experiment. This may indicate that our RD simulations do not sufficiently account for the concentration and distribution of ${{V}}_{\mathrm{Zn}}^{2-}$ (equation (3)) in this somewhat moderate n-type conductive sample, and further experiments may be necessary to resolve it. Nonetheless, the extraction of ${D}_{{\rm{eff}}}$ (equation (5)) will not be affected and the determination of E_a(In) = 2.2 eV is still valid. It should be noted that since we could not determine D_p, similar results (${D}_{{\rm{eff}}}$) can also be obtained using a much simpler approach, e.g. Fair's model (see [17]), instead of using a system of RD equations.

On the other hand, the overall diffusion processes can be approximated through [17]

$$\begin{eqnarray}&&{E}_{a}{({\rm{In}})={E}_{{f}}{({V}_{\mathrm{Zn}}^{2-})-{E}_{b}({{\rm{In}}}_{\mathrm{Zn}}{V}_{\mathrm{Zn}})}^{-}+{E}_{m}({{\rm{In}}}_{\mathrm{Zn}}{V}_{\mathrm{Zn}})}^{-}+{k}_{B}T,\end{eqnarray} \tag{ 6 }$$

where the term k_BT is about 0.1 eV at the studied temperatures. _F is positioned about 0.2 eV below the conduction band edge, which combined with the DFT-results in figure 3 yield ${E}_{f}({V}_{\mathrm{Zn}}^{2-})$ values between 2.0–2.1 eV at the plateau of the In profiles during diffusion. Thus, using our extracted value of E_a = 2.2 eV (together with E_b = 1.19 eV) this gives ${E}_{m}{({{\rm{In}}}_{\mathrm{Zn}}{V}_{\mathrm{Zn}})}^{-}=1.2\pm 0.3\,\mathrm{eV}$, which is in excellent agreement with that predicted by DFT calculations (1.22 eV) [16], demonstrating self-consistency in the theoretical predictions by [16].

Interestingly, our In results differ from those reported previously by Thomas [15] (E_a(In) = 3.16 eV), Sakaguchi et al [14] (E_a(In) = 2.68 eV) and Nakagawa et al [13] (E_a(In) = 1.17 eV). It is, however, not obvious why the activation energy obtained by Sakaguchi et al [14] (2.68 eV) differs so much (1.5 eV) from that found by Nakagawa et al [13] (1.17 eV) performing nominally the identical experiment. It can also be noted that Nakagawa et al [13] report a surprisingly low pre-exponential factor, D₀, of 2.9 × 10⁻⁷ cm² s⁻¹.

To summarize, we have studied the diffusion of In in single crystalline ZnO and found that it can be well described as mediated by ${{V}}_{\mathrm{Zn}}^{2-}$, through the intermediate formation and dissociation of (${\mathrm{In}}_{\mathrm{Zn}}{{V}}_{\mathrm{Zn}}$)⁻ pairs. An activation energy of 2.2 ± 0.2 eV, with a diffusion constant of 4 × 10⁻² cm² s⁻¹ is obtained for the diffusion of In at temperatures 1000 °C up to 1150 °C. This is lower than that reported from previous experiments in the literature [14, 15], but is in good agreement with recent DFT results [16] when considering the somewhat moderate n-type conductivity in the present sample with a local Fermi-level position at ∼0.2 eV below the conduction band edge. Moreover, the theoretical predictions of E_f(${{V}}_{\mathrm{Zn}}^{2-}$), ${E}_{b}{({{\rm{In}}}_{\mathrm{Zn}}{{V}}_{\mathrm{Zn}})}^{-}$ and ${E}_{m}{({{\rm{In}}}_{\mathrm{Zn}}{{V}}_{\mathrm{Zn}})}^{-}$ made in [16] are found to be self-consistent with our experimental In diffusion profiles.

We gratefully acknowledge the financial support from the Research Council of Norway for funding of the DYNAZOx-project (no. 221992), Salient (no. 239895) and NORTEM (197405/F50), and the University of Oslo for funding through FUNDAMENT (no.151131, Fripro Toppforsk program), and the Norwegian Micro- and Nano-Fabrication Facility (NorFab 245963).