DFT Study of Oxygen Ion Migration in Mg-Doped Cubic Zirconia

1. Introduction

Zirconium dioxide (ZrO₂) is a promising oxide material widely used in solid oxide fuel cells, lithium-ion batteries, photocatalysis, and optoelectronics due to its high thermal and chemical stability, ionic conductivity, and mechanical strength [1,2,3,4]. ZrO₂ nanoparticles exhibit activity in the photodegradation of organic pollutants [2] and are employed in batteries to enhance anode stability and prevent degradation [3,4]. In addition, ZrO₂ is in demand in catalysis and as a dielectric material in thin-film structures [5,6,7,8], highlighting its versatility and importance among modern functional materials [9,10,11,12,13,14,15,16]. The ZrO₂ is also one of the most promising materials for nuclear engineering [17,18,19,20,21,22,23], and when doped, it is an effective phosphor for many luminescent applications [24,25,26].

ZrO₂ is a wide-bandgap semiconducting oxide with a polymorphic structure that exists under atmospheric pressure in three allotropic modifications [27,28]: monoclinic (m-ZrO₂) [29], tetragonal (t-ZrO₂), and cubic (c-ZrO₂). The monoclinic phase (baddeleyite) is thermodynamically stable at room temperature, while the transitions to the tetragonal and cubic modifications occur at 1478 K and 2650 K, respectively [30,31]. Controlling the phase stability of ZrO₂ is crucial for its application in ceramics, including solid oxide fuel cells (SOFCs), dental materials, and jewelry [32,33,34,35,36,37,38,39,40,41].

Pure zirconium dioxide is characterized by a low concentration of oxygen vacancies, which limits its ionic conductivity. The activation energy for ion transport is determined by the sum of the vacancy formation energy and the migration energy barrier. Various dopants are used to enhance conductivity and stabilize specific phases. The most common dopant elements include Yb [42], Ce [43,44], and Y [45], as well as Al, Ca, Eu, Mg, Na, Sc, and Tb [46]. For example, the addition of 10 mol% Sc stabilizes the cubic c-ZrO₂ phase, resulting in high ionic conductivity [47]. Doping with Al allows tuning of the bandgap width and enhances the absorption coefficient in the tetragonal phase of ZrO₂ [48], while lanthanide ions not only stabilize the ZrO₂ phase but also modulate its optical emission properties [49,50,51,52,53].

Doping ZrO₂ with magnesium oxide (MgO) leads to the formation of oxygen vacancies, which enhances ionic conductivity. Studies show that MgO doping affects the polymorphic transformations and mechanical properties of ZrO₂. It has been established that low concentrations of MgO stabilize the tetragonal phase, while higher concentrations promote the transition to the cubic phase. The strengthening of the material is associated with phase transitions, increased dislocation density, and the formation of intergranular MgO inclusions [54].

Previous DFT studies have shown that sulfur doping increases the formation energy of oxygen vacancies and raises the migration barrier, resulting in reduced ionic conductivity [55]. The effect of MgO on the radiation resistance of ZrO₂ was also investigated in [56], where it was demonstrated that an optimal concentration of 0.10 mol% MgO leads to the formation of a substitutional phase Zr₀.₉Mg₀.₁O₂, which reduces deformation and enhances the material’s strength. Studies on the defect structure of ZrO₂ [57] revealed that polyalanines form chains of oxygen vacancies, which may act as conductive filaments in ReRAM devices, enabling resistive switching through oxygen migration. Moreover, doping ZrO₂ with impurity ions affects the system’s total energy, contributing to phase stabilization and improved ionic conductivity, which is important for applications in solid oxide fuel cells [58].

Despite numerous studies on zirconia modeling, the question of how the spatial arrangement of Mg atoms relative to the oxygen migration path affects ionic conductivity remains unresolved. It is unclear whether substitution along the migration pathway enhances conductivity through direct interaction with vacancies or, conversely, creates obstacles due to local structural distortion.

