Wide band gap semiconductors like gallium oxide are promising materials for high-power optoelectronic device applications. We show here a combined density functional theory and molecular dynamics study of diffusion pathways for different defects in *β*-Ga_{2}O_{3}. Molecular dynamics simulations result in a smaller equilibrium volume compared to density functional theory, but the overall lattice remains relatively unchanged even with the inclusion of defects, outside of the local distortions that occur to accommodate the presence of a defect. Slight thermal expansion occurs with elevated temperature and a combination of electron localization function and Bader charge analysis reveals that the oxygen interstitial is the most mobile defect as temperature is increased. However, interstitial cations may diffuse at elevated temperature due to a relatively small amount of charge transfer between the defect and lattice. The mobile oxygen defects are shown to increase the mobility of oxygen ions from the lattice, which can be beneficial for electrochemical applications when controlled through annealing processes.

## I. INTRODUCTION

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,^{1} which make them suitable for high-power, high-chemically stable environments.^{2,3} Specifically, the monoclinic beta phase of gallium oxide, *β*-Ga_{2}O_{3}, 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^{2} V^{−1} s^{1}, has shown immense promise over the last decade of research. Also, numerous dopants have been studied that show the tunability of *β*-Ga_{2}O_{3} through defect engineering.^{4,5} Recently, we have shown that indium is a stable isoelectronic dopant for gallium in *β*-Ga_{2}O_{3},^{6} 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_{2}O_{3}, resulted in an increase in n-type conduction, which was validated through first principles calculations. This highlights the versatility of *β*-Ga_{2}O_{3} 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.^{7} Defects in Ga_{2}O_{3} have been studied for common point defects like Si using hybrid density functional theory (DFT) which highlights the advancement of *ab-initio* simulations.^{8} Our previous study also using hybrid DFT^{6} 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.^{9} 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_{2}O_{3} 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+\delta 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+ \delta 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_{2}O_{3}. 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.

## II. COMPUTATIONAL METHODS

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_{Ga}In_{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.

^{3}. The thermalization phase comprised three stages, each lasting 500 picoseconds (ps): the first stage in the microcanonical ensemble (NVE), the second stage in the canonical ensemble (NVT), and the third stage in the isobaric and isothermal ensemble (NPT) with the isotropic pressure control.

^{17}The chosen simulation times have been proven to be sufficient to reach equilibrium, since the energies rapidly converged. All statistical data have been gathered during the production phase, which lasted for 5000 ps, under NPT conditions, employing a Nosé-Hoover thermostat and barostat.

^{18}A time step of 1 femtosecond (fs) has been used for both phases with the Ewald summation for long-range electrostatic interactions and the Verlet algorithm for atomic motions.

^{19}Short-range interactions between atoms have been described using the Buckingham potential within the Born model for ionic solids.

^{20–22}Equation (1) defines this potential, where

*r*represents the distance between ions

*i*and

*j*, and

*A*,

*ρ*, and

*C*are empirically fit parameters, such that

*α*) is represented by the core shell model, where the ions are divided into two kinds: one that represents the nuclei with positive charge and practically all the atomic mass (∼90%), and the valence electrons that have a negative charge [

*Y*in Eq. (2)] and only about 10% of the mass. In practice, we use 90% of the mass in the core and 10% of the mass in the shell. Both parts are connected by a harmonic spring of elastic constant

*k*. The relation is presented in Eq. (2), such that

The interatomic interaction is described in pairs. The self-interaction between O^{2−} ions and the interaction between O^{2−} and Ga^{3+} have been used to structurally describe Ga_{2}O_{3}.^{23} The In^{3+} – O^{2−} 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_{2}O_{3}. 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.^{28}

Short range interaction . | Core shell model . | |||||
---|---|---|---|---|---|---|

Interaction . | A (eV)
. | ρ (Å)
. | C (eV⋅Å^{6})
. | . | Y (e)
. | k (eV⋅Å^{−2})
. |

Ga^{3+} – O^{2}^{−} | 1625.72 | 0.3019 | 0.0 | – | – | – |

In^{3+} – O^{2}^{−} | 1949.49 | 0.3187 | 0.0 | – | – | – |

O^{2}^{−} – O^{2}^{−} | 9547.96 | 0.2192 | 32.00 | O^{2}^{−} | −2.04 | 6.30 |

Short range interaction . | Core shell model . | |||||
---|---|---|---|---|---|---|

Interaction . | A (eV)
. | ρ (Å)
. | C (eV⋅Å^{6})
. | . | Y (e)
. | k (eV⋅Å^{−2})
. |

Ga^{3+} – O^{2}^{−} | 1625.72 | 0.3019 | 0.0 | – | – | – |

