Computational study and ion diffusion analyses of native defects and indium alloying in β-Ga₂O₃ structures

Wide band gap oxides – with a gap greater than 3 eV – are part of the next generation of materials to revolutionize the semiconductor industry. These materials have multiple tunable device parameters,¹ which make them suitable for high-power, high-chemically stable environments.^2,3 Specifically, the monoclinic beta phase of gallium oxide, β-Ga₂O₃, with a band gap of 4.8 eV and the ability to be n-type doped, along with a relatively high electron mobility of 200 cm² V⁻¹ s¹, has shown immense promise over the last decade of research. Also, numerous dopants have been studied that show the tunability of β-Ga₂O₃ through defect engineering.^4,5 Recently, we have shown that indium is a stable isoelectronic dopant for gallium in β-Ga₂O₃,⁶ where substitutional and interstitial indium defect complex formation results in impurity band conduction and an increase in n-type conductivity. The inclusion of small concentrations of indium, used to help lattice match NiO with Ga₂O₃, resulted in an increase in n-type conduction, which was validated through first principles calculations. This highlights the versatility of β-Ga₂O₃ as a tunable wide band gap semiconductor through the understanding of defects in the bulk material.

Defects are ubiquitous in devices and must be included when studying material properties from first principles. Both intrinsic (vacancy, host-interstitial, and/or antisite) and extrinsic (substitutional, foreign-interstitial) defects may be studied in the dilute limit using the supercell approach.⁷ Defects in Ga₂O₃ have been studied for common point defects like Si using hybrid density functional theory (DFT) which highlights the advancement of ab-initio simulations.⁸ Our previous study also using hybrid DFT⁶ highlighted the stability of a near valence band In_Ga defect state, and a near conduction band In_Ga + In_i defect complex, with a transition occurring within 0.5 eV of the conduction band minimum; here we use the standard Kroger-Vink notation.⁹ While static properties may be used to study likely formation of defects, dynamic properties require more elaborate thermodynamic methodologies.

To understand the nature of transient properties like defect diffusion in β-Ga₂O₃ using computational methods, one may turn to molecular dynamics (MD) simulations.^10,11 One evolves Newton's equations of motion, solving for atomic trajectories using periodic boundary conditions along with appropriate statistical ensembles. The general schema is to (1) evolve the system from t to $t+δt$ and update the positions of each particle at the evolved time, (2) compute the atomic velocities at the midpoint of the time step, i.e., $t+ δ t 2$⁠, (3) compute the forces and accelerations of each particle at the evolved time, and (4) evolve the velocity to the full time step. This process is iterative and results in the time evolution of each defective system which allows one to study defect trajectories.

In this work, we use a combination of DFT and MD simulations to study the diffusion of defects in β-Ga₂O₃. The paper is outlined as follows: Sec. II outlines the DFT and MD parameters and potentials used, Sec. III details the results for structural, ionic, and electronic simulations, and Sec. IV summarizes our findings.

The results presented in this study were obtained from a theoretical investigation by integrating two methodologies: density functional theory (DFT) calculations and molecular dynamics (MD) simulations. The ab initio calculations, based on DFT,^12,13 have been performed using the Vienna Ab Initio Simulation Package (VASP)^14,15 and the projector augmented wave (PAW) pseudopotentials as implemented in VASP. The MD simulations were carried out using the open-source software Large-scale Atomic/Molecular Massively Parallel Simulator (LAMMPS).^16,17

Gallium oxide was studied in the monoclinic (β phase) crystal structure in eight different configurations including both intrinsic and indium-based defects: pure, interstitial oxygen (O_i), oxygen vacancy (V_O), interstitial gallium (Ga_i), gallium vacancy (V_Ga), interstitial indium (In_i), substitutional indium for gallium (In_Ga) and the impurity complex (In_GaIn_i). The initial positions of all considered configurations are displayed in Fig. 1. For DFT simulations, supercells of 120 atoms (1 × 3 × 2) were used, considering 48 gallium atoms and 72 oxygen atoms. On the other hand, for the MD calculations, the supercell has 960 atoms (2 × 6 × 4), composed of 384 gallium atoms and 576 oxygen atoms. The number of defects for each simulation is one for DFT calculations and eight for MD calculations; in this way, we guarantee the same concentration in both methods, which is 8.33% per cell. The lattice angles are kept fixed at α = 90°, β = 103.78° and γ = 90°. It was observed that the presence of an In_i in the largest Ga-O cage coupled with an In_Ga occupying an octahedral cation site are more energetically favorable.