In this work, a theoretical investigation of ionic conductivity in cubic zirconia (c-ZrO₂) stabilized with magnesium was carried out using density functional theory (DFT) with the hybrid B3LYP functional. The study focuses on the influence of oxygen vacancies and the spatial positioning of magnesium dopants on the energy profile of oxygen ion migration. The scientific novelty of this work lies in the comparative analysis of two doping strategies and in identifying the decisive role of the local coordination environment in forming effective oxygen migration pathways, which can be applied in the design of next-generation oxide electrolytes.

2. Materials and Methods

The ideal crystal of cubic zirconium dioxide (c-ZrO₂) has a face-centered cubic structure with space group symmetry Fm-3m (225). The ZrO₂ crystal structure corresponds to a fluorite-type arrangement formed in the cubic space group Fm̅3m. Zr⁴⁺ atoms are coordinated by eight O²⁻ atoms in a body-centered cubic configuration, with a Zr–O bond length of 2.20 Å. Each O²⁻ ion is bonded to four equivalent Zr⁴⁺ atoms, forming an OZr₄ tetrahedral coordination polyhedral with shared angles and edges [59].

In this study, density functional theory (DFT) was employed using the hybrid B3LYP functional implemented in the CRYSTAL23 software package [60]. A 2×2×2 supercell of cubic zirconia (c-ZrO₂) containing 96 atoms was used as the model system (see Figure 1). Quantum chemical calculations were performed on a Lenovo ThinkStation P620 workstation equipped with an AMD Ryzen™ Threadripper™ PRO 5995WX processor.

To calculate the structural and electronic properties, total energy calculations were performed without geometry optimization. The initial lattice parameter was set to a = 5.10 Å [61], and the calculated band gap was 5.2781 eV, which is comparable to the experimental value of 6.1 eV [62].

To accurately model ionic conductivity in cubic zirconia (c-ZrO₂), full geometry optimization was performed, including both atomic positions and lattice parameters. This process ensured an accurate structure necessary for further studies of ionic transport.

The system optimization was carried out using density functional theory (DFT), with the following convergence criteria: an energy threshold of 10⁻⁶ eV, a gradient threshold of 4.5·10⁻⁴ eV/Å, and a maximum atomic displacement of 1.8·10⁻³ Å. These settings enabled high precision in calculating the structural and energetic characteristics of the system.

Oxygen vacancies were introduced during the calculations to simulate defects that play a key role in ionic conductivity. These vacancies were created by removing oxygen atoms from the lattice, thereby forming possible migration paths for oxygen ions through the material.

In addition, magnesium dopants were introduced into the system by substituting some of the zirconium atoms, allowing the investigation of magnesium’s influence on ionic conductivity. Magnesium in the structure can shift the positions of vacancies and create new defect pathways, affecting oxygen ion migration and, consequently, the material’s conductivity.

These calculations involving magnesium doping and oxygen vacancies provide a basis for further analysis of ionic conductivity in c-ZrO₂, which is important for developing advanced materials for solid oxide fuel cells and other applications. The analysis of the system’s total energy as a function of oxygen vacancy concentration and magnesium content allowed us to examine the influence of defects on the phase stability and electronic properties of zirconia.

3. Results

3.1. Investigation of the Electronic Properties of Magnesium-Stabilized Cubic Zirconia

It is well known that the creation of oxygen vacancies in ZrO₂ contributes to structural stabilization and enhances ionic conductivity. These vacancies are modeled from first principles by removing one or more oxygen atoms from the supercell.

The modeling is carried out using standard density functional theory (DFT) methods. Total and partial densities of electronic states (TDOS and PDOS) are calculated using the hybrid B3LYP functional. The valence electron configurations are 4d²5s² for zirconium and 2s²2p⁴ for oxygen. Oxygen vacancies are modeled within a cubic ZrO₂ supercell containing 96 atoms. The geometric parameters of the nearest atoms are not optimized.

As shown in Figure 2, in pristine ZrO₂, the valence band is primarily composed of oxygen 2p orbitals, with a maximum around −4.5 eV, while the conduction band starts above the Fermi level and includes contributions from Zr 4d states. Upon the creation of an oxygen vacancy, an additional narrow localized state appears near the bottom of the conduction band, associated with the remaining unpaired electron.

In the PDOS spectra (Figure 2), the formation of a single oxygen vacancy in ZrO₂ leads to the appearance of a localized state approximately 0.1 eV below the Fermi level. At the same time, the vacancy contributes negligibly to the valence band. The main contributions to the valence and conduction bands originate from oxygen (O 2p) and zirconium (Zr 4d) atoms, respectively. These findings align with the general model of defect state formation in oxide dielectrics and demonstrate the potential role of vacancies in enabling ionic conductivity.

With an increase in the concentration of oxygen vacancies in ZrO₂ to three per supercell, the total density of states retains the typical oxide profile with a wide band gap (Figure 3). However, in the PDOS projected onto the vacancy sites, a localized state appears within the band gap (approximately 0.1–0.5 eV below the Fermi level), associated with the defects.

The height and width of the peaks indicate an accumulation of defect states that maintain a discrete character, suggesting that the vacancies are spatially separated. This confirms their influence on ionic conductivity.

Figure 4 presents the total and partial density of states (DOS and PDOS) for a zirconia structure containing one magnesium substitutional atom and one oxygen vacancy. The upper part of the valence band is primarily formed by oxygen 2p electrons, while the conduction band is mainly composed of zirconium 4d electrons. As shown in the PDOS, magnesium contributes minimally to the electronic structure of the system, which is consistent with its Mg²⁺ electronic configuration and the absence of unoccupied d-states.

It is clearly observed that the defect states associated with the oxygen vacancy are in the lower part of the conduction band (approximately 1.3 eV above the Fermi level). The position of the oxygen vacancy state has shifted compared to the undoped system.

Figure 5 shows the total and partial density of states (DOS and PDOS) for a zirconia system containing two oxygen vacancies and two magnesium atoms. As in previous cases, the valence band is mainly formed by oxygen 2p orbitals, while the conduction band consists primarily of zirconium 4d states. The presence of oxygen vacancies leads to the formation of defect levels in the upper part of the band gap, which in this system are located closer to the conduction band edge—approximately 4.5–5.5 eV above the Fermi level.

These defect states exhibit greater intensity than in the case with a single vacancy, indicating a cumulative effect of electron localization as the number of vacancies increases. The contribution of magnesium atoms near the Fermi level remains minimal, confirming their role as compensating dopants that influence charge neutrality but do not directly participate in the formation of the band structure. Thus, double magnesium substitution stabilizes the vacancies, shifting their associated states into the conduction band region.

Figure 6 presents the total and partial density of electronic states (DOS and PDOS) for a zirconia system containing three oxygen vacancies and three magnesium atoms. The total DOS shows a characteristic separation between the valence and conduction bands with a wide band gap. As seen in the PDOS, the contribution from oxygen vacancies forms a clearly defined narrow state located near the bottom edge of the conduction band, approximately 5.5–6.0 eV above the Fermi level.

This indicates that as the concentration of vacancies increases, the energy level of the associated states shifts upward while maintaining a localized character. Additional magnesium substitution compensates for the excess charge, but the contribution from Mg remains minimal in both the valence and conduction bands. Thus, an increase in the number of defects enhances localization effects and may promote the formation of new channels for ionic conduction.

An analysis of the total and partial density of states (DOS and PDOS) was conducted for ZrO₂ systems with one and three oxygen vacancies and varying numbers of Mg atoms. It was shown that the electronic structure is strongly dependent on the defect configuration.

In the system with a single vacancy, a localized state appears below the Fermi level, which may hinder ionic conductivity. The addition of one Mg atom shifts this state above the Fermi level and reduces its intensity, indicating charge compensation and state delocalization, thereby promoting improved oxygen mobility.

With an increase in the number of Mg atoms and vacancies to 2 and 3, further suppression of localized states is observed, while the contribution of Mg to the valence and conduction bands remains minimal, confirming its role as a charge compensator.