In^{3+} – O^{2}^{−} | 1949.49 | 0.3187 | 0.0 | – | – | – |

O^{2}^{−} – O^{2}^{−} | 9547.96 | 0.2192 | 32.00 | O^{2}^{−} | −2.04 | 6.30 |

*n*(

*r*) is the average number of particles per shell within a radius

*r*from the reference particle and

*ρ*is the particle density over the entire system:

^{16}

*α*) is intrinsically connected to the lattice parameters and can be computed using Eq. (4) for an isotropic system,

^{29}where

*V*represents the lattice volume,

*V*

_{0}denotes the reference lattice volume, and

*T*stands for temperature, such that

*β*-Ga

_{2}O

_{3}structures, we compute the mean square displacement (MSD) of the atoms during the production phase of the MD simulations from the Eq. (5), where $ \u27e8 r 2 ( t ) \u27e9$ represents the time-dependent MSD,

*t*is the time,

*N*is the number of particles and $ r \u2192$ is the position vector. The MSD is connected to molecular self-diffusion via the Stokes-Einstein equation shown in Eq. (6), where

*d*is the dimensionality. The self-diffusion coefficient (

*D*) for species

_{i}*i*is determined by the gradient of the MSD plot using Eq. (6):

*E*, through Eq. (7), where

_{a}*D*

_{0}represents the temperature-independent pre-exponential factor,

*k*

_{B}is the Boltzmann constant, and

*T*is the temperature. The value of

*E*is consequently derived from the gradient of ln(

_{a}*D*) plotted against 1/

_{i}*T*in Arrhenius plots.

In the DFT calculations we used the generalized gradient approximation (GGA) within the Perdew-Burke-Erzenhof (PBE) framework^{30} for the electronic exchange-correlation potential. The following valence electronic distributions were used for each PAW pseudopotential: Ga – 3d^{10} 4s^{2} 4p^{1}, In – 4d^{10} 5s^{2} 5p^{1}, and O – 2s^{2} 2p^{4}. 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\xd71\xd71$ Monkhorst–Pack^{31} 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),^{43} 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.

## III. RESULTS

### A. Structural properties

The structural properties of both pure *β*-Ga_{2}O_{3} 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_{2}O_{3} molecule. The minimum energy $( E 0 \u2217)$ 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.^{34} 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.^{43}

Structure . | DFT . | MD . | ||||||
---|---|---|---|---|---|---|---|---|

Volume (Å^{3})
. | a (Å)
. | b (Å)
. | c (Å)
. | Volume (Å^{3})
. | a (Å)
. | b (Å)
. | c (Å)
. | |

Pure | 212.0 | 12.295 | 3.052 | 5.817 | 208.8 | 12.232 | 3.036 | 5.787 |

O_{i} | 213.6 | 12.327 | 3.060 | 5.832 | 209.2 | 12.240 | 3.038 | 5.791 |

V_{O} | 212.5 | 12.305 | 3.054 | 5.822 | 208.0 | 12.217 | 3.033 | 5.780 |

Ga_{i} | 216.8 | 12.387 | 3.075 | 5.861 | 211.7 | 12.289 | 3.051 | 5.814 |

V_{Ga} | 212.1 | 12.297 | 3.052 | 5.818 | 205.7 | 12.172 | 3.021 | 5.759 |

In_{i} | 214.5 | 12.343 | 3.064 | 5.840 | 213.6 | 12.326 | 3.060 | 5.832 |

In_{Ga} | 212.0 | 12.295 | 3.052 | 5.817 | 210.3 | 12.262 | 3.044 | 5.802 |

In_{Ga}In_{i} | 216.6 | 12.384 | 3.074 | 5.859 | 214.0 | 12.332 | 3.061 | 5.835 |

Structure . | DFT . | MD . | ||||||
---|---|---|---|---|---|---|---|---|

Volume (Å^{3})
. | a (Å)
. | b (Å)
. | c (Å)
. | Volume (Å^{3})
. | a (Å)
. | b (Å)
. | c (Å)
. | |

Pure | 212.0 | 12.295 | 3.052 | 5.817 | 208.8 | 12.232 | 3.036 | 5.787 |

O_{i} | 213.6 | 12.327 | 3.060 | 5.832 | 209.2 | 12.240 | 3.038 | 5.791 |

V_{O} | 212.5 | 12.305 | 3.054 | 5.822 | 208.0 | 12.217 | 3.033 | 5.780 |

Ga_{i} | 216.8 | 12.387 | 3.075 | 5.861 | 211.7 | 12.289 | 3.051 | 5.814 |

V_{Ga} | 212.1 | 12.297 | 3.052 | 5.818 | 205.7 | 12.172 | 3.021 | 5.759 |

In_{i} | 214.5 | 12.343 | 3.064 | 5.840 | 213.6 | 12.326 | 3.060 | 5.832 |

In_{Ga} | 212.0 | 12.295 | 3.052 | 5.817 | 210.3 | 12.262 | 3.044 | 5.802 |

In_{Ga}In_{i} | 216.6 | 12.384 | 3.074 | 5.859 | 214.0 | 12.332 | 3.061 | 5.835 |

The bond angles and bond lengths of the β-Ga_{2}O_{3} – 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.^{43} At very low temperatures, one can differentiate the peaks associated to tetrahedral angles at interval $ 100 \xb0\u2212 120 \xb0$ and to octahedral angles at intervals $\u223c 80 \xb0\u2212 100 \xb0$ and $\u223c 160 \xb0\u2212 180 \xb0$ [see Fig. S3(a)].^{43} These results are coherent with the ADF obtained by DFT calculations, highlighting that both methods well predict the octahedral and tetrahedral bond angles.

. | . | Ga-O (Å) . |
---|---|---|

DFT | Ga-O (1) | 1.83–1.87 |

Ga-O (2) | 1.93–2.09 | |

MD | Ga-O (1) | 1.78–1.88 |

Ga-O (2) | 1.91–2.10 | |

Exp.^{35} | Ga-O (1) | 1.80–1.85 |

Ga-O (2) | 1.95–2.08 |

. | . | Ga-O (Å) . |
---|---|---|

DFT | Ga-O (1) | 1.83–1.87 |

Ga-O (2) | 1.93–2.09 | |

MD | Ga-O (1) | 1.78–1.88 |

Ga-O (2) | 1.91–2.10 | |

Exp.^{35} | Ga-O (1) | 1.80–1.85 |

Ga-O (2) | 1.95–2.08 |

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).^{43}

Figure 4 depicts the calculated cell volume versus temperature for pure *β*-Ga_{2}O_{3}. We observe a slight increase in the lattice parameters *a*, *b*, *c*, and volume when increasing the temperature (see Table S1),^{43} 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^{−6} and 5.99 × 10^{−6} K^{−1}, respectively, as shown in Fig. 4. Information regarding the fitting is presented in Table S2.^{43} 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}

### B. Ions-diffusion

The mechanism involved in the process of ions diffusion in *β*-Ga_{2}O_{3} 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^{43} show the Mean Square Displacement (MSD) of Ga, O, and In ions versus time (*t*) for pure and with native defects *β*-Ga_{2}O_{3} 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),^{43} which shows the diffusion of oxygen ions in *β*-Ga_{2}O_{3}-O_{i}. Despite Figs. S4(g) and S4(i)^{43} plus Figs. S5(d), S5(f) and S5(i)^{43} 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).^{43} 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^{43} 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*means low diffusion, otherwise low

_{a}*E*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 (

_{a}*E*) for

_{a}*β*-Ga

_{2}O

_{3}-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,

^{39}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

^{−16}Å

^{2}/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

^{43}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.

### C. Electron localization function and bader charge

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_{2}O_{3} and in the defective O_{i}, In_{i} and Ga_{i} in Ga_{2}O_{3} 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^{40} 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_{2}O_{3}, 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.^{41} For O_{i}, as shown in Fig. 8(b), we see that the interstitial oxygen breaks the bond between Ga_{24} and O_{29}. It also forms 2 bonds, one ionic with Ga_{24} and the other covalent with O_{29}. 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^{42} for oxygen are depicted in Fig. 9(a). We note that the oxygen atoms assume a negative charge and taking the pure *β*-Ga_{2}O_{3} 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.

## IV. CONCLUSIONS

Based on ab initio DFT and MD simulations, this work presents the structural properties for intrinsic defects in indium-alloyed *β*-Ga_{2}O_{3} 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_{2}O_{3} that depend on electrochemical processes. It also helps to understand the influence of native defects and impurities involving indium in In-Ga_{2}O_{3} crystals.

**ACKNOWLEDGMENTS**

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.

## AUTHOR DECLARATIONS

### Conflict of Interest

The authors have no conflicts to disclose.

### Author Contributions

**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).

## DATA AVAILABILITY

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

## REFERENCES

*Encyclopedia of Condensed Matter Physics*

*Encyclopedia of Condensed Matter Physics*

_{2}O

_{3}crystal structures, structural angles, MSD results, residence time, oxygen diffusion coefficient, volume and lattice parameters for different temperatures and thermal expansion coefficient.