The interatomic interaction is described in pairs. The self-interaction between O²⁻ ions and the interaction between O²⁻ and Ga³⁺ have been used to structurally describe Ga₂O₃.²³ The In³⁺ – O²⁻ interatomic potential parameters were derived separately. To obtain these interatomic potential parameters, we have used the general utility lattice program (GULP) code^24–27 to fit them according to the structural data for orthorhombic In₂O₃. The full set of interatomic potentials used is presented in Table I. The total charge state of our MD simulations is always neutral. Therefore, in the presence of interstitial or vacancy defects the system charge was neutralized by compensating the charges of the ion removed or added to the gallium and indium atoms of the structure. This method allows to simulate only low defect concentrations.²⁸ 

In the DFT calculations we used the generalized gradient approximation (GGA) within the Perdew-Burke-Erzenhof (PBE) framework³⁰ for the electronic exchange-correlation potential. The following valence electronic distributions were used for each PAW pseudopotential: Ga – 3d¹⁰ 4s² 4p¹, In – 4d¹⁰ 5s² 5p¹, and O – 2s² 2p⁴. A value of 500 eV has been set for the cutoff energy in the plane wave expansions. The electronic and ionic relaxations have been performed until the residual forces on the ions were found to be <10 meV/Å. A $1×1×1$ Monkhorst–Pack³¹ k-points mesh has been adopted for Brillouin zone (BZ) integrations. Details of the convergence criteria for the kpoints-mesh and energy cutoff are presented in Support Information, as shown in Figs. S1(a) and S1(b),⁴³ respectively.

Chemical bonding within the structures were examined through an analysis of the charge density, the electron localization function (ELF), and the Bader charge. The Bader charge was computed using the Henkelman's code,^32,33 which was integrated into VASP.

The structural properties of both pure β-Ga₂O₃ and those containing defects outlined in Fig. 1 were initially examined through their cell optimization. Figure 2(a) shows the total energy versus volume for all structures obtained by DFT and MD calculations. Each structure has been isotropically compressed and expanded and its atomic positions have been fully relaxed to the minimum energy. For the sake of comparison between methods, the energy values were adjusted to the unit cell, which contains 4 units of the β-Ga₂O₃ molecule. The minimum energy $( E 0 ∗)$ found for each method was set to zero. The optimized volume, corresponding to the minimal energy, as well as the related lattice constants a, b and c, are displayed in Table II. It is observed that both methods present very similar results, and in a good agreement with experimental data, where a, b, c parameters and β angle are equal to 12.21 Å, 3.04 Å, 5.80 Å and 103.8°, respectively.³⁴ Nevertheless, when compared to DFT, MD yields slightly smaller optimized volume values, in the case of the pure structure, there was a 1.5% difference in volume and 0.5% variation in lattice parameters. The DFT relaxed structures are shown in Fig. S1.⁴³ 

The bond angles and bond lengths of the β-Ga₂O₃ – pure structure revealed good agreement between the methods employed and the experimental values. The MD results have been obtained after system thermalization, in three steps employing a Nosé-Hoover thermostat and barostat for a total time of 1500 ps, at a temperature of 10 K (for more details see the description in Sec. II). The gallium atoms assume two different configurations in the unit cell, corresponding to tetrahedral and octahedral sites (see the illustration in Fig. S2).^43, Table III shows the calculated bond lengths Ga-O values, gathered from the RDF function in MD simulations, for the two different sites in good agreement with those reported in the literature. The O-Ga-O and Ga-O-Ga bond angles were analyzed by the angular distribution functions (ADF) displayed in Fig. S3.⁴³ At very low temperatures, one can differentiate the peaks associated to tetrahedral angles at interval $100 °− 120 °$ and to octahedral angles at intervals $∼ 80 °− 100 °$ and $∼ 160 °− 180 °$ [see Fig. S3(a)].⁴³ These results are coherent with the ADF obtained by DFT calculations, highlighting that both methods well predict the octahedral and tetrahedral bond angles.

To verify the thermal stability of both pure and defective structures, we have performed a series of MD simulations by varying the temperature from 200 to 2000 K at ambient pressure. Figures 3(a), 3(c) and 3(e) depict the ADF for the Ga-O-Ga and O-Ga-O angles, and the RDF for the Ga-O pair, respectively, for the pure structure at different temperatures. Both distributions were obtained after thermalization and averaged over the 5000 ps of simulation. Notably, the crystalline arrangement is not significantly affected by temperature variations, although there is an increase in dispersion, which is expected given the greater kinetic energy of the atoms. Figures 3(b), 3(d) and 3(f) present a similar analysis, but now considering both pure and defective structures at 1000 K. Neither bond angles nor Ga-O pair distance show significant alteration by the introduction of defects, indicating that despite localized distortions created to accommodate additional ion atoms, the crystalline structure remains practically the same (see illustrations in Fig. S1).⁴³ 