In the case of three vacancies without magnesium, pronounced unoccupied states are detected above the Fermi level—an indication of electronic traps. Their disappearance upon Mg introduction highlights the stabilizing effect of substitution.

Thus, configurations with a Mg:V_O ratio close to 1:1 are the most effective for enhancing ionic conductivity by eliminating localized states and creating favorable pathways for oxygen migration. The results are consistent with literature data on perovskites and emphasize the importance of targeted defect engineering in designing ZrO₂-based ion-conducting materials.

3.2. Investigation of Ionic Conductivity

Defects significantly influence the chemical and physical properties of metal oxides. Point and bulk defects often determine their catalytic, electrical, optical, and mechanical activity. Surface oxygen vacancies play a key role in heterogeneous catalysis. The oxygen vacancy formation energy on the surface (Eₒᵥₐc) serves as an important parameter for predicting the catalytic properties of oxide materials.

To model an oxygen vacancy, one oxygen atom was removed from the supercell. The vacancy formation energy was calculated using the following expression:

$$E^F=E^d+E_O-E^p$$

(1)

where E^F- is the oxygen vacancy formation energy, E^d- is the total energy of the defective system, E_O- is the energy of an isolated oxygen atom, and E^p - is the total energy of the pristine (defect-free) system.

All calculations in this study were performed using the B3LYP hybrid functional within the framework of density functional theory (DFT), as implemented in the CRYSTAL23 software package.

The oxygen vacancy formation energy on the surface (Eₒᵥₐc) is an important parameter that determines the catalytic activity of metal oxides. Estimating these energies can support the development of promising catalysts using data-driven approaches. In the present study, Eₒᵥₐc values were determined for c-ZrO₂. It is known that Eₒᵥₐc is strongly influenced by the band gap width, lattice formation energy, and electron affinity.

The oxygen vacancy formation energy was calculated as the difference between the total energies of the defective supercell and the ideal crystal, considering the chemical potential of an oxygen atom, defined as the energy of a single O atom in its ground (doublet) state.

The calculated oxygen vacancy formation energy is 12.6 eV, which is consistent with the typical overestimation of bond strength inherent to hybrid functionals. Since all subsequent calculations (including substitutional defects and defect complexes) were performed using the B3LYP functional, the use of this value ensures methodological consistency and allows for a reliable comparison of the relative energetic characteristics of different configurations.

3.2.1. Ionic Conductivity via Oxygen Vacancy-Mediated Oxygen Ion Migration

In this study, we examined the mechanism of oxygen migration in the defective structure of c-ZrO₂. A 96-atom supercell was constructed, in which one oxygen atom was initially removed to create a vacancy. Subsequently, an oxygen atom from the upper layer was gradually displaced in 0.2 Å steps until it filled the oxygen vacancy at approximately 2.5 Å. This process was repeated for other layers of the structure (Figure 7, curve 1). At each step, the total energy of the system was calculated.

In the graph presented here and in subsequent analyses, the vertical axis represents the relative total energy of the system. It was found that the system’s energy reaches a minimum when the oxygen atom occupies the vacancy and increases when it is near zirconium atoms. The dependence of the total energy on the oxygen atom’s position exhibits a sinusoidal pattern, indicating the presence of periodic energy barriers and stable configurations within the structure.

Figure 7 shows the energy profile along the oxygen ion migration path in defective cubic zirconia (c-ZrO₂) through vacancy sites. Curve 1 in Figure 7 corresponds to the migration of an oxygen ion into a single vacancy site (Figure 8a), while curve 2 represents migration through three oxygen vacancies (Figure 8b). The energy maxima observed at distances of 1.2–1.4 Å, 3.7–3.9 Å, and 6.3 Å indicate energy barriers that arise when the oxygen atom passes through a zirconium plane, interacting with Zr atoms and increasing the system’s energy.

The energy minimum, observed at 2.5 Å, 5 Å, and 7.5 Å, corresponds to positions where the oxygen atom sequentially occupies vacancy sites, leading to system stabilization and energy minimization. A comparison of the migration barriers for an oxygen ion in a structure with one vacancy (2.7 eV) and along a path with three vacancies (1.32 eV) reveals a reduction of 1.4 eV.

It should be noted that aligning multiple oxygen vacancies along a single trajectory is statistically unlikely. Nevertheless, the results confirm that the formation of oxygen vacancies significantly lowers the energy barriers for oxygen migration, thereby enhancing the material’s ionic conductivity.

3.2.2. Ionic Conductivity of ZrO₂ Stabilized with MgO Aligned Along the Migration Path

To assess the effect of doping on the stability and ionic conductivity of zirconia, the energetic characteristics of systems containing one, two, and three magnesium atoms substituting zirconium atoms near a pre-existing oxygen vacancy were calculated. The total energy of the system consistently decreased with an increasing number of magnesium atoms: −6391.30 a.u. for one Mg, −6469.26 a.u. for two, and −6547.24 a.u. for three, indicating an increase in the thermodynamic stability of complex defects Zr₃₂₋ₓMgₓO₆₄₋ₓ with higher dopant concentrations. These values are significantly lower than the energy of the system with a single vacancy (−6237.92 a.u.) and the undoped system (−6313.33 a.u.), confirming the energetic favorability of forming complex defects.

The substituting Mg²⁺ atoms aligned along the oxygen migration path lie in the same (100) crystallographic plane. This arrangement contributes to stabilizing the system containing oxygen vacancies.

In Figure 9a, the profile of the relative total energy during oxygen migration in the Zr₃₁Mg₁O₆₃ system (with one Mg atom aligned along the migration path) shows a shift in the positions of energy minima and maxima compared to the undoped structure. The graph illustrates the energy variation as a function of the oxygen ion’s displacement coordinate.

At the points corresponding to vacancies (2.5 Å, 5.0 Å, and 7.5 Å), where deep minima were ideally expected, local maxima are instead observed: 0.462 eV, 2.264 eV, and 2.424 eV, respectively. This indicates that the Mg²⁺ ion located near the trajectory disrupts the symmetry of the potential curve and distorts the local electrostatic field. In the initial segment (0–2 Å), the energy gradually decreases, reaching a minimum of 0.067 eV at 2.0 Å, but then instead of further decreasing, it rises. At 3.7–3.9 Å and 6.3–6.4 Å, where the oxygen ion passes near zirconium, the energy increases to 1.26 and 2.37 eV, respectively. Similarly, when passing through the Mg-Zr planes (1.2–1.4 Å), the energy is approximately ~0.35 eV, which corresponds to the energy maximum in the undoped system with oxygen vacancies.

From the dependence of the system’s total energy on the distance with one magnesium atom, it is seen that the total energy is lower compared to pure ZrO₂, indicating the thermodynamic stability of the defect structure. In the migration within the initial structure with a single oxygen vacancy, the maximum barriers reach ~2.7 eV. Thus, magnesium doping leads to a decrease in the total energy of the system and a reduction in migration barriers, although it is accompanied by a shift in critical points and a disruption of the migration path symmetry.

In Figure 9b, the dependence of the total energy of the system on the oxygen ion migration distance in the Zr₃₀Mg₂O₆₂ structure, where two zirconium atoms are replaced by magnesium along the migration path, is shown.

In the interval from 0.2 to 2.5 Å, a monotonic decrease in energy from 1.98 to 0.27 eV is observed with a local minimum of 0.048 eV at 2.2 Å and an absolute minimum (0 eV) at 2.9 Å. The minima at 2.2–2.9 Å are explained by the stable localization of oxygen near the vacancy.

Values at 1.4 Å and 3.9 Å, corresponding to the passage of oxygen through the Mg–Zr and Zr–Zr coordination planes, are characterized by local energy bends: at 1.4 Å, the energy is 0.774 eV, and at 3.9 Å — 1.18 eV. This indicates the interaction of oxygen with nearby magnesium and zirconium atoms, creating local distortions in the electrostatic field of the lattice.

At 5.0 Å (the second vacancy), there is a sharp energy rise to 3.03 eV, indicating an increase in the potential energy of the electrostatic field. The oxygen ion exits the potential well corresponding to the vacancy and crosses the zirconium plane (at 5.4 Å — energy 2.39 eV). The maximum barrier of 4.62 eV is observed at 7.5 Å, which is associated with the completion of the migration path.

In Figure 10c, the dependence of the reduced total energy of the Zr₂₉Mg₃O₆₁ system corresponding to the oxygen ion migration along a trajectory with three magnesium ions is shown. In the initial region (0.2–2.5 Å), where the ion passes through the magnesium–zirconium planes, the energy gradually decreases from 2.27 eV to 0.96 eV. From 2.7 Å to 4.7 Å, a deep potential well is observed with energy ranging from 0.01664 eV to 0.03447 eV, which corresponds to the ion passing near Mg ions. However, when exiting this region (r > 5 Å), as the oxygen begins to leave the Mg influence zone and approaches the zirconium planes, the energy sharply increases. At 5.2 Å, the energy is 1.32 eV and continues to rise to 3.31 eV at 7.5 Å. Despite the relatively high final barrier, the system with three magnesium atoms exhibits the lowest total energy among all studied configurations, indicating its high stability. This emphasizes that with increasing Mg concentration, the system becomes more stable, but the electrostatic field becomes less symmetric, and the barriers may increase again in certain regions.

The analysis revealed that local substitution of zirconium atoms with magnesium near an oxygen vacancy has a significant impact on both the thermodynamic stability of the structure and the energy profile of oxygen ion migration. As the number of Mg atoms increases, the total energy of the system decreases consistently, indicating the formation of more stable configurations: −6391.30 a.u. (1 Mg), −6469.26 a.u. (2 Mg), −6547.24 a.u. (3 Mg).

However, the migration energy profiles exhibit a complex dependence on the number and arrangement of the dopant atoms. With one Mg atom, the energy barriers shift but remain moderate compared to the undoped system (2.7 eV). In the case of two Mg atoms, deep potential wells are formed, but the end of the migration path shows a sharp energy increase up to 4.62 eV, indicating hindered migration. The third configuration, with three Mg atoms, offers the best stabilization of oxygen near the vacancy (energies ~0 eV), but is also accompanied by a rise in the energy barrier to 3.3 eV at the later stage of migration.

Thus, optimal ionic conductivity is achieved by balancing the number and spatial arrangement of dopant atoms: excessive magnesium can enhance local stabilization but may disrupt the uniformity of the migration channel, increasing energy barriers at the final stages of migration.

Table 1. Energy barriers for oxygen migration in c-ZrO₂ with different numbers and arrangements of Mg atoms.

Doping configuration	Number of Mg	Min. barrier, eV	Max. barrier, eV
Zr₃₁Mg₁O₆₃	1	0.4	2.4
Zr₃₀Mg₂O₆₂	2	0.5	3.1
Zr₂₉Mg₃O₆₁	3	0.8	1.32

Table 1. Energy barriers for oxygen migration in c-ZrO₂ with different numbers and arrangements of Mg atoms.

Doping configuration	Number of Mg	Min. barrier, eV	Max. barrier, eV
Zr₃₁Mg₁O₆₃	1	0.4	2.4
Zr₃₀Mg₂O₆₂	2	0.5	3.1
Zr₂₉Mg₃O₆₁	3	0.8	1.32

3.2.3. Ionic Conductivity of ZrO₂ Stabilized with MgO with Spatially Distributed Mg Ions

As a continuation of the investigation into the effect of magnesium doping on ionic conductivity in cubic zirconia (c-ZrO₂), a series of models were examined in which Mg atoms substituted zirconium atoms outside the oxygen ion migration trajectory. In the structures with Mg doping positioned away from the migration path, the magnesium atoms replace zirconium in various planes, primarily oriented along (001) and (100). Two Mg²⁺ ions are localized within the same (001) plane at z = 2.5645 Å, while the first Mg atom is displaced to the plane at z = 7.6936 Å, corresponding to another parallel (001) plane. This arrangement of magnesium ions, removed from the migration trajectory, eliminates their direct interaction with oxygen vacancies and limits their effect on the total system energy to background electrostatic influences (Figure 11).

Three configurations were selected for calculation: Zr₃₁Mg₁O₆₃, Zr₃₀Mg₂O₆₂, and Zr₂₉Mg₃O₆₁, in which one, two, and three Mg atoms, respectively, substitute zirconium atoms outside the oxygen ion migration path, thus avoiding direct interaction with vacancies. The total energies of these systems were −6391.2362 Hartree, −6469.1623 Hartree, and −6547.2159 Hartree, respectively. This indicates that increasing MgO concentration contributes to thermodynamic stabilization.

From Figure 12a, corresponding to the Zr₃₁Mg₁O₆₃ configuration, the overall shape of the energy profile resembles the energy variation during oxygen ion migration via vacancies. The minimum energy is observed at 2.85 Å and equals 0 eV, indicating the position of the oxygen vacancy. At other points, such as 1.2–1.4 Å, the energy reaches 0.92 and 0.90 eV, respectively, while near 3.65–3.85 Å, maxima of 0.99–1.01 eV are observed, despite the absence of direct interaction from the Mg atom. Further along, at 5.85–6.25 Å, the energy rises again to 0.94–1.12 eV.

This behavior indicates that the Mg atom, located outside the migration path, does not exert a significant stabilizing effect on the migration trajectory but does alter the electrostatic field. In contrast to the case where Mg ions are aligned along the migration path, the graph displays a sinusoidal pattern, and the peak energy values do not exceed 1.12 eV.

Figure 12b shows the dependence of the reduced total energy of the system corresponding to the Zr₃₀Mg₂O₆₂ configuration, in which two magnesium atoms substitute zirconium atoms outside the migration path. At distances corresponding to the position of the vacancy, small potential barriers are observed. The minimum energy is reached at 2.9 Å and equals 0 eV, while the maximum occurs at distances of 1.2 Å (1.3 eV), 3.8 Å (0.85 eV), 5.0 Å (1.1 eV), 6.4 Å (1.3 eV), and 7.5 Å (1.1 eV). The placement of Mg²⁺ outside the migration path reduces the height of the energy barrier compared to when Mg²⁺ is located along the path.

Figure 12c shows the dependence of the reduced total energy of the system corresponding to the Zr₂₉Mg₃O₆₁ configuration. Here, three magnesium atoms substitute for zirconium outside the oxygen migration path. The energy along the trajectory fluctuates within 0.0–1.05 eV, with an absolute minimum at 2.9 Å (0 eV), corresponding to the position of the oxygen vacancy. Additional energy minima are observed at 0.2 Å (0.05 eV), 2.2 Å (0.13 eV), and 5.4–5.6 Å (~0.26 eV), indicating stabilization of the oxygen ion near these points. Maxima are observed at 1.2 Å (1.01 eV), 2.5 Å (0.5 eV), 3.7 Å (0.9 eV), 5.0 Å (1.05 eV), 6.4 Å (0.85 eV), and 7.5 Å (1.05 eV), indicating moderate energy barriers. Thus, introducing three Mg atoms outside the migration path preserves the symmetry of the electrostatic field and lowers the energy amplitude, providing a stable oxygen migration pathway.

The analysis of the influence of magnesium substitutions on the total energy of the oxygen migration system in the ZrO₂ structure has shown that the nature and effectiveness of the modification directly depend on both the number of Mg²⁺ ions and their spatial arrangement relative to the migration path.

Table 2. Energy Barriers for Oxygen Migration in c-ZrO₂ with Different Numbers and Positions of Mg Atoms.

Doping configuration	Number of Mg	Min. barrier, eV	Max. barrier, eV
Zr₃₁Mg₁O₆₃	1	0.4	1.12
Zr₃₀Mg₂O₆₂	2	0.5	1.34
Zr₂₉Mg₃O₆₁	3	0.5	1.05

Table 2. Energy Barriers for Oxygen Migration in c-ZrO₂ with Different Numbers and Positions of Mg Atoms.

Doping configuration	Number of Mg	Min. barrier, eV	Max. barrier, eV
Zr₃₁Mg₁O₆₃	1	0.4	1.12
Zr₃₀Mg₂O₆₂	2	0.5	1.34
Zr₂₉Mg₃O₆₁	3	0.5	1.05

When magnesium ions are located outside the migration trajectory, the variation of the normalized total energy becomes smooth, and the maximum energy barriers are 0.9 eV and 1.3 eV, respectively. The results of the performed calculations show that the spatial positioning of the implanted magnesium ions outside the migration path significantly reduces the height of the potential barriers. Overall, increasing the number of Mg²⁺ ions, when optimally positioned either along or outside the trajectory, contributes to a decrease in the energy barriers.

4. Discussion

The creation of a single oxygen vacancy in cubic ZrO₂ leads to the formation of defect levels within the band gap and a reduction in the energy barrier for oxygen ion migration. A particularly significant factor is the presence of lattice oxygen atoms along the migration pathway—when these are absent, the barrier is considerably reduced, facilitating ionic conductivity. These findings are consistent with theoretical studies that emphasize the critical role of vacancies in zirconia-based oxide systems [63,64].

When zirconium atoms are substituted with magnesium along the oxygen migration pathway, stable complex defects are formed, as evidenced by the consistent decrease in the system’s total energy with increasing Mg content. However, the migration energy profiles exhibit a nonlinear dependence on the number and position of dopant atoms along the channel. For instance, with one Mg atom, a moderate reduction of the barrier to 1.27 eV is observed; with two Mg atoms, a sharp increase in the primary barrier to 3.03 eV occurs midway through the path; and with three Mg atoms, the most balanced behavior is achieved—a continuous channel is formed with nearly zero energy and moderate barriers of 0.96–1.32 eV, providing optimal conditions for migration. These differences are explained by local lattice distortions and charge redistribution induced by different doping configurations along the channel [65].

In contrast, substituting zirconium with magnesium outside the oxygen migration pathway also leads to a decrease in the system’s total energy with increasing Mg content but does not cause significant changes in the migration energy landscape. With two and three substitutions outside the migration channel, the barriers remain moderate (with a maximum of 0.96–1.3 eV), the profile maintains symmetry, and no sharp local distortions are observed. This indicates that such doping exerts only a background electrostatic effect and smooths the potential landscape. Similar conclusions are presented in studies on the role of vacancies and their correlation with doping in oxide systems [64].

A comparative analysis of the oxygen migration energy profiles for zirconium substitution with magnesium along and outside the migration pathway reveals different impacts on ionic conductivity. Mg atoms positioned outside the pathway reduce the overall system energy (Table 3) but do not significantly lower the migration barriers. In contrast, Mg substitution along with the migration trajectory can either enhance or hinder conductivity, depending on the number and distribution of the substitutions. Considering the experimental impossibility of aligning Mg substitutions strictly along a specific migration channel direction, the most effective approach is to introduce the maximum possible amount of Mg substitutions throughout the crystal volume in combination with oxygen vacancies.

This strategy enables the formation of a continuous, energetically favorable path for oxygen migration. These conclusions align with computational results on perovskites, where the importance of coordinated defect structure tuning for achieving high ionic conductivity is emphasized [66], as well as with modeling studies demonstrating the impact of rhodium and silver doping on oxygen behavior in oxide and perovskite systems [67,68,69], and the significance of the local electronic environment near vacancies and color centers in fluoride [70,71] and oxides [72,73,74].

5. Conclusions

The obtained results are of great importance for understanding the mechanisms of crystalline structure stabilization in ZrO₂ and can be applied in the development of materials for solid oxide fuel cells (SOFCs) and oxygen sensors.

Thus, the main objective of this study has been achieved. A fundamental relationship has been established between the number and spatial distribution of magnesium substitutions and the shape of the energy landscape for oxygen migration. The findings confirm that controlled defect distribution (vacancies and substitutions) is a key factor in optimizing ionic conductivity in zirconia, opening new prospects for the targeted design of solid oxide electrolytes and sensing materials.