Molecular Dynamics Approach to the Physical Mixture of In2O3 and ZrO2: Defect Formation and Ionic Diffusion

Recent research on the use of physical mixtures In2O3-ZrO2 has raised interesting questions as to how their combination enhances catalytic activity and selectivity. Specifically, the relationship between oxygen diffusion and defect formation and the epitaxial tension in the mixture should be further investigated. In this study, we aim to clarify some of these relationships through a molecular dynamics approach. Various potentials for the two oxides are compared and selected to describe the physical mixture of In2O3 and ZrO2. Different configurations of each single crystal and their physical mixture are simulated, and oxygen defect formation and diffusion are measured and compared. Significant oxygen defect formation is found in both crystals. In2O3 seems to be stabilized by the mixture, while ZrO2 is destabilized. Similar results were found for the ZrO2 doping with In and ln2O3 doping with Zr. The results explain the high activity and selectivity catalyst activity of the mixture for the production of isobutylene from ethanol.


Introduction
Indium oxide is a versatile material that arouses interest in several applications, such as optoelectronics [1] and catalysis [2,3]. Recent researches have shown that the interest in this material for catalysis is due partially to the strong capacity of the cubic phase of In 2 O 3 to form oxygen vacancies, especially in the presence of monoclinic ZrO 2 [4][5][6][7]. The causes of this phenomenon are still not well understood and have been the target of theoretical and experimental studies for its elucidation. These materials seem to form an epitaxial alignment that favor defects' formation, such as. oxygen vacancies, providing an interesting catalytic behavior [6,7].
Further, the oxygen dynamics of ZrO 2 and In 2 O 3 favor some reactions of industrial interest, like the syntheses of isobutene and acetone from ethanol [5,6] and also methanol and methane from CO 2 [8,9]. These reactions depend on the oxygen lattice mobility, according to the Mars Van Krevelen mechanism [10].
The optimization and development of new In 2 O 3 catalysts depend on understanding the phenomenon of point defects formation and mobility. In this sense, computer simulations help predict the behavior and properties of these materials. [11] Some authors have used Density Functional Theory (DFT) to investigate the morphology and stability of different faces of cubic In 2 O 3 [12,13]. Walsh et al. [14] found that the stability of In 2 O 3 surfaces can be described in the following order: (111) > (110) > (100). In addition, Zhang et al. [13] verified that the face (111) of In 2 O 3 presents greater stability under oxidizing conditions, while the face (100) is more stable under oxygen-lean conditions. Through theoretical and experimental study, Ziemba et al. [15] used DFT calculations to simulate experimental Raman results and investigate the oxygen dynamics for cubic In 2 O 3 under oxidizing and reducing conditions. From the calculations, they determined that the energy of formation 2 of 18 of an oxygen vacancy in the bulk of an In 2 O 3 unit cell was 3.40 eV, consistent with the values found in the literature [1,16]. In this study, the authors found that, under reducing conditions at temperatures above 120 • C, oxygen migrates from the bulk to the surface until it reaches an equilibrium condition, which is more favorable than maintaining oxygen vacancies on the surface of the oxide [15]. Nørskov et al. [17] calculated the energy of formation of oxygen vacancies on the surface (111) of In 2 O 3 considering a supported catalyst model represented by In 2 O 3 (111) layers on ZrO 2 (111). The thickness of the In 2 O 3 layer was varied (from 5 to 2 layers) to simulate different conditions of dispersion of In 2 O 3 on the support. They found that reducing the In 2 O 3 layer thickness increased the formation energy for oxygen vacancies, suggesting that increasing this oxide dispersion makes defect formation more difficult. For the more dispersed model (2 layers), the vacancy formation energy was 0 eV, while for the 5-layer model, the calculated value was −0.1 eV. The small vacancy formation energies show a relative stabilization of the In 2 O 3 lattice in the physical mixture, resulting in a minor defect formation and a higher oxygen diffusion activation energy.
Molecular dynamics also contributed to the understanding and better description of materials but used empirical interatomic potentials to describe the structure of pure cubic In 2 O 3 [1] or doped with elements as Sn [18,19]. Such potentials have often been used in molecular dynamics simulations to investigate the formation of defects as vacancies, Frenkel and Schottky [16,19,20], and to assess the stability of different faces of cubic In 2 O 3 [14]. Warschkow et al. [18] proposed a Buckingham potential for In 2 O 3 , based on studies by Bush et al. [21]. Although this potential could provide adequate values for the structural properties of cubic In 2 O 3 , it barely fitted the dielectric parameters. Walsh et al. [1] reported a more accurate potential that simulates both the structural and physical properties of the material. Furthermore, such potential proved adequate for studies of intrinsic defects such as vacancies, Schottky, and Frenkel in this oxide.
Recently, Hou et al. [19] developed a set of paired interatomic potentials, combining Buckingham, Lennard-Jones, Polynomial, Polynomial Harmonic Potentials, and Shell models. This new potential proved accurate to simulate structural and physical properties and understand the properties resulting from the formation of defects in the oxide.
For ZrO 2 , most reported interatomic potentials are more appropriate for cubic or tetragonal phases and may present some deviations when used for simulations involving the monoclinic phase. Schelling et al. [22] developed a Buckingham potential to investigate the phenomenon of transition from cubic and tetragonal phases to cubic yttria-stabilized zirconia (YSZ). Kilo et al. [23] investigated oxygen diffusion in the same material over a finite temperature range, modifying the literature's potentials. Lau & Dunlap [24] used diluted concentrations of Y 2 O 3 in ZrO 2 and described a suitable potential for the temperature range of 300-1400 K, using the rigid ion model approximation. Through a different line of study, Duin et al. [25] developed a ReaxFF reactive force field to search into the intergranular diffusion of perovskites like Y-doped BaZrO 3 (BYZ).
Despite these studies, there is not enough literature describing the interaction of In 2 O 3 with ZrO 2 based on molecular dynamics. Therefore, some questions are raised: would it be possible to use these potentials to describe the behavior of In 2 O 3 and ZrO 2 catalysts? Furthermore, how the presence of ZrO 2 could affect the stability and oxygen mobility of In 2 O 3 ?
This work uses molecular dynamics and empirical potentials to evaluate the interaction between In 2 O 3 and ZrO 2 oxides in proximity, the defect formation in the In 2 O 3 oxide, and the diffusion of ionic species in mixtures of these oxides. The present results are discussed considering the known catalytic properties of the In 2 O 3 /ZrO 2 system and its effectiveness in catalyzing the isobutylene synthesis from ethanol.

Methodology
The used interatomic potentials are a sum of a Coulomb term describing electrostatic interactions; and a Buckingham potential, which empirically describes the quantum interactions of the Pauli repulsion and the van der Waals attraction [26].
where r ij is the distance between atoms i and j, q i and q j are their charges. The parameters A, ρ, and C were empirically derived. The simulations were performed using the LAMMPS software, version 3 march 2020 [27] with the metal style. The cutoff radius was 10 Å for the Buckingham and 25 Å for the Coulomb interactions [1]. The long-range interactions were calculated using Ewald summations and long-range corrections were not used. The potentials' parameters are presented in Table 1 with their respective references. These potentials presented the best match for oxygen-oxygen interactions and thus were best suited for combination. Further analysis on potential selection is presented in other sections of this work.

In 2 O 3 Model Verification and Selection
This oxide's crystal structure is described by a unitary cell containing 80 atoms and pertains to the Ia-3 spatial group [1]. The input files with the atomic positions for the rigid ion framework were obtained through a Python code using the Atomic Simulation Environment (ASE) library [28]. A simple Python algorithm produced the input files with the atomic positions for the core-shell framework from the rigid ion files.
Preliminary experiments were conducted on two existing potentials to determine which better describes the crystal. Furthermore, for each potential, a rigid ion and a core-shell case were investigated to check for a significant variation.
These experiments consisted of building a lattice parameter versus energy graph for each potential choice using a 4 × 4 × 4 In 2 O 3 supercell. Runs for different lattice parameters were performed in NVT ensembles at room temperature for 4 ps using 0.2 fs time steps. The resultant medium potential energy per atom was registered ( Figure S1). After that, a third-degree polynomial equation was fitted to the simulated points to help identify the local minimum which corresponds to the equilibrium lattice parameter. The bulk modulus, B, was extracted from the fitted polynomial through the following equation [29,30]: V is the volume at constant temperature, P is the pressure, N is the number of atoms in the unit cell, a 0 is the equilibrium lattice parameter, a is the lattice parameter, and E is the potential energy described by the fitted third-degree polynomial equation.
For the linear expansion coefficient calculation, an NVT run with 4 ps of equilibration plus 0.2 ps steps of measurement was performed for different temperatures using a time step of 0.2 fs. The lattice parameter for each temperature was obtained and used to build a graph of elongation versus the difference in temperature to determine the linear expansion coefficient ( Figure S2). The results were compared with experimental data and the original calculated values by their authors [1,18]. The adequate potential and case-core-shell or rigid ion-were chosen.

Defect Energy
The energies of various defects on the In 2 O 3 crystal and its lattice energy were calculated and compared with Walsh et al. [1]. The rigid ion and the core-shell models were used, measured, and compared. Two different supercell sizes were considered: 3 × 3 × 3 and 4 × 4 × 4.
For each model and each supercell size, input files were crafted depicting the following cases: crystal with no defects; 8b indium atom vacancy; 24d indium atom vacancy; 48e oxygen atom vacancy; 8a interstitial indium atom; 16c interstitial indium atom; 24d interstitial indium atom; 16c interstitial oxygen atom; 24d interstitial oxygen atom; and 8a interstitial oxygen atom. The letters correspond to the Wyckoff positions [31,32]. Ten different input files with randomly positioned defects were generated for each defect type and interatomic potential. A minimization run was executed for each file, and the total energy average was calculated for each defect type. These values were subtracted from the total energy of the defect-free lattice with the same ion model and supercell size.
For the Frenkel and Schottky defect energies, a calculation using the lower energy vacancies and interstitial atoms energies [1] was performed: For the lattice energies, the total potential energy of the In 2 O 3 crystal without defects was divided by the number of atoms, then multiplied by 5, the number of atoms in the In 2 O 3 formula.

ZrO 2 Model Selection
Bandura et al. compared different interatomic potentials to describe zirconia cubic, monoclinic and tetragonal phases [33]. Three interatomic potentials using the Buckingham potential were chosen for further analysis to allow compatibility with Walsh's potential for In 2 O 3 , especially regarding the O −2 -O −2 interaction, since the crystals were expected to exchange O −2 ions. These candidate potentials did not have the monoclinic phase as their most stable. Instead presented an artificial rutile-like structure as their most stable phase. Stability tests were performed for these 3 potentials using 4 × 4 × 4 ZrO 2 supercells, and NPT runs for 80 ps with 1fs time steps for stabilization plus 60 ps of measurement. The radial distribution functions of their zirconium atoms present in the stabilized phase were calculated and compared to the RDFs of other ZrO 2 phases: artificial rutile-like, monoclinic and tetragonal. The cubic structure was not considered for it was the least stable phase for all 3 candidate potentials.
To further investigate compatibility between ZrO 2 and In 2 O 3 potentials, a graph of potential energy versus distance between oxygen ions was determined and compared.