Figure 4 depicts the calculated cell volume versus temperature for pure β-Ga₂O₃. We observe a slight increase in the lattice parameters a, b, c, and volume when increasing the temperature (see Table S1),⁴³ which relates to thermal expansion. The thermal expansion coefficients are derived using linear regression, as described in Eq. (4), while considering the isotropic behavior artificially imposed to the system. The system presented two different behaviors, one between 200 and 1000 K and the other between 1200 and 2000 K, and two values for the thermal expansion coefficient obtained were equal to 2.74 × 10⁻⁶ and 5.99 × 10⁻⁶ K⁻¹, respectively, as shown in Fig. 4. Information regarding the fitting is presented in Table S2.⁴³ Although in this work the system is under an isotropic dilation constrain, these values are comparable to the experimental ones found in the literature, as can be seen in Table S3.^36–38,43

The mechanism involved in the process of ions diffusion in β-Ga₂O₃ is still not well known. Therefore, in this work we analyze the possibility of ionic diffusion in structures with different types of defects and analyze which defects promote diffusion and what mechanisms are involved in the process. Figures S4 and S5⁴³ show the Mean Square Displacement (MSD) of Ga, O, and In ions versus time (t) for pure and with native defects β-Ga₂O₃ for temperatures between 1000 and 2000 K, obtained by Eq. (5). A diffusive pattern occurs when the MSD grows linearly over time, thus no relevant diffusion is observed for most cases in the interval of 5 ns. The only case in which diffusion is more consistent is in Fig. S4(d),⁴³ which shows the diffusion of oxygen ions in β-Ga₂O₃-O_i. Despite Figs. S4(g) and S4(i)⁴³ plus Figs. S5(d), S5(f) and S5(i)⁴³ showed small disturbances in their MSDs of gallium and indium ions, they appear to stem from occasional displacements rather than a diffusive pattern, as will be discussed later.

For a closer look at the oxygen interstitial diffusion, we made a one-dimensional MSD calculation, presented in Fig. 5, and this reveals that the O_i moves preferentially in the y axis direction, inside the channel formed by tetrahedral and octahedral complexes (Fig. S2).⁴³ To confirm this statement, we have tracked the position of each interstitial atom throughout the MD simulation. Figure 6 shows all the positions visited by each O_i, Ga_i and In_i interstitial atoms during the MD simulations at 1500 K. Figure 6(a) shows that the interstitial oxygen really diffuses through the structure preferentially in the y direction. From Fig. 6(b), one can notice that the interstitial oxygen also assumes the lattice positions in the structure. To investigate this behavior, we analyzed the mean residence time, which is defined as the time span during which the oxygen remains in its lattice position. Figure S6⁴³ shows the mean residence time versus temperature, indicating that the oxygen can vacate its lattice position while exchanging with the interstitial oxygen. Figures 6(c) and 6(d) show the snapshot for interstitial indium and Figs. 6(e) and 6(f) for interstitial gallium. From these images, it is evident that these structures do not present any diffusive behavior within the duration of the simulation.

The activation energy (E_a) of self-diffusion can be understood as the energy barrier necessary for diffusion to occur. High E_a means low diffusion, otherwise low E_a leads to high diffusion. Activation energy is related to the diffusion coefficient through Eq. (7), which was obtained from linear regression of the calculated data. Figure 9 shows the Arrhenius plots and the diffusion activation energy (E_a) for β-Ga₂O₃-O_i. The obtained activation energy was 0.94 eV. Table IV shows the oxygen diffusion coefficient values found, and we can see that the Di values are small when compared to other materials such as some perovskites,³⁹ this indicates a slow dynamic, for this reason we carried out simulations for temperatures above 1000 K, with the aim of accelerating the process of diffusion and be able to better understand the associated mechanisms. However, from the fitted values for Eq. (7) we can estimate the diffusion coefficient at 300 K, which was 1.2 × 10⁻¹⁶ Ų/ps, an extremely low value, which suggests that there is negligible oxygen ions diffusion at room temperature, even with O_i native defects. Figure S7⁴³ shows the behavior of the diffusion coefficient with temperature, as well as a plot of Eq. (7) with the values obtained by fitting the Arrhenius plot presented in Fig. 7.

To understand the mechanisms involved in the diffusion process, a careful study was carried out to obtain a complete understanding of the structure and characteristics of the chemical bonds that occur in pure β-Ga₂O₃ and in the defective O_i, In_i and Ga_i in Ga₂O₃ structures. Thus, charge density, electron location function (ELF) and Bader charge were obtained to analyze the chemical bonding in the structures as studied using DFT calculations.

The electron localization function was introduced by Becke and Edgecombe⁴⁰ as a function that indicates the spherically mean probability of finding a given electron with a defined spin in the vicinity of another electron with the same spin. In this way, the ELF determines the position of a pair of opposite spin electrons. For the normalized function 0 ≤ ELF ≤ 1, the value close to 1 means a high probability of finding the electron localization, while when close to 0 a low probability is expected.

In Fig. 8 we show the ELF plots; the blue regions mean low probability for localized electrons while red means high probability. For pure β-Ga₂O₃, shown in Fig. 8(a), we observe mostly ionic bonds since the electrons are only localized around oxygen atoms It's observed an ionic bond, with the transference of electrons of Ga atoms to O atoms, that present high function ELF in your vicinity.⁴¹ For O_i, as shown in Fig. 8(b), we see that the interstitial oxygen breaks the bond between Ga₂₄ and O₂₉. It also forms 2 bonds, one ionic with Ga₂₄ and the other covalent with O₂₉. Figure 8(c), for Ga_i and Fig. 8(d) for In_i show very similar results. In them we can see that the interstitial atom of Ga or In forms 4 bonds with the neighbor oxygens, two in the (010) Miller plane and two in the (001) plane, all of them of ionic character. It is also notable that although they create distortions in the crystal, they do not break any bonds in the original lattice.

With these results we see that interstitial cations adhere more strongly to the crystal lattice, as they form four ionic bonds, while interstitial oxygen makes only two, requiring less energy to disconnect from the lattice and diffuse. The oxygen close to O_i is also less trapped in the crystalline lattice, which makes it possible to exchange places with the interstitial atoms in virtue of thermal effects.

The results for the Bader charges⁴² for oxygen are depicted in Fig. 9(a). We note that the oxygen atoms assume a negative charge and taking the pure β-Ga₂O₃ structure as a reference, we can also note that only the interstitial oxygen and the closest oxygen O (29) in the O_i structure undergo significant variations in their charge states. This corroborates the fact that the two oxygens form a covalent bond and are weakly linked to the lattice. In Fig. 9(b) we see that the Ga and In atoms assume a positive charge and that only the structures with Ga_i and In_i undergo significant changes in the charge states, mainly for the interstitial atom. This allows us to conclude that even though the interstitial cation makes four covalent bonds with the structure, the exchange of charges is less intense and therefore eventual displacements can occur with increasing temperature.

Based on ab initio DFT and MD simulations, this work presents the structural properties for intrinsic defects in indium-alloyed β-Ga₂O₃ structures. MD simulations showed the effects of temperature in the crystal structures as well as in the ion diffusion coefficient; the crystalline structure remained stable, in the range temperature studied (200–2000 K). Furthermore, the presence of interstitial oxygen has been shown to improve the diffusion of oxygen ions, being the only structure that presented a consistent diffusion (for T > 1000 K) due to fewer connections with the crystal lattice compared to other interstitial atoms. Also, it was observed that the diffusion of oxygen ions is greater in the [010]-direction, because of the channel created within the structure by the octahedral and tetrahedral sites in this direction. In summary, the results presented in this work provide evidence for discovering new applications for β-Ga₂O₃ that depend on electrochemical processes. It also helps to understand the influence of native defects and impurities involving indium in In-Ga₂O₃ crystals.

We would like to acknowledge the Army Research Office for funding under Grant Number W911NF-20-1-0298. Calculations were done using the High Performance Computing LEAP (Learning, Exploration, Analysis, and Processing) cluster at Texas State University. D.D.B thanks the Center for Computational Engineering and Sciences (FAPESP/CEPID Grant # 2013/08293-7) for the computational resources.

The authors have no conflicts to disclose.

Nathan Rabelo Martins: Data curation (equal); Formal analysis (equal); Investigation (equal); Methodology (equal); Software (equal); Validation (equal); Visualization (equal); Writing – original draft (equal). Luiz Augusto Ferreira de Campos Viana: Data curation (equal); Formal analysis (equal); Investigation (equal); Methodology (equal); Writing – original draft (equal). Alan Antônio das Graças Santos: Data curation (equal); Formal analysis (equal); Investigation (equal); Methodology (equal); Writing – original draft (equal). Daiane Damasceno Borges: Conceptualization (equal); Data curation (supporting); Formal analysis (supporting); Funding acquisition (equal); Investigation (supporting); Methodology (supporting); Project administration (equal); Resources (equal); Software (equal); Supervision (equal); Writing – original draft (equal). Eric Welch: Data curation (equal); Formal analysis (equal); Investigation (equal); Methodology (equal); Writing – original draft (equal). Pablo Damasceno Borges: Conceptualization (lead); Formal analysis (supporting); Funding acquisition (supporting); Investigation (equal); Methodology (equal); Project administration (equal); Resources (equal); Software (equal); Supervision (equal); Validation (equal); Writing – original draft (equal). Luisa Scolfaro: Conceptualization (equal); Formal analysis (equal); Funding acquisition (equal); Project administration (equal); Supervision (equal).

The data that support the findings of this study are available upon request from the authors.