Defective In 2 O 3 and ZrO 2
For the defective In 2 O 3 experiments, 4 types of defects were investigated: Schottky defects, oxygen Frenkel defects, oxygen vacancies, and substitutional defects with Zr atoms occupying In sites. For defective ZrO 2 experiments, Schottky defects, Frenkel defects, oxygen vacancies, and substitutional defects with In atoms occupying Zr sites. For each defect, input files representing different initial conditions were prepared using the ASE library. A 3 × 3 × 3 supercell was used to investigate defects in In 2 O 3 . 13 Schottky defects were created by adding 26 In vacancies and 39 O vacancies to the base supercell. For Frenkel defects, 65 oxygen Frenkel defects were added to the supercell. Sixty-five oxygen atoms were eliminated from the supercell for pure oxygen vacancy defects. So the total number of vacancies was always 65, which is 3% of the total number of atoms in the perfect crystal.
The effect of Zr doping in the In 2 O 3 crystal was investigated with two different files containing Zr doped In 2 O 3 . One had 78 Zr atoms, and another had 195 Zr atoms, substituted at random 8a In positions, as this is the indium atom site with the lowest vacancy formation energy. Vacancies for indium atoms were generated to equilibrate the crystals' charge. Thus, the files with 78 Zr atoms contained the same number of In vacancies of the Schottky defect file, 26; and the files with 195 Zr atoms contained the same number of total vacancies of the Schottky defect file 65.
Similar files were generated to investigate defects on ZrO 2 using a 5 × 5 × 5 base supercell. So to create 15 Schottky defects, 15 Zr vacancies and 30 O vacancies were created in the base supercell. For Frenkel defects, 45 oxygen Frenkel defects were added to the supercell. For pure oxygen vacancies, 45 oxygen atoms were taken from the supercell. So for each defect type, the total number of vacancies was 45, that is, 3% of the total number of atoms in the perfect crystal.
The effect of In doping in the ZrO 2 crystal was investigated with two different files containing In using the 5 × 5 × 5 base supercell. One had 60 In atoms, and another had 90 In atoms, substituted at random Zr positions, as this crystal only possesses one type of site for Zr atoms. Vacancies for oxygen atoms were generated to equilibrate the crystals' charge. Thus, the files with 60 In atoms contained the same number of oxygen vacancies of the Schottky defect file, 30; and the files with 90 In atoms contained the same number of total vacancies of the Schottky defect file, 45.
The defective crystals were built using the rigid model allowing an easier integration between different potentials. They were subject to a cg style hydrostatic pressure minimization at zero pressure, followed by a 50 ps equilibration hydrostatic NPT run at 1.01325 bar and a 1000 ps hydrostatic pressure NPT measuring run at various temperatures and 1.01325 bar. The simulation time and the time step of 1 fs were chosen to be similar to values used to simulate oxygen diffusion in Kilo et al. [23]. The temperatures were 850, 900, 950, and 1000 • C. The O −2 -O −2 potential used for the doping simulations was that of the host crystal instead of Walsh's potential for all O −2 -O −2 interactions employed in the description of the physical mixtures. The mean square displacement of each atom type was computed during these experiments.

In 2 O 3 -ZrO 2 Physical Mixture
Input files describing physical mixtures arrangements of layers of In 2 O 3 and ZrO 2 were crafted using the ASE library to generate slabs of each crystal in different planes and then processing these slabs with Python scripts to appropriately position the ZrO 2 slabs on top of a larger In 2 O 3 slab. As the ASE code outputted diamond-shaped slabs, this shape was conserved for the input data files. The used slabs were 3 × 3 × 3 In 2 O 3 (111), 5 × 5 × 4 ZrO 2 (111), 3 × 7 × 5 ZrO 2 (102), and 3 × 6 × 6 ZrO 2 (211). The 3 × 3 × 3 In 2 O 3 (111) slab dimensions were taken as a starting point to accommodate the ZrO 2 chosen surface, and the ZrO 2 slabs dimensions for the three different planes were chosen to cover most of the In 2 O 3 (111) slab surface; their length and width not being larger than the In 2 O 3 (111) slab surface; and the similarity in area and number of atoms between them. This last item is important because it allows a reasonable comparison between the different studied interfaces. The dimensions and number of atoms of the ZrO 2 slabs are as follows: 5 × 5 × 4 ZrO 2 (111)-1,180,372 A 2 and 1200 atoms; 3 × 7 × 5 ZrO 2 (102)-1266,318 A 2 and 1260 atoms; 3 × 6 × 6 ZrO 2 (211)-1210,655 A and 1296 atoms. Figure 1a,b shows an example of the initial file used to simulate the physical mixture system with an In 2 O 3 base. It is important to note that the generated input files have the In 2 O 3 crystal touching the lateral edges of the periodic box, while this does not happen for the ZrO 2 crystal present in these files. As this study uses periodic boundaries, which mirror the crystal structure in all directions, this leads to the formation of a continuous layer of the base crystal, leading to it experiencing a higher stabilization. Figure 1c shows the same system with the periodic structure replicated to show the formation of a continuous layer for the base oxide. a higher stabilization. Figure 1c shows the same system with the periodic structure replicated to show the formation of a continuous layer for the base oxide. For each In2O3/ZrO2 combination, a twin file containing approximately 3% vacancies due to Schottky defects in each crystal was also crafted.
This file was manipulated into a rectangular box to fit LAMMPS limitations on simulation box angle and could not accept the box angle of the ZrO2 base used.
For comparison's sake, a 6 × 6 × 5 ZrO2(111) slab was crafted as the base with its laterals touching the periodic boundary, and the 3 × 3 × 2 In2O3(111) slab placed on its top. If there are no significant differences between this slab's arrangement and the previous, then the periodic boundary effect to stabilize the base crystal is insignificant. This setup file was manipulated into a rectangular box to fit LAMMPS limitations on simulation box angle and could not accept the box angle of the ZrO2 base used.
All files described were built using the rigid ion model. They were subject to a cg style triclinic minimization at 1.01325 bar, followed by a 50 ps equilibration triclinic NPT run at 1.01325 bar and a 1000 ps measuring triclinic NPT run at 850, 900, 950, 1000, and 1050 °C, and 1.01325 bar. The simulation time and the time step of 1 fs were chosen to be similar to values used to simulate oxygen diffusion in Kilo et al. [25]. During these experiments, the following information was acquired: the potential energy of each atom; the mean square displacement of each atom type, distinguishing the ZrO2 and In2O3 oxygen atoms; and the radial distribution function of each atom type.

In2O3-ZrO2 Interface
Lewis [34] estimated the changes in cation coordination number on the repulsive term of the pair potential. They showed that the ionic radii for the tetrahedral site are about 0.94 the ionic radii for the octahedral site. This difference results in a reduction of 14% of the parameter A in Equation (1). The In2O3 and ZrO2 structures used in the present calculation have the same coordination number, so the same Zr +4 -O −2 and In +3 -O −2 pair potentials may be used inside the In2O3 and ZrO2 structures. Concerning the metallic ions at the interface region, one expects that the coordination number stays mostly the same since the oxygen ions must keep the electric neutrality of the crystal. Therefore, the same Zr +4 -O −2 and In +3 -O −2 pair potentials were used for both crystals and the interface. Following Lewis [34], the O −2 -O −2 interaction is the same throughout the mixture, and the cationcation interaction is purely Coulombic.
Before the mixture simulations, the chosen In2O3 and ZrO2 interfaces were placed at a short distance from each other, and the system was allowed to relax in all directions at 1.01325 bar (1 atm) and the desired temperature with a maximum volume change of 0.001 in any one direction at each interaction. After the equilibration, the timestep was reset, For each In 2 O 3 /ZrO 2 combination, a twin file containing approximately 3% vacancies due to Schottky defects in each crystal was also crafted.
This file was manipulated into a rectangular box to fit LAMMPS limitations on simulation box angle and could not accept the box angle of the ZrO 2 base used.
For comparison's sake, a 6 × 6 × 5 ZrO 2 (111) slab was crafted as the base with its laterals touching the periodic boundary, and the 3 × 3 × 2 In 2 O 3 (111) slab placed on its top. If there are no significant differences between this slab's arrangement and the previous, then the periodic boundary effect to stabilize the base crystal is insignificant. This setup file was manipulated into a rectangular box to fit LAMMPS limitations on simulation box angle and could not accept the box angle of the ZrO 2 base used.
All files described were built using the rigid ion model. They were subject to a cg style triclinic minimization at 1.01325 bar, followed by a 50 ps equilibration triclinic NPT run at 1.01325 bar and a 1000 ps measuring triclinic NPT run at 850, 900, 950, 1000, and 1050 • C, and 1.01325 bar. The simulation time and the time step of 1 fs were chosen to be similar to values used to simulate oxygen diffusion in Kilo et al. [25]. During these experiments, the following information was acquired: the potential energy of each atom; the mean square displacement of each atom type, distinguishing the ZrO 2 and In 2 O 3 oxygen atoms; and the radial distribution function of each atom type.

In 2 O 3 -ZrO 2 Interface
Lewis [34] estimated the changes in cation coordination number on the repulsive term of the pair potential. They showed that the ionic radii for the tetrahedral site are about 0.94 the ionic radii for the octahedral site. This difference results in a reduction of 14% of the parameter A in Equation (1). The In 2 O 3 and ZrO 2 structures used in the present calculation have the same coordination number, so the same Zr +4 -O −2 and In +3 -O −2 pair potentials may be used inside the In 2 O 3 and ZrO 2 structures. Concerning the metallic ions at the interface region, one expects that the coordination number stays mostly the same since the oxygen ions must keep the electric neutrality of the crystal. Therefore, the same Zr +4 -O −2 and In +3 -O −2 pair potentials were used for both crystals and the interface. Following Lewis [34], the O −2 -O −2 interaction is the same throughout the mixture, and the cation-cation interaction is purely Coulombic.
Before the mixture simulations, the chosen In 2 O 3 and ZrO 2 interfaces were placed at a short distance from each other, and the system was allowed to relax in all directions at 1.01325 bar (1 atm) and the desired temperature with a maximum volume change of 0.001 in any one direction at each interaction. After the equilibration, the timestep was reset, and an initial equilibration was run for 50,000 steps before the final simulation was run for 1,000,000 steps.

In 2 O 3 -ZrO 2 Physical Mixture-ZrO 2 Monolayer
ZrO 2 was also simulated as a monolayer over the In 2 O 3 crystal. Three different surface planes over an In 2 O 3 crystal were crafted. The In 2 O 3 dimensions were the same as before but the ZrO 2 height corresponded to a unit cell perpendicular to the chosen plane. The number of atoms for each monolayer is ZrO 2 (111)-300 atoms; ZrO 2 (102)-252 atoms; and ZrO 2 (211)-324 atoms.
For each In 2 O 3 /ZrO 2 (monolayer) combination, a twin file containing approximately 3% vacancies due to Schottky defects in each crystal was also crafted. Figure 2 shows an example of the In 2 O 3 /ZrO 2 (monolayer) combination. and an initial equilibration was run for 50,000 steps before the final simulation was run for 1,000,000 steps.

In2O3-ZrO2 Physical Mixture-ZrO2 Monolayer
ZrO2 was also simulated as a monolayer over the In2O3 crystal. Three different surface planes over an In2O3 crystal were crafted. The In2O3 dimensions were the same as before but the ZrO2 height corresponded to a unit cell perpendicular to the chosen plane. The number of atoms for each monolayer is ZrO2(111)-300 atoms; ZrO2(102)-252 atoms; and ZrO2(211)-324 atoms.
For each In2O3/ZrO2(monolayer) combination, a twin file containing approximately 3% vacancies due to Schottky defects in each crystal was also crafted. Figure 2 shows an example of the In2O3/ZrO2(monolayer) combination.
All of the files were built using the rigid ion model. The files for the monolayers were subjected to the same simulation script used on the bulk physical mixtures.

MSD and Diffusion
The diffusion coefficient was computed from the atomic mean square displacement (MSD) determined every 100 timesteps for each group of atoms. The MSD shows a transient region that was not considered, and the stationary region was fitted with a linear function. The slope was used to calculate the diffusivity as a function of temperature, according to the following equations [35]: where ri(t) is the position of the atom i at the time t, D is the diffusion coefficient, t is the time.
The diffusivity, D, and the activation energy, EA, were obtained according to the following equations: Due to the relatively small equilibration run and low temperatures studied, many MSD curves only stabilized to a linear-like shape at some later point during the measurement run. Thus, the function for determining diffusivity was only fitted during this interval. For all other MSD curves that presented a linear shape from the beginning of the measurement run, the diffusivity function was fitted from steps 200 thousand to 1 million. All of the files were built using the rigid ion model. The files for the monolayers were subjected to the same simulation script used on the bulk physical mixtures.

MSD and Diffusion
The diffusion coefficient was computed from the atomic mean square displacement (MSD) determined every 100 timesteps for each group of atoms. The MSD shows a transient region that was not considered, and the stationary region was fitted with a linear function. The slope was used to calculate the diffusivity as a function of temperature, according to the following equations [35]: where r i (t) is the position of the atom i at the time t, D is the diffusion coefficient, t is the time. The diffusivity, D, and the activation energy, E A , were obtained according to the following equations: Due to the relatively small equilibration run and low temperatures studied, many MSD curves only stabilized to a linear-like shape at some later point during the measurement run. Thus, the function for determining diffusivity was only fitted during this interval. For all other MSD curves that presented a linear shape from the beginning of the measurement run, the diffusivity function was fitted from steps 200 thousand to 1 million.

Defect Formation
The total number of interstitial oxygen atoms was also studied in the simulation runs of defective In 2 O 3 , defective ZrO 2 , and all physical mixtures at the lowest measured temperature, 850 • C. Because many studied cases already contained vacancies from the beginning of the simulation, only oxygen interstitials were measured to determine the oxygen vacancy formation through Frenkel defects in each system in a comparable manner. The time series and time averaging functions of the visualization software OVITO [36] were used to average the number of interstitial atoms over the total number of state 'snapshots' (LAMMPS dumps) taken after the stabilization of the system. The average interstitial count was divided by the total number of oxygen sites in the perfect crystal of the corresponding In 2 O 3 or ZrO 2 slab to normalize the different crystal sizes.
As with the MSD function fitting, the average of interstitial oxygen atoms was taken from step 200 thousand to one million as a rule of thumb. For systems in which MSD stabilized only later in the measurement run, this average was taken in the same interval used for the diffusivity function fitting.

Radial Distribution Function (RDF)
In some simulation runs where the coordination of atoms was measured, radial distribution function (RDF) measures were computed for that atom group through a LAMMPS compute every 200 timesteps. These measurements were then outputted every 20 thousand steps and analyzed through a Python script, producing graphs of RDF vs distance from the atom. The first undulation represents the first coordination shell of the atom, the first neighbors. The area under this first coordination shell is the atom's coordination number. The coordination number of the first coordination shell was taken directly from the LAMMPS compute.

In 2 O 3 Model Verification and Selection
The rigid ion model with Walsh's potential provided the best description for the lattice parameter and the linear expansion coefficient experimental values, as evidenced in Table 2. So it was chosen for the present calculations. Lau & Dunlap [24] had already argued in their work that this approach is appropriate in MD simulations when no electric field is applied to the solid. The choice of the rigid ion model simplifies the integration between different Buckingham potentials. Table 3 shows all considered potential parameters for each crystal.

Rigid Ion Model Results
Lattice parameter (Å) [ Table 4 shows that for both the 3 × 3 × 3 and the 4 × 4 × 4 supercell, the rigid ion model was overall more accurate when compared to the data presented by Walsh. The small differences between our simulations and Walsh's results indicate that the supercell method used by us and the Mott-Littleton approach used by Walsh are comparable, which further validates the proposed rigid ion model.

ZrO 2 Model Selection
The results presented in the literature for the various ZrO 2 potentials [33] were organized to represent only those potentials capable of accurately modeling the ZrO 2 crystal through the rigid-ion model. Table 5 compares candidate potentials for describing the ZrO 2 crystal. The properties evaluated are the lattice parameters a, b, and c; the unit cell angle β; and the bulk modulus. This analysis concluded that the best three candidate potentials were Schelling's, Lau's, and Kilo's. It was found that indeed none of them presents the monoclinic phase under the proposed simulated conditions. In the temperature range 473.15 K to 1273.15 K, both Lau's and Kilo's potential presented a zirconium coordination number of 6, as shown in Figure 3, while Schelling's coordination number for zirconium was 8. Thus, it can be concluded that Lau's and Kilo's potentials developed an orthorhombic, rutile-like cell structure, while Schelling's potential presented a tetragonal structure. Kilo's potential was chosen to describe the ZrO 2 compound. It results in an orthorhombic structure with a small deviation from the true monoclinic structure (less than 10 degrees in beta) and has an O −2 -O −2 potential that resembles Walsh's potential used to describe the In 2 O 3 crystal, as shown in Figure 4. The similarity of the O −2 -O −2 potentials suggests the use of only one potential description for both In 2 O 3 and ZrO 2 phases. When comparing the potential energy for oxygen ions interaction for the three best ZrO 2 potentials and the two In 2 O 3 potentials, it is clear that the potentials of Walsh and Kilo show great compatibility. So, Walsh's potential was chosen to describe anionic interactions in the physical mixtures, as noted in Table 1. It was found that indeed none of them presents the monoclinic phase u posed simulated conditions. In the temperature range 473.15 K to 1273.15 and Kilo's potential presented a zirconium coordination number of 6, as sh 3, while Schelling's coordination number for zirconium was 8. Thus, it can that Lau's and Kilo's potentials developed an orthorhombic, rutile-like while Schelling's potential presented a tetragonal structure. Kilo's potenti to describe the ZrO2 compound. It results in an orthorhombic structure with ation from the true monoclinic structure (less than 10 degrees in beta) and potential that resembles Walsh's potential used to describe the In2O3 crysta     Table 6 shows the activation energy for oxygen diffusion considering diffe fects in both crystals. First, oxygen exhibits higher diffusion activation energy in Z all intrinsic defects, probably due to the smaller interstitial size in this crystal. I observed that the diffusion activation energy increases following the same sequ intrinsic defects for both the zirconia and indium oxide. Schottky has the lowest followed by Frenkel and oxygen vacancies. To understand this phenomenon, OVI used. It was found that oxygen vacancies tended to distribute themselves symm due to their effective repelling charges, which seem to stabilize the lattice and diffusion to an extent. For Frenkel defects, it was observed that interstitial atoms to recombine with vacancies during the equilibration run. This way, by the time th urement was done, there were significantly fewer Frenkel defects, and the crys almost perfect, thus making diffusion more difficult. Doped ZrO2 exhibits the lowest diffusion activation energy for the substituti fects, while doped In2O3 oxygen exhibits the highest, which may be attributed to of defect sites generated by the doping and their effective charges. In doped ZrO ates oxygen vacancies and substitutional indium sites. Oxygen vacancies are m ceptible to aid in diffusion than cationic vacancies due to their effective positive and their matching size with the oxygen atom. Also, the effective negative charg  Table 6 shows the activation energy for oxygen diffusion considering different defects in both crystals. First, oxygen exhibits higher diffusion activation energy in ZrO 2 for all intrinsic defects, probably due to the smaller interstitial size in this crystal. It is also observed that the diffusion activation energy increases following the same sequence of intrinsic defects for both the zirconia and indium oxide. Schottky has the lowest energy, followed by Frenkel and oxygen vacancies. To understand this phenomenon, OVITO was used. It was found that oxygen vacancies tended to distribute themselves symmetrically due to their effective repelling charges, which seem to stabilize the lattice and prevent diffusion to an extent. For Frenkel defects, it was observed that interstitial atoms tended to recombine with vacancies during the equilibration run. This way, by the time the measurement was done, there were significantly fewer Frenkel defects, and the crystal was almost perfect, thus making diffusion more difficult. Doped ZrO 2 exhibits the lowest diffusion activation energy for the substitutional defects, while doped In 2 O 3 oxygen exhibits the highest, which may be attributed to the type of defect sites generated by the doping and their effective charges. In doped ZrO 2 generates oxygen vacancies and substitutional indium sites. Oxygen vacancies are more susceptible to aid in diffusion than cationic vacancies due to their effective positive charge and their matching size with the oxygen atom. Also, the effective negative charge of the substitutional indium atoms repels oxygen, avoiding trapping the anions in this way. Figure 5 shows the path of some interstitial oxygen atoms compared with the substitutional indium atoms position, showing that these cations tend to be avoided by diffusing anions. These phenomena help lower the diffusion activation energy for oxygen in zirconia. However, Zr doped In 2 O 3 creates indium vacancies and substitutional zirconium sites, making oxygen diffusion more energetically costly. The created indium vacancies are larger than oxygen atoms and have an effective charge of −3, making them repel oxygen atoms, thus not favoring their diffusion. Further, substitutional zirconium atoms have an effective charge of +1, attracting and trapping oxygen atoms in their vicinity, and slowing their diffusion. Figure 6 shows the trapping effect of Zr substitutional atoms on oxygen diffusion. oxygen atoms and have an effective charge of -3, making them repel oxygen atoms, thus not favoring their diffusion. Further, substitutional zirconium atoms have an effective charge of +1, attracting and trapping oxygen atoms in their vicinity, and slowing their diffusion. Figure 6 shows the trapping effect of Zr substitutional atoms on oxygen diffusion. The data for oxygen defect formation in Table 7 does not show many interesting patterns. The ZrO2 crystal showed higher overall defect formation. Also, as mentioned before, the recombination of interstitials and vacancies during the equilibration step reduces the number of Frenkel defects by the time of the measurement run. No obvious correlation between diffusion activation energy and defect formation was identified.   oxygen atoms and have an effective charge of -3, making them repel oxygen atoms, thus not favoring their diffusion. Further, substitutional zirconium atoms have an effective charge of +1, attracting and trapping oxygen atoms in their vicinity, and slowing their diffusion. Figure 6 shows the trapping effect of Zr substitutional atoms on oxygen diffusion. The data for oxygen defect formation in Table 7 does not show many interesting patterns. The ZrO2 crystal showed higher overall defect formation. Also, as mentioned before, the recombination of interstitials and vacancies during the equilibration step reduces the number of Frenkel defects by the time of the measurement run. No obvious correlation between diffusion activation energy and defect formation was identified.  The data for oxygen defect formation in Table 7 does not show many interesting patterns. The ZrO 2 crystal showed higher overall defect formation. Also, as mentioned before, the recombination of interstitials and vacancies during the equilibration step reduces the number of Frenkel defects by the time of the measurement run. No obvious correlation between diffusion activation energy and defect formation was identified.  Table 8 shows the oxygen diffusion activation energy with positive values for all attempted combinations. The mean square displacement of cations was very low, and for this reason, their diffusion activation energies were not calculated ( Figure S3). It should also be noted that some calculations have near-zero or even negative oxygen diffusion activation energies due to the very low diffusivity (Figures S4 and S5). It is believed that the computation time was too short or the temperature too low.  Table 8 presents interesting trends. In all physical mixtures, the activation energy for oxygen diffusion in no-defects indium oxide was greater than the value observed for the 3% Schottky non-mixture In 2 O 3 crystal ( Table 6). The physical mixture shows the same qualitative effect on oxygen diffusion activation energy as the doping of these oxides with one another-ZrO 2 elevates oxygen activation energy in In 2 O 3 , and In 2 O 3 lowers the oxygen activation energy in ZrO 2 . The activation energy for oxygen diffusion in no-defect zirconia in the physical mixture was lower than for all intrinsic defects in non-mixture ZrO 2 . However, the lowest energy for all ZrO 2 oxygens was still found for the case of substitutional indium defects in a non-mixture ZrO 2 crystal, which points to an important result.

In 2 O 3 -ZrO 2 Physical Mixture
The mechanical, physical mixtures with the associate epitaxial piling of the ZrO 2 and In 2 O 3 nanocrystals affect oxygen diffusion equivalent to a doping process. It results from either small doping at the two materials interface or, most probably, the lattice tensions caused by the epitaxy between the two crystals. Yamamoto et al. [41] observed from DFT calculations that compressive stress caused a reduction in the oxygen diffusion coefficient in monoclinic and tetragonal ZrO2. These lattice tensions are illustrated in Figures 7 and 8, which show how the interaction between crystals creates patterns of potential energies. For the In 2 O 3 crystal, it seems these tensions stabilize oxygen atoms, inhibiting their diffusion even more than doping with zirconium, which agrees with the calculations from Nørskov et al. [17], that showed that oxygen vacancy formation energy is elevated in In 2 O 3 by its contact with ZrO 2 . For ZrO 2 , the effect seems to be destabilization of the atom, although doping with indium still stimulates diffusion more significantly.  For the ZrO2 oxygens present in bulk ZrO2 over In2O3 and for In2O3 oxygens present in In2O3 over ZrO2, another interesting pattern is the lower diffusion activation energy of the setups without intrinsic defects in comparison to those with them. This behavior also points to the great influence of lattice tensions on oxygen diffusion in the studied physical mixtures. It is argued that for crystals with intrinsic defects, there is more room for lattice tensions to be dissipated through lattice distortion, while in perfect crystals there is not an easy outlet for this accumulated energy apart from atoms jumping lattice sites.
Lastly, activation energies for ZrO2 oxygens diffusion varied little between physical mixtures, presenting values between 0.51 and 0.88 eV, which is likely due to the stronger attraction oxygen atoms have towards zirconium atoms.
The results for defect formation in physical mixtures in Table 9 do not show any obvious correlation between defect formation and diffusion activation energy. Nevertheless, some general statements can be made about them. For all cases, physical mixtures had higher defect formation than their respective crystals with Schottky defects. So the mixture promotes defect formation in both crystals. These results strengthen earlier rationalizations in catalysis, which linked the physical mixture's catalytic activity to increased oxygen vacancies in the In2O3 lattice [5][6][7]9,17]. Beyond that, physical mixtures with Schottky defects exhibited higher Frenkel defect formation percentages when compared to their non-defective counterpart. For In2O3 as the substrate, a monolayer of ZrO2 increases defects formation in both crystals, although this effect is more pronounced in the crystal which forms the monolayer.
Moreover, the ZrO2 crystal tended to form more Frenkel defects, just as in the defective crystal cases, except when it was the base crystal with edges touching the periodic boundary. It suggests that whether or not a crystal's edges touch the periodic boundary does have significant implications for the systems studied. It should also be noted that  For the ZrO2 oxygens present in bulk ZrO2 over In2O3 and for In2O3 oxygens present in In2O3 over ZrO2, another interesting pattern is the lower diffusion activation energy of the setups without intrinsic defects in comparison to those with them. This behavior also points to the great influence of lattice tensions on oxygen diffusion in the studied physical mixtures. It is argued that for crystals with intrinsic defects, there is more room for lattice tensions to be dissipated through lattice distortion, while in perfect crystals there is not an easy outlet for this accumulated energy apart from atoms jumping lattice sites.
Lastly, activation energies for ZrO2 oxygens diffusion varied little between physical mixtures, presenting values between 0.51 and 0.88 eV, which is likely due to the stronger attraction oxygen atoms have towards zirconium atoms.
The results for defect formation in physical mixtures in Table 9 do not show any obvious correlation between defect formation and diffusion activation energy. Nevertheless, some general statements can be made about them. For all cases, physical mixtures had higher defect formation than their respective crystals with Schottky defects. So the mixture promotes defect formation in both crystals. These results strengthen earlier rationalizations in catalysis, which linked the physical mixture's catalytic activity to increased oxygen vacancies in the In2O3 lattice [5][6][7]9,17]. Beyond that, physical mixtures with Schottky defects exhibited higher Frenkel defect formation percentages when compared to their non-defective counterpart. For In2O3 as the substrate, a monolayer of ZrO2 increases defects formation in both crystals, although this effect is more pronounced in the crystal which forms the monolayer.
Moreover, the ZrO2 crystal tended to form more Frenkel defects, just as in the defective crystal cases, except when it was the base crystal with edges touching the periodic boundary. It suggests that whether or not a crystal's edges touch the periodic boundary does have significant implications for the systems studied. It should also be noted that For the ZrO 2 oxygens present in bulk ZrO 2 over In 2 O 3 and for In 2 O 3 oxygens present in In 2 O 3 over ZrO 2 , another interesting pattern is the lower diffusion activation energy of the setups without intrinsic defects in comparison to those with them. This behavior also points to the great influence of lattice tensions on oxygen diffusion in the studied physical mixtures. It is argued that for crystals with intrinsic defects, there is more room for lattice tensions to be dissipated through lattice distortion, while in perfect crystals there is not an easy outlet for this accumulated energy apart from atoms jumping lattice sites.
Lastly, activation energies for ZrO 2 oxygens diffusion varied little between physical mixtures, presenting values between 0.51 and 0.88 eV, which is likely due to the stronger attraction oxygen atoms have towards zirconium atoms.
The results for defect formation in physical mixtures in Table 9 do not show any obvious correlation between defect formation and diffusion activation energy. Nevertheless, some general statements can be made about them. For all cases, physical mixtures had higher defect formation than their respective crystals with Schottky defects. So the mixture promotes defect formation in both crystals. These results strengthen earlier rationalizations in catalysis, which linked the physical mixture's catalytic activity to increased oxygen vacancies in the In 2 O 3 lattice [5][6][7]9,17]. Beyond that, physical mixtures with Schottky defects exhibited higher Frenkel defect formation percentages when compared to their non-defective counterpart. For In 2 O 3 as the substrate, a monolayer of ZrO 2 increases defects formation in both crystals, although this effect is more pronounced in the crystal which forms the monolayer. Moreover, the ZrO 2 crystal tended to form more Frenkel defects, just as in the defective crystal cases, except when it was the base crystal with edges touching the periodic boundary. It suggests that whether or not a crystal's edges touch the periodic boundary does have significant implications for the systems studied. It should also be noted that defect formation was disproportionately higher in ZrO 2 in the cases where it was not the substrate crystal compared to In 2 O 3 when it was not the substrate crystal. So, even if the stabilization effect of the periodic boundary is important, the interaction between ZrO 2 and In 2 O 3 plays a major role in promoting or limiting defect formation. More specifically, the interaction promotes more oxygen vacancies for the ZrO 2 lattice more than for the In 2 O 3 . These results are in contrast to experimental evidence [5,7]. Nonetheless, the studies that yielded this experimental evidence used a much smaller weight ratio of indium to zirconium oxide when preparing their catalysts than the weight ratio of the simulated physical mixtures.
Although the diffusion of cations was not high enough to allow calculation of its activation energies, it was high enough to be detected by a simple visual inspection in OVITO, as shown in Figure 9. So Zr atoms in In 2 O 3 have been found to stabilize it, thus limiting excessive vacancy formation on its surface [13]. Furthermore, the diffused cations would form a very diluted solution in the other crystal, with a concentration of approximately one in thousands in the receiving crystal, but this could have important implications for catalysis. Catalysts of ZrO 2 doped by Zn in similar concentrations have already been found effective in catalyzing isobutylene synthesis from ethanol [42]. defect formation was disproportionately higher in ZrO2 in the cases where it was not the substrate crystal compared to In2O3 when it was not the substrate crystal. So, even if the stabilization effect of the periodic boundary is important, the interaction between ZrO2 and In2O3 plays a major role in promoting or limiting defect formation. More specifically, the interaction promotes more oxygen vacancies for the ZrO2 lattice more than for the In2O3. These results are in contrast to experimental evidence [5,7]. Nonetheless, the studies that yielded this experimental evidence used a much smaller weight ratio of indium to zirconium oxide when preparing their catalysts than the weight ratio of the simulated physical mixtures. Although the diffusion of cations was not high enough to allow calculation of its activation energies, it was high enough to be detected by a simple visual inspection in OVITO, as shown in Figure 9. So Zr atoms in In2O3 have been found to stabilize it, thus limiting excessive vacancy formation on its surface [13]. Furthermore, the diffused cations would form a very diluted solution in the other crystal, with a concentration of approximately one in thousands in the receiving crystal, but this could have important implications for catalysis. Catalysts of ZrO2 doped by Zn in similar concentrations have already been found effective in catalyzing isobutylene synthesis from ethanol [42]. Though they do not constitute conclusive evidence, these observations about the zirconia crystal all point in the direction that it may have a more important role in the catalysis of ethanol to isobutylene than previously proposed when in a physical mixture with In2O3. Sun et al. [43] demonstrated in their theoretical DFT and experimental studies that the presence of ZrO2 played a key role in stabilizing the surface oxygen atoms of In2O3 in In2O3/ZrO2-based catalysts in the hydrogenation of CO2 to methanol. Our study provides complementary results, that ZrO2 may likewise play an important role in the reaction to obtain isobutene. Increased oxygen diffusivity and defect formation within the crystal are good indicators of its higher ability to contribute to the acetone generation via Mars Van Krevelen mechanism one important intermediate of the isobutene synthesis [10]. Besides that, the diffusion of indium atoms to the ZrO2 lattice is significant as the use of doped Though they do not constitute conclusive evidence, these observations about the zirconia crystal all point in the direction that it may have a more important role in the catalysis of ethanol to isobutylene than previously proposed when in a physical mixture with In 2 O 3 . Sun et al. [43] demonstrated in their theoretical DFT and experimental studies that the presence of ZrO 2 played a key role in stabilizing the surface oxygen atoms of In 2 O 3 in In 2 O 3 /ZrO 2 -based catalysts in the hydrogenation of CO 2 to methanol. Our study provides complementary results, that ZrO 2 may likewise play an important role in the reaction to obtain isobutene. Increased oxygen diffusivity and defect formation within the crystal are good indicators of its higher ability to contribute to the acetone generation via Mars Van Krevelen mechanism one important intermediate of the isobutene synthesis [10]. Besides that, the diffusion of indium atoms to the ZrO 2 lattice is significant as the use of doped ZrO 2 with another low valence dopant [42,44,45], Zn, has already been proved effective for the isobutylene synthesis from ethanol by acting on the rate-limiting step of the generation of acetone, intermediate of this olefin synthesis, from ethanol) [44,45]. Thus a complete understanding of this cascade reaction could involve both In 2 O 3 and ZrO 2 acting on different steps along its way. Nonetheless, the ZrO 2 -In 2 O 3 system seems to follow the moderation principle suggested by McFarland [46] which states that the oxygen must not be easily removed from the catalyst surface since it would be difficult to have it back and close the Mars Van Krevelen mechanism.

Conclusions
The results regarding oxygen diffusion pointed out a distinct pattern whereby the physical mixture of one oxide on another had the same qualitative effect as their dopingthe activation energy for oxygen diffusion was lowered in ZrO 2 and became higher in In 2 O 3 . This behavior is most likely due to the lattice tensions generated by the epitaxy between the crystals.
Further, the modified physical mixture promotes oxygen Frenkel defect formation in both crystals at higher concentrations in the case of the ZrO 2 crystal, which confirms the observations from past studies that oxygen vacancies formation is promoted in In 2 O 3 by the presence of ZrO 2 .
The present calculations also indicate the presence of very small cation diffusion between both crystals in the physical mixture, which will change the defect concentration in the indium and zirconium oxides.
For the In 2 O 3 crystal in the physical mixture, the present results matched well the vacancy formation, oxygen stabilization, and higher diffusion activation energy. In contrast, for the ZrO 2 crystal in the modified physical mixture, higher overall oxygen vacancy formation, the lower activation energy for oxygen diffusion, and the diffusion of indium cations into its structure all pointed to a role not previously considered in the literature for this catalytic system. Furthermore, the changes observed for ZrO 2 in our computational analyzes contrasted with some experimental data obtained previously. Indeed, published previous investigations show higher concentrations of ZrO 2 than those of this work. Considering that the techniques usually employed in catalyst characterization do not show lateral resolution, it was not possible to observe these ZrO 2 changes.
Thus, the physical mixture's role as a catalyst is due to the role of two modified materials and the interactions resulting from the interdiffusion of In 2 O 3 and ZrO 2 crystals.
Supplementary Materials: The following supporting information can be downloaded at: https: //www.mdpi.com/article/10.3390/ijms24032426/s1, Figure S1: -average potential energy per atom vs lattice parameter graph for In 2 O 3 described by Walsh's rigid ion model; Figure S2: elongation vs difference in temperature for In 2 O 3 described by Walsh's rigid ion model; Figure S3: MSD over time for indium atoms in the ZrO 2 (111) over In 2 O 3 (111) physical mixture with no initial defects; Figure S4: MSD over time for In 2 O 3 oxygen atoms in the ZrO 2 (211) over In 2 O 3 (111) physical mixture with no initial defects; Figure S5: ln(D) over T -1 for In 2 O 3 oxygen atoms in the ZrO 2 (211) over In 2 O 3 (111) physical mixture with no initial defects.

Conflicts of Interest:
The authors declare that they have no known competing financial interests or personal relationships that could have appeared to influence the work reported in this paper.