First-principles prediction of high oxygen-ion conductivity in trilanthanide gallates Ln₃GaO₆

ABSTRACT

We systematically investigated trilanthanide gallates (Ln₃GaO₆) with the space group Cmc2₁ as oxygen-ion conductors using first-principles calculations. Six Ln₃GaO₆ (Ln = Nd, Gd, Tb, Ho, Dy, or Er) are both energetically and dynamically stable among 15 Ln₃GaO₆ compounds, which is consistent with previous experimental studies reporting successful syntheses of single phases. La₃GaO₆ and Lu₃GaO₆ may be metastable despite a slightly higher energy than those of competing reference states, as phonon calculations predict them to be dynamically stable. The formation and the migration barrier energies of an oxygen vacancy (V_O) suggest that eight Ln₃GaO₆ (Ln = La, Nd, Gd, Tb, Ho, Dy, Er, or Lu) can act as oxygen-ion conductors based on V_O. Ga plays a role of decreasing the distances between the oxygen sites of Ln₃GaO₆ compared with those of Ln₂O₃ so that a V_O migrates easier with a reduced migration barrier energy. Larger oxygen-ion diffusivities and lower migration barrier energies of V_O for the eight Ln₃GaO₆ are obtained for smaller atomic numbers of Ln having larger radii of Ln³⁺. Their oxygen-ion conductivities at 1000 K are predicted to have a similar order of magnitude to that of yttria-stabilized zirconia.

Graphical abstract

KEYWORDS:

• Oxygen-ion conductor
• Ln₃GaO₆
• lanthanide gallate
• first-principles material design

CLASSIFICATION:

• 50 Energy Materials
• 207 Fuel cells / Batteries / Super capacitors
• 401 1st principle calculations
• 107 Glass and ceramic materials
• 404 Materials informatics / Genomics

1. Introduction

Oxygen-ion conductors with high oxygen-ion conductivities (σ_O) have been developed for important applications [1,2], such as electrolytes of solid oxide fuel cells (SOFC), oxygen separation membranes, and gas sensors. Currently, yttria-stabilized zirconia (YSZ) is widely used because of its advantages, such as abundance, chemical stability, non-toxicity, and low cost. This material is known to exhibit a σ_O of ~10⁻² S/cm at a high temperature of 1000 K [2,3]. For useful industrial applications, it is necessary to be further improved to have similar σ_O values at lower temperatures or a higher σ_O at a similar temperature. Indeed, there are several oxides such as Gd-doped CeO₂ (GDC) [4] and pure or Er-doped δ-phase of Bi₂O₃ [5,6], which have been reported to have higher σ_O at the same temperature than YSZ. However, there have been great efforts for the development of new oxygen-ion conductors that can substitute for YSZ with enough merits in practice.

Lanthanide gallates (Ln–Ga–O) have also been investigated as oxygen ionic conductors. Actually, Sr- and Mg-doped LaGaO₃ (LSGM) [7–11] with a perovskite structure have been widely developed as an electrolyte for intermediate temperature (~500–750 °C) [12] SOFCs because their σ_O is higher than that of YSZ at this temperature. In addition, trilanthanide gallates in another composition ratio (Ln₃GaO₆) have also been introduced as oxygen-ion conductors. Purohit et al. [13] reported that pure and Sr- or Ca-doped Nd₃GaO₆ achieved a σ_O of 2.7 × 10⁻⁴ and 0.7 × 10⁻² S/cm at 800 °C, respectively. Iakovleva et al. [14] reported that a Sr-doped Gd₃GaO₆ achieved a total conductivity of 0.9 × 10⁻² S/cm (with σ_O of 0.3 × 10⁻² S/cm and the remaining conductivity likely originating from protons) at 800 °C. Their σ_O values are not significantly different from that of YSZ at ~1000 K. Some Ln₃GaO₆ (Ln = Nd, Sm, Eu, Gd, Tb, Dy, Ho, or Er) were synthesized as a single phase with the space group Cmc2₁ [15–17]. The abovementioned experimental reports suggest that other Ln₃GaO₆ compounds may show similar (or higher) σ_O values.

In contrast to experimental reports on Ln₃GaO₆, first-principles calculations on these ternary oxides are very rare. However, for binary lanthanide sesquioxides (Ln₂O₃) [18–21], the locations of the f-levels inside the O 2p–Ln 5d band gap (E_g) have been predicted using first-principles calculations (see Section 3.2). Whether the localized f-levels exist inside or outside of E_g(O 2p–Ln 5d) depending on the type of Ln (number of f-electrons) changes the characteristics of the highest occupied molecular orbital (HOMO) and the lowest unoccupied molecular orbital (LUMO), and the E_g value. Therefore, it is useful to investigate the electronic structures of Ln₃GaO₆ by the first-principles calculations.

In this study, we used the systematic first-principles calculations to determine the structural and energetic properties based on the crystal and electronic structures of Ln₃GaO₆ with all types of Ln. We classified the type of Ln₃GaO₆ according to the location of f-levels in the band gap. Then, for the Ln₃GaO₆ type, whose electronic structures are less affected by the f-levels, we further investigated the energetic and dynamical stabilities, and the properties related to σ_O, namely, the formation energy (E_v) and migration barrier energy (E_m) of the oxygen vacancy (V_O). Finally, we compared the computed oxygen-ion diffusivity (D_O) among the Ln₃GaO₆ types including the V_O by the first-principles molecular dynamics (FPMD) and predicted their σ_O values.

2. Methodology

2.1. Crystal structure

Figure 1 shows the crystal structure of Ln₃GaO₆ in the orthorhombic crystal system with the space group Cmc2₁. This crystal structure consists of two, one, and four types of Ln, Ga, and O sites, respectively. The Wyckoff letters with the multiplicity of La(I), La(II), Ga, O(I), O(II), O(III), and O(IV) are 8a, 16b, 8a, 16b, 16b, 8a, and 8a, respectively, when we consider them in the crystal structure shown in Figure 1. The coordination numbers (CNs) of Ln [both of La(I) and La(II)] and Ga with respect to O are seven and four, respectively. Note that the CN of Ln in the binary Ln₂O₃ (both of cubic and hexagonal crystal systems) is six. O(I), O(II), O(III), and O(IV) have three, four, four, and three chemical bonds with Ln, respectively. In addition, O(I), O(III), and O(IV) have one chemical bond with Ga, whereas O(II) has no chemical bond with Ga, but only with Ln. In total, we obtain 15 types of Ln₃GaO₆ (Ln = La, Ce, Pr, Nd, Pm, Sm, Eu, Gd, Tb, Dy, Ho, Er, Tm, Yb, or Lu) in this crystal structure.

Figure 1. Crystal structure of Ln₃GaO₆ with space group Cmc2₁. The supercell (1 × 1 × 2 conventional orthorhombic unit cell) with 24 Ln, 8 Ga, and 48 O atoms is shown. Different types of Ln and O sites are presented by different colors.

2.2. First-principles calculations

All first-principles calculations were performed using the project augmented wave (PAW) [22,23] method implemented in the Vienna Ab-initio Simulation Package (VASP) [24,25] within the framework of a generalized gradient approximation (GGA) of Perdew–Burke–Ernzerhof (PBE) form [26]. As the valence electrons, 5p, 5d, and 6s for Ln (5s is also included for La, Ce, Pr, Nd, Pm, and Sm), 3d, 4s, and 4p for Ga, and 2s and 2p for O were considered. In addition, on-site Coulomb interaction [27] with an effective U−J of 5 eV (GGA+U) for the f-orbitals of Ln was partially employed in the calculations. The value of U−J was selected by referring to previous theoretical calculations [19,28–30] on Ln₂O₃, which suggested an adequate range of the value as 5–7 eV with crosschecks on experimental reports. For the GGA+U method, 4f, 5s, 5p, 5d, and 6s electrons of Ln were also considered as the valence electrons. Henceforth, we refer to the former and latter methods as the GGA/w.o.f (without f-valences but in the frozen core) and GGA+U/w.f methods (with f-valences), respectively.

The structural relaxations of primitive unit cells with the associated changes in lattice constants and atomic coordinates were performed until the interatomic force on each atom was reduced to within 0.005 eV/Å. The cutoff energy was set as 500 eV. Brillouin zone was sampled by using Γ-centered 4 × 4 × 4 meshes. Gaussian smearing for the Brillouin zone integrations was employed with widths of 0.1 eV. Electronic spin polarizations were basically turned on except for the case where the up spin and down spin were confirmed to be compensated by each other.

We also computed E_v and E_m because their low values are favorable for achieving a higher σ_O [31,32]. The E_v values were computed using the supercell method [33] as follows:

(1) $E_v^q=E\left(\mathrm{L}{\mathrm{n}}_3\mathrm{G}\mathrm{a}{\mathrm{O}}_6:V_{\mathrm{O}}^q\right)-E\left(\mathrm{L}{\mathrm{n}}_3\mathrm{G}\mathrm{a}{\mathrm{O}}_6\right)+{\mu }_{\mathrm{O}}+q\left(E_{\mathrm{H}\mathrm{O}\mathrm{M}\mathrm{O}}+E_{Fermi}\right)+\mathrm{\Delta }{E}_{\mathrm{L}\mathrm{Z}}$(1)

where E(Ln₃GaO₆:V_O^q) is the energy of a supercell including a V_O, E(Ln₃GaO₆) is the energy of a supercell of pure Ln₃GaO₆, μ_O is the chemical potential of O, q is the charge of V_O, E_HOMO is the eigenvalue of the HOMO level mainly formed by O 2p, E_Fermi is the Fermi level as a variable, and ΔE_LZ is the correction term suggested by Lany and Zunger [34,35] to compensate for the image charge interactions between supercells for the charged V_O. The calculations for the supercell with a V_O were performed until the interatomic forces on each atom were reduced to within 0.02 eV/Å for fixed lattice constants. The sizes of the supercells (80 atoms) were set as 1 × 1 × 2 orthorhombic conventional unit cells as shown in Figure 1. The Γ-centered 2 × 2 × 2 meshes were used for the k-space sampling. The transformations between primitive and orthorhombic conventional unit cells were performed by SPGLIB implemented in PHONOPY [36,37]. The dielectric constants for the correction term of the image charge interaction were obtained using primitive cells with the density functional perturbation theory (DFPT) [38]. When V_O are charged, neutralizing background charges were added.

The phonon calculations of Ln₃GaO₆ were performed by using PHONOPY [36,37]. To calculate the phonon frequencies, we combined the supercell approach and DFPT [38]. Force constants of the oxides were calculated by applying DFPT at the Γ point to the supercell model (the same size supercell with 80 atoms) of the oxides. Phonon frequencies were then calculated from the force constants. Ln₃GaO₆ compounds were judged to be dynamically stable if they did not have imaginary phonon frequencies in their phonon densities of states (DOSs).

The E_m values were calculated using the climbing image nudged elastic band (CI-NEB) method [39,40] with three intermediate images. The CI-NEB calculations were performed until the forces decreased below 0.03 eV/Å with a spring constant of 5 eV/Ų between the images. We employed a doubly charged oxygen vacancy (V_O²⁺) for E_m assuming that it was formed by the reduction in cation valence upon doping.

2.3. First-principles molecular dynamics

In order to investigate D_O, we performed FPMD using the same size supercell as for the computations of E_v and E_m including one V_O²⁺. For computational efficiency, the cutoff energy was decreased to 300 eV, only Γ-point was sampled for the Brillouin zone, and the valence electrons were considered as 5p, 5d, and 6s for Ln (5s is also included for La and Nd), 4s and 4p for Ga, and 2s and 2p for O. For O, a soft pseudopotential with a reduced maximum cutoff energy (288 eV) was employed. The GGA/w.o.f method was used without spin polarizations.

The D value of each chemical element at high temperatures can be obtained from the mean square displacement (MSD) after a long simulation time t by the following equation [41,42]:

(2) $D\left(t\right)={lim}_{t\to \mathrm{\infty }}\frac{\mathrm{M}\mathrm{S}\mathrm{D}}{6t}$(2)

Here, MSD(t) can be obtained by

(3) $\mathrm{M}\mathrm{S}\mathrm{D}\left(t\right)=\frac{1}{n}\sum _i{\left|r_i\left(t\right)-r_i\left(0\right)\right|}^2$(3)

where n is the number of atoms of each chemical element, and r_i(0) and r_i(t) are the positions of the ith atom at the time origin and new time t, respectively. To improve the statistics, the MSD is smoothed by averaging MSD(t, t₀) using multiple time origins t₀ by Equation (4)(4) $\mathrm{M}\mathrm{S}\mathrm{D}\left(t,t_0\right)=\frac{1}{n}\sum _i{\left|r_i\left(t+t_0\right)-r_i\left(t_0\right)\right|}^2$(4) .

(4) $\mathrm{M}\mathrm{S}\mathrm{D}\left(t,t_0\right)=\frac{1}{n}\sum _i{\left|r_i\left(t+t_0\right)-r_i\left(t_0\right)\right|}^2$(4)

The interval of 2 fs for each cycle of atomic movements was used. For the first 1000 cycles, the temperature was elevated from 573 K to the target temperature of 1873 K. Then, FPMD calculations were performed on the basis of Nose thermostat [43,44] with a canonical ensemble for 150,000 cycles taking the average over time from 4 to 300 ps to obtain the MSD with enough movements of O atoms without considering the first 2000 cycles. D is averaged for three trials of FPMD, with different initial positions of V_O²⁺ at O(I), O(II), and O(III) sites.

3. Results and discussion

3.1. Crystal structure

Figure 2(a) shows the O–Ln and O–Ga bond lengths of the 15 types of Ln₃GaO₆ in this study. We separate the two types of O–Ln according to the CN of O. Calculated and experimental values [15], including Ln–O and Ga–O bond lengths, are provided in Supplementary Table S1. Note that when we obtain the mean bond lengths, the O–Ln bond lengths are counted as the distances between O and its three or four neighboring Ln, while the Ln–O bond lengths are counted as the distances between Ln and its seven neighboring O. We only considered the bond lengths among the distances between O and cations within a cutoff radius of 1.3 × the maximum peak obtained from the radial distribution functions (see Supplementary Figure S1). Our computational results obtained by both the GGA+U/w.f and GGA/w.o.f methods are consistent with the experimental reports. The computed bond lengths are slightly longer than the experimental bond lengths with an error of 3.5%. The differences of the bond lengths between the GGA+U/w.f and GGA/w.o.f methods are very small.

Figure 2. (a) Mean bond lengths of 15 Ln₃GaO₆ and (b) volume (eV/atom, left vertical axis) of 15 Ln₃GaO₆ and Ln₂O₃, and Ln³⁺ radius (right vertical axis). The experimental values of volume and the radius of Ln³⁺ are from Refs. [15] and [45], respectively.

In general, the Ln–O bond lengths decrease with increasing atomic numbers of Ln. For the GGA+U/w.f method, the O–La (CN of O = 4) and O–La (CN of O = 5) bond lengths are longer than O–Lu (CN of O = 4) and O–Lu (CN of O = 5) bond lengths by 0.26 and 0.29 Å, respectively. However, the O–Ga bond lengths show a little variation with the atomic number of Ln. The O–Ga bond lengths of La₃GaO₆ and Lu₃GaO₆ differ from each other only by 0.02 Å.

Figure 2(b) shows the volumes of Ln₃GaO₆ and the radii of Ln³⁺ (r_Ln3+). The calculated and experimental values [15] of volumes and lattice constants are mutually consistent, as shown in Supplementary Table S2.

As with the O–Ln bond lengths, in general, the volumes of Ln₃GaO₆ decrease with increasing atomic numbers of Ln. In addition, these values have a strong correlation with r_Ln3+. Considering that the O–Ga bond lengths change little among Ln₃GaO₆, these results suggest that the size of Ln₃GaO₆ depends on r_Ln3+.

In addition, we also compare the volumes of Ln₃GaO₆ and binary Ln₂O₃ (space group Ia-3 in the cubic crystal system) to determine the volume difference due to the existence of Ga as shown in Figure 2(b). As with Ln₃GaO₆, the decreasing tendency of the volume of Ln₂O₃ with increasing atomic number of Ln is found, which is consistent with a previous study by Petit et al. [18]. The volumes per atom of Ln₃GaO₆ are smaller than those of Ln₂O₃ by 9–15%. This may be ascribed to the considerably smaller ionic radius of Ga³⁺ (0.62 Å) than those of Ln³⁺ (0.86–1.03 Å) [45].

3.2. Electronic structure and location of f-levels

Before proceeding to our computational results of Ln₃GaO₆, it is worth reviewing previous studies on the locations of localized f-levels inside or outside of the E_g of Ln₂O₃. Without f-levels, HOMO and LUMO are mainly generated by O 2p and Ln 5d orbitals, respectively. However, because the localized f-levels locate inside E_g(O 2p–Ln 5d), the E_g values and the characteristics of orbitals forming HOMO and LUMO strongly depend on Ln and its f-orbital filling. To correct the location of f-levels, Jiang et al. [19,20] employed GW-level calculations based on the wavefunctions of the LDA+U method and obtained similar E_g values to experimental values. Gillen et al. [21] used hybrid functionals, namely, HSE03 and HSE06, and screened exchange local-density-approximations (sx-LDA). Altman et al. [46] observed the dependence of the location of f-levels on the type of Ln by X-ray absorption spectroscopy. These studies show a common trend: La₂O₃, Gd₂O₃, and Lu₂O₃, which have the earliest, the middle, and the last f-numbers of Ln, respectively, have less dependence on E_g(O 2p–Ln 5d) for the f-levels because their f-levels are relatively far from the O 2p-HOMO and Ln 5d-LUMO. Ce₂O₃ and Tb₂O₃, which have one larger f-number of Ln than La₂O₃ and Gd₂O₃, respectively, have occupied f-levels that are located much higher than those of La₂O₃ and Gd₂O₃. Then, from Ce₂O₃ to Gd₂O₃, and from Tb₂O₃ to Lu₂O₃, the locations of the occupied and unoccupied f-levels monotonically get shifted to lower levels with increasing atomic number of Ln. In the result, Ce₂O₃ and Pr₂O₃ have E_g(Ln 4f–Ln 5d) with the occupied f-levels located above the O 2p levels, whereas Ln₂O₃ (Ln = Pm, Sm, Eu, Tm, or Yb) have E_g(O 2p–Ln 4f) with unoccupied f-levels located below Ln 5d levels. For the other Ln₂O₃, the locations of the f-levels are shallow or out of E_g(O 2p–Ln 5d).

GW-level calculations can yield more accurate electronic structures; however, with current computational resources, it is not easy to systematically apply them to various materials that have large unit cells. Alternatively, we analyzed the location of the f-levels using two steps; first, we directly analyzed the f-levels using the projected electronic DOSs (PDOSs) by the GGA+U/w.f method. Then, we compared the E_g values of the GGA/w.o.f and GGA+U/w.f methods.

We analyzed whether the f-levels are formed inside E_g using the PDOSs obtained by the GGA+U/w.f method as shown in Figure 3. Without the f-levels, the HOMO and LUMO are mainly generated by O 2p and Ga 4s orbitals, respectively, which differs from that of Ln₂O₃. Unoccupied Ln 5d-levels are slightly higher than the Ga 4s-LUMO (also see the PDOS obtained by the GGA/w.o.f method in Supplementary Figure S2).

Figure 3. Electronic PDOS of 15 Ln₃GaO₆ obtained by the GGA+U/w.f method. Only Ln 4f (red) and 5d (yellow), Ga 4s (blue), O 2p (black) levels are drawn for simplicity. Energy is shifted by the computational Fermi level. Occupied and unoccupied levels are drawn as solid and dashed lines, respectively. Positive and negative values of PDOS denote the up-spin and down-spin states of electrons, respectively. For easier viewing, the unoccupied La 5d, Ga 4s, and O 2p levels are 2, 30, and 2 times magnified, respectively, whereas the Ln 4f-levels are 0.15 times decreased.

For La₃GaO₆, the unoccupied and occupied f-levels are generated far from the O 2p-HOMO and Ga 4s-LUMO, respectively, and spin polarizations are almost compensated. For Ce₃GaO₆, the occupied f-levels at the up-spin state are clearly located over the O 2p and form the HOMO decreasing the E_g value. Then, from Pr₃GaO₆ to Gd₃GaO₆, the occupied and unoccupied f-levels at the up-spin states move downward together with increasing atomic number of Ln. Gd₃GaO₆ has a clear E_g(O 2p–Ga 4s) with the occupied f-levels at the up-spin state and unoccupied f-levels at the down-spin state far from the O 2p-HOMO and Ga 4s-LUMO, respectively. For Tb₃GaO₃6, the occupied f-level at the down-spin state is generated just below the O 2p-HOMO. Then, from Tb₃GaO₆ to Lu₃GaO₆, the occupied and unoccupied f-levels at the down-spin state move downward again with increasing atomic number of Ln. Finally, Lu₃GaO₆ has all occupied f-levels under O 2p-LUMO with compensated spin polarizations.

In addition, we compared the E_g values obtained by the GGA+U/w.f and GGA/w.o.f methods. Figure 4 shows the E_g values of the GGA+U/w.f and GGA/w.o.f methods. The E_g(O 2p–Ga 4s) obtained by the GGA/w.o.f method are similar for all Ln₃GaO₆ with values of 3.47–3.83 eV. On the other hand, the E_g values obtained by the GGA+U/w.f method significantly change depending on the types of Ln because of the location of f-levels. As suggested by Lal and Gaur [47], we can classify Ln₃GaO₆ into four groups according to the types of orbitals that form the LUMO and HOMO. We also consider the numerical difference of E_g between the GGA+U/w.f and GGA/w.o.f methods (|ΔE_g|).

Figure 4. E_g values of 15 Ln₃GaO₆ obtained by the GGA+U/w.f and GGA/w.o.f methods. The E_g values of eight Ln₃GaO₆ (Ln = La, Nd, Gd, Tb, Ho, Dy, Er, or Lu) are not significantly different from each other because occupied or unoccupied f-levels are not formed deeply inside E_g(O 2p–Ga 4s).

We summarize the classification of Ln₃GaO₆ in Figure 5. For the first group [Figure 5(a)], Ln₃GaO₆ (Ln = La, Nd, Gd, Tb, Dy, Ho, Er, or Lu) have E_g(O 2p–Ga 4s), which are less affected by the location of the f-levels (|ΔE_g| < 0.3 eV). For the second group [Figure 5(b)], Ce₃GaO₆ and Pr₃GaO₆ have smaller E_g(Ln 4f–Ga 4s) obtained by the GGA+U/w.f method than those obtained by the GGA/w.o.f method owing to the occupied f-levels over the O 2p-LUMO. For the third group [Figure 5(c)], Ln₃GaO₆ (Ln = Pm, Sm, Eu, Tm, or Yb) also have smaller E_g(O 2p–Ln 4f) obtained by the GGA+U/w.f method than those obtained by the GGA/w.o.f method because of the unoccupied f-levels below the Ga 4s-LUMO inside E_g. These computational results suggest that the f-electrons should be considered as valence electrons for the second and third groups of Ln₃GaO₆. For the last group [Figure 5(d)], we do not find any Ln₃GaO₆, which has both occupied and unoccupied f-levels inside E_g(O 2p–Ga 4s). The dependence of the location of f-levels for Ln₃GaO₆ on the type of Ln is generally consistent with that for Ln₂O₃, which is confirmed by GW-level calculations [19,20] or experiments [46].

Figure 5. Schematic for the classification of Ln₃GaO₆ according to the location of occupied (red solid line) and unoccupied (red dashed line) f-level of Ln inside E_g(O 2p–Ga 4s). Black solid and blue dashed curves denote the HOMO (O 2p) and LUMO (Ga 4s), respectively. (a) The first group of Ln₃GaO₆ has E_g(O 2p–Ga 4s) less affected by the location of f-levels of Ln. (b) The second group of Ln₃GaO₆ has E_g(Ln 4f–Ga 4s). (c) The third group of Ln₃GaO₆ has E_g(O 2p–Ln 4f). (d) The last group of Ln₃GaO₆ has E_g(Ln 4f – Ln 4f), but they are not found.

Different features of f-levels depending on the type of Ln of Ln₃GaO₆ suggest different types of applications. The deep unoccupied f-levels of the third group of Ln₃GaO₆ (Ln = Pm, Sm, Eu, Tm, or Yb) may have a smaller excitation energy between the HOMO and the unoccupied f-level than other Ln₃GaO₆. For example, (Gd_1−xEu_x)₃GaO₆ compounds can be used as phosphors [48]. On the other hand, the first group of Ln₃GaO₆ (Ln = La, Nd, Gd, Tb, Dy, Ho, Er, or Lu) may have common properties of wide-gap compounds, whose optical absorption, transport properties, and magnetic properties are not affected by the f-levels of Ln. Henceforth, we focus on the first group of eight Ln₃GaO₆ (Ln = La, Nd, Gd, Tb, Dy, Ho, Er, or Lu) using only the GGA/w.o.f calculations, which are much faster than the GGA+U/w.f calculations because the latter often suffers from slow convergence due to the localized f-levels and spin polarizations.

3.3. Energetic and dynamical stability

Figure 6 shows the energies of Ln₃GaO₆ (Ln = La, Nd, Gd, Tb, Dy, Ho, Er, or Lu) compared with those of two competing reference states. For the competing reference states, first, the coexisting state of binary sesquioxides of Ln₂O₃ and Ga₂O₃ (space group C2/m) is employed. Second, the coexisting state of Ln₄Ga₂O₉ (space group P2₁/c) and Ln₂O₃ is employed. According to phase diagrams in experimental studies [49–51], the coexisting state of La₄Ga₂O₉ (space group P2₁/c) and La₂O₃ can be formed at the composition of La₃GaO₆.

Figure 6. Energy differences (eV/atom) of eight Ln₃GaO₆ (Ln = La, Nd, Gd, Tb, Ho, Dy, Er, or Lu) with respect to their competing reference states, namely, the coexisting states of ¾(Ln₂O₃) + ¼(Ga₂O₃) and ¾(Ln₄Ga₂O₉) + ¼(Ln₂O₃).

In the result, seven Ln₃GaO₆ except for Lu₃GaO₆ are energetically more stable than their competing reference states of Ln₂O₃ and Ga₂O₃. The absolute energy differences tend to decrease with increasing atomic number of Ln. Lu₃GaO₆ is energetically slightly less stable than their corresponding coexisting states of Ln₂O₃ and Ga₂O₃.

Compared with the coexisting state of Ln₄Ga₂O₉ and La₂O₃, the seven Ln₃GaO₆ series, except for La₃GaO₆, are energetically more stable. The energy of La₃GaO₆ is slightly higher (0.002 eV/atom) than that of the coexisting state of La₄Ga₂O₉ and La₂O₃. This small energy difference implies a possibility of synthesis of La₃GaO₆ as a single phase depending on the synthesis process, such as extrinsic doping and changing of thermodynamic variables.

To summarize, the six Ln₃GaO₆ (Ln = Nd, Gd, Tb, Dy, Ho, or Er) are more stable than the two reference states. This computational result on the energetic stability is consistent with the experimental reports. These six oxides on top of Sm₃GaO₆ and Eu₃GaO₆ were synthesized as single phases by sintering binary sesquioxide powders at 1400 °C [15] or by crystallization of reactive high temperature hydroxide melt [17].

We also investigate the dynamical stability of the eight Ln₃GaO₆ (Ln = La, Nd, Gd, Tb, Dy, Ho, Er, or Lu) by confirming whether imaginary phonon frequencies on the phonon DOS exist or not. As Figure 7 shows, no imaginary phonon frequencies are found for all the eight Ln₃GaO₆, although we confirmed that two of them (La₃GaO₆ and Lu₃GaO₆) are not energetically in the ground-state. The dynamical stabilities of these oxides suggest that they can also exist as metastable states if they can be once synthesized.

Figure 7. Phonon DOS of eight Ln₃GaO₆ (Ln = La, Nd, Gd, Tb, Ho, Dy, Er, or Lu). The vertical dashed line denotes zero phonon frequency. No imaginary phonon frequencies are found, which indicates that these eight oxides are dynamically stable.

To investigate the stability of these Ln₃GaO₆, we consider Pauling’s radius ratio rule [52]. This rule suggests that favorable CNs of cations in the crystal structure depend on the ionic radius ratio (r_Ln3+/r_O2−) of cations (Ln in this study) and anions (O in this study, 1.4 Å [53]). According to this rule, the conditions of 0.41 ≤ r_Ln3+/r_O2− < 0.59 and 0.59 ≤ r_Ln3+/r_O2− < 0.73 correspond to favorable CNs of Ln of six and seven, respectively [54]. As shown in Figure 2(b), in general, r_Ln3+ decreases with increasing atomic number of Ln with values of 0.86 (Lu)–1.03 (La) Å. These values can be converted to r_Ln3+/r_O2− of 0.61 (Lu)–0.74 (La). This suggests that as r_Ln3+ is larger (the atomic number of Ln is smaller), it prefers to form the CN of cation of seven rather than of six. The CNs of Ln of Ln₃GaO₆ and Ln₂O₃ are seven and six, respectively. Therefore, as r_Ln3+/r_O2− value is larger (the atomic number of Ln is smaller), Ln₃GaO₆ is preferred to be formed rather than binary oxides, which can be confirmed by the energy differences between Ln₃GaO₆ and binary oxides as shown in Figure 6. Lu₃GaO₆ becomes less stable than the binary oxide because its r_Ln3+/r_O2− (0.61) is decreased to near the limit of 0.59.

3.4. Formation energies of V_O

Experimental reports [13,14] on Nd₃GaO₆ and Gd₃GaO₆ doped with acceptor dopants (aliovalent cations with lower valences) such as Sr and Ca (for example, Nd_3−xSr_xGaO_6−2/x) imply that Ln₃GaO₆ have the diffusion mechanism of O based on V_O²⁺. Therefore, it is worthwhile to investigate the possibility of generating V_O²⁺ assuming such type of doping condition with acceptor dopants in Ln₃GaO₆. When the negatively charged defects generated by substitution of acceptor dopants for Ln³⁺ or Ga³⁺ are introduced, holes can be increased with shifting the Fermi level from the center of band gap to the HOMO and positively charged defects such as V_O²⁺ (for semiconductors, these defects are sometimes called ‘hole killers’ [55]) can be generated to compensate for the charge imbalance. In this study, we mainly focus on the case where the Fermi level is shifted to HOMO because a few fractions of extrinsic dopants are commonly substituted (high concentration of approximately 10²⁰–10²²/cm³) for constituent cations in oxygen-ion conductors [2,13,14].

The E_v values for the eight Ln₃GaO₆ (Ln = La, Nd, Gd, Tb, Dy, Ho, Er, or Lu) are summarized in Table 1. We employed μ_O at the pressure and temperature of 1 atm and 1500 K, respectively, based on the sintering conditions for syntheses of Ln₃GaO₆ (1250–1450 °C) [13–15]. The μ_O value at a certain temperature (T) and pressure (p) can be obtained as the half energy of O₂ gas molecule minus 1.75 eV by the following equation:

(5) $\text{{mu \_{rm{O}}}left({T,p} right) = {1 over 2}Eleft({{{rm{O}}\_2}} right) + {mu \_{rm{O}}}left({T,{p^^circ }} right) + {1 over 2}kT{rm{ln}}left({{p over {{p^^circ }}}} right)}$(5)

Table 1. E_v of eight Ln₃GaO₆ (Ln = La, Nd, Gd, Tb, Ho, Dy, Er, or Lu). The computational results in this table are obtained using the GGA/w.o.f method. The Fermi level is located at the HOMO (VBM). The values are in eV.

O site		La3GaO6	Nd3GaO6	Gd3GaO6	Tb3GaO6	Dy3GaO6	Ho3GaO6	Er3GaO6	Lu3GaO6
	Eg	3.83	3.62	3.59	3.57	3.56	3.54	3.51	3.47
O(I)	Ev(VO^0)	2.89	2.91	2.92	2.91	2.91	2.91	2.91	2.91
	Ev(VO^+)	1.25	1.64	1.50	1.51	1.52	1.53	1.54	1.55
	Ev(VO^2+)	−0.49	−0.09	−0.01	0.00	0.01	0.02	0.04	0.04
	ε(2+/0)^a	1.69	1.50	1.46	1.46	1.45	1.45	1.44	1.43
O(II)	Ev(VO^0)	4.92	4.67	4.81	4.83	4.84	4.85	4.87	4.92
	Ev(VO^+)	1.65	2.00	1.86	1.89	1.90	1.92	1.94	1.99
	Ev(VO^2+)	−0.77	−0.44	−0.40	−0.39	−0.39	−0.39	−0.38	−0.40
	ε(1+/0)^a	3.27	2.67	2.95	2.94	2.94	2.93	2.93	2.93
	ε(2+/1+)^a	2.42	2.44	2.26	2.28	2.29	2.31	2.32	2.39
O(III)	Ev(VO^0)	2.68	2.82	2.85	2.86	2.86	2.86	2.86	2.86
	Ev(VO^+)	1.07	1.53	1.40	1.41	1.41	1.42	1.43	1.42
	Ev(VO^2+)	−0.71	−0.25	−0.18	−0.17	−0.17	−0.17	−0.14	−0.17
	ε(2+/0)^a	1.69	1.54	1.52	1.51	1.51	1.51	1.50	1.52
O(IV)	Ev(VO^0)	2.79	2.86	2.89	2.89	2.89	2.89	2.90	2.90
	Ev(VO^+)	1.50	1.90	1.81	1.83	1.84	1.85	1.88	1.90
	Ev(VO^2+)	0.55	1.00	1.16	1.20	1.23	1.26	1.31	1.36
	ε(1+/0)^a	1.29	0.96	1.08	1.06	1.05	1.04	1.02	1.00
	ε(2+/1+)^a	0.95	0.90	0.65	0.63	0.61	0.59	0.57	0.54

^a Locations are counted from the O 2p-HOMO.

where μ_O(T, p°) can be taken from thermodynamic tables (see the values in other temperatures in Supplementary Table S3) at the standard pressure of 1 atm [56,57]. In addition, the dielectric constants of these eight La₃GaO₆, which are used for the image charge correction [34,35], range from 14.6 to 19.4 as summarized in Supplementary Table S4.

First, we analyze E_v as a function of the Fermi level of La₃GaO₆ as shown in Figure 8. The three types of E_v(V_O⁰) for the O(I), O(III), and O(IV) sites, which have both O–Ln and O–Ga bonds, are all similar at 2.68–2.89 eV, whereas the E_v(V_O⁰) from the O(II) site, which has only O–Ln bonds, is much larger at 4.92 eV. When the Fermi level is at the center of E_g or shifted into the LUMO, the E_v(V_O⁰) values for the O(I), O(III), and O(IV) sites are much lower than those of E_v(V_O⁺) and E_v(V_O²⁺). In contrast, when the Fermi level is shifted into the HOMO (assuming that acceptor dopants are introduced to substitute for cations, so that positive charges try to compensate for the charge imbalance), the E_v(V_O²⁺) value for the O(II) site is lower than those of the other sites by 0.06–1.04 eV at the same Fermi level, although the other sites also have lower E_v(V_O²⁺) than E_v(V_O⁺) and E_v(V_O⁰) at this region. When the Fermi level is near the HOMO, the E_v(V_O²⁺) values for the three types of O sites except for O(IV) site become even negative values, which imply a spontaneous generation of V_O²⁺.

Figure 8. E_v as a function of Fermi level of La₃GaO₆. μ_O is set to be ½E(O₂) − 1.75 eV to describe the pressure and temperature at 1 atm and 1500 K, respectively. The other Ln₃GaO₆ (Ln = Nd, Gd, Tb, Dy, Ho, Er, or Lu) series showed a similar trend (see Supplementary Figures S3). The gradient of the lines denotes the charge of the defect.

The thermodynamic transition points, which the charge of V_O changes, ε(2+/1+), ε(1+/0), and ε(2+/0) are all deeply formed. The highest thermodynamic transition point, which is ε(1+/0) for the O(II) site, forms a deep level of ~0.6 eV below the LUMO.

We compare the E_v values among the eight Ln₃GaO₆. The E_v values as a function of the Fermi level for the eight Ln₃GaO₆ are also shown in Supplementary Figure S3. The trends in the change of stable charge of V_O according to the Fermi level are similar to that of La₃GaO₆ as abovementioned. When the Fermi level is located at the center of E_g, the most stable is V_O⁰ at the O(III) site with E_v values of 2.68–2.86 eV, which are similar to each other for Ln₃GaO₆. When the Fermi level is close to the LUMO and HOMO, V_O⁰ and V_O²⁺ are energetically favorable, respectively. The most energetically stable V_O²⁺ is preferred for the O(II) site, which has only O–Ln bonds. As shown in Supplementary Figure S4, the top of the HOMO is mainly formed by the O(II) site. Therefore, the lower E_v value for the O(II) site than those of the other O site is reasonable because the E_v value indicates the energy required for breaking the O–cation bonds and removing O out of the oxide. When the Fermi level is located at the center of E_g, the E_v values of V_O²⁺ at the O(II) sites range from 3.06 to 3.19 eV, which are also similar to each other for Ln₃GaO₆. When the Fermi level is located at the HOMO, the E_v value of V_O²⁺ at the O(II) site of La₃GaO₆ ranges from 0.33 to 0.38 eV. This is lower than those of the other seven Ln₃GaO₆ because the E_g value of La₃GaO6 is ~0.2 eV larger so that the decrease in the E_v value with the gradient of −2 is additional in the range of larger E_g. The E_v(V_O²⁺) values for the O(I), O(II), and O(III) sites of Ln₃GaO₆ are negative or close to zero, which implies a spontaneous generation of V_O²⁺. This suggests that Ln₃GaO₆ can easily generate V_O²⁺ when acceptor dopants substitute for Ln³⁺ or Ga³⁺. In addition, when the pressure of O becomes lower, the E_v value for V_O further decreases as indicated by Equations (1(1) $E_v^q=E\left(\mathrm{L}{\mathrm{n}}_3\mathrm{G}\mathrm{a}{\mathrm{O}}_6:V_{\mathrm{O}}^q\right)-E\left(\mathrm{L}{\mathrm{n}}_3\mathrm{G}\mathrm{a}{\mathrm{O}}_6\right)+{\mu }_{\mathrm{O}}+q\left(E_{\mathrm{H}\mathrm{O}\mathrm{M}\mathrm{O}}+E_{Fermi}\right)+\mathrm{\Delta }{E}_{\mathrm{L}\mathrm{Z}}$(1) ) and (5(5) $\text{{mu \_{rm{O}}}left({T,p} right) = {1 over 2}Eleft({{{rm{O}}\_2}} right) + {mu \_{rm{O}}}left({T,{p^^circ }} right) + {1 over 2}kT{rm{ln}}left({{p over {{p^^circ }}}} right)}$(5) ). Positive charges such as holes and V_O²⁺ compete with each other to compensate for the limited concentration of negative charged defects. Therefore, the generation of V_O²⁺ can be increased by decreasing the pressure of the oxygen source.

3.5. Migration barrier energy (E_m) of V_O

There are various migration paths of V_O²⁺ in Ln₃GaO₆. Instead of investigating all types of paths for all types of Ln₃GaO₆, we chose La₃GaO₆ as a representative sample and obtained E_m for all possible migration paths considering the nearest-neighboring O sites. Table 2 presents the E_m values for possible migration paths of La₃GaO₆. We named the path after the type of V_O²⁺ site in the initial state for the CI-NEB in the order of O–O distances. The O(II) site is energetically the most stable for V_O²⁺, so that this site is considered as a reference. When we consider the E_m value on the migration paths between two O sites other than the O(II) site, we added the energy difference (ΔE_ref) between the more stable O site (between initial and final states) on the migration path and the O(II) site.

Table 2. E_m for various paths of eight La₃GaO₆. The computational results in this table are obtained using the GGA/w.o.f method. IS, TS, and FS denote the initial, transition, and final states for CI-NEB, respectively.

VO^2+ site of IS^a	Path	VO^2+ site of FS	E(FS) − E(IS) (eV)	O site distance (Å)	Em (eV)^b	E(TS) − E(IS) with VO^2+[O(II)] (eV)^c	Index
O(I)	I-1	O(I)	0.00	2.880	0.62	0.89
[+0.27 eV]	I-2	O(III)	−0.16	2.950	0.51	0.63
	I-3	O(II)	−0.27	2.987	1.03	1.03
	I-4	O(IV)	1.02	3.187	1.24	1.51
	I-5	O(IV)	1.02	3.298	1.89	2.16
	I-6	O(IV)	1.02	3.304	2.12	2.38
	I-7	O(I)	0.00	3.407	1.02	1.29
	I-8	O(II)	−0.27	3.462	0.94	0.94
O(II)	II-1	O(II)	0.00	2.867	0.59	0.59	The lowest Em
[+0.00 eV]	II-2	O(I)	0.27	2.987			Same as path I-3
	II-3	O(II)	0.00	3.004	1.02	1.02
	II-4	O(III)	0.11	3.269	0.76	0.76
	II-5	O(I)	0.27	3.462			Same as path I-8
O(III)	III-1	O(I)	0.16	2.950			Same as path I-2
[+0.11 eV]	III-2	O(III)	0.00	3.048	0.74	0.85
	III-3	O(IV)	1.18	3.086	1.67	1.78
	III-4	O(II)	−0.11	3.269			Same as path II-4
	III-5	O(IV)	1.18	3.271	1.66	1.77
O(IV)	IV-1	O(III)	−1.18	3.086			Same as III-3
[+1.18 eV]	IV-2	O(I)	−1.02	3.187			Same as I-4
	IV-3	O(III)	−1.18	3.271			Same as path III-5
	IV-4	O(I)	−1.02	3.298			Same as path I-5
	IV-5	O(I)	−1.02	3.304			Same as path IV-5

^a ΔE_ref, which can be defined as E(IS) − E(IS) with V_O²⁺[O(II)], with a unit of eV are in the brackets.

^b Larger value between E(TS) − E(IS) and E(TS) − E(FS).

^c E(IS) with V_O²⁺[O(II)] is the most stable among four types of O sites.

We identify five migration paths with the lowest E_m value in La₃GaO₆, namely, paths II-1, II-4, II-5 (I-8), III-1 (I-2), and III-2, which are between O(II)–O(II), O(II)–O(III), O(II)–O(I), O(III)–O(I), and O(III)–O(III) sites, respectively. Then, we also compute the E_m values on the five migration paths for the other seven Ln₃GaO₆ (Ln = Nd, Gd, Tb, Dy, Ho, Er, or Lu). These five types of migration paths are displayed in Supplementary Figure S5. Because the O(IV) site has ΔE_ref of ~1 eV, which is higher than various E_m, V_O²⁺ may rarely reach this site.

Figure 9(a) shows the E_m values for the five migration paths of the eight Ln₃GaO₆. As abovementioned, we additionally considered the ΔE_ref for the paths III-1 and III-2 because V_O²⁺ at the O(III) or O(I) sites are less stable than at the O(II) site. Note that those energy differences are slightly larger than those summarized in Table 1 because the image charge correction is not applied for the CI-NEB method.

The lowest E_m value among the five migration paths is found to be path II-1, which is between O(II)–O(II) (see Supplementary Figure S5). In addition, path II-1 is delocalized and connected to the neighboring cell; therefore, this can be the dominant migration path. In contrast, paths II-4, II-5, and III-1 with higher (E_m+ΔE_ref) than the E_m value of path II-1 are localized and not connected; therefore, complex combinations of paths for migrating V_O²⁺ to reach the neighboring cell are required. Path III-2 is delocalized and connected to the neighboring cell; however, this path is formed between the less stable O(III) sites. To observe total ionic movements on various migration paths with time scale, other methods such as kinetic Monte Carlo [58,59] and molecular dynamics are necessary. We used FPMD in this study; thus, we will discuss the D_O that is obtained from the total ionic movements by FPMD in the next subsection.

La₃GaO₆ has the lowest E_m value on migration path II-1 among the eight Ln₃GaO₆. Then, as the atomic number of Ln increases, the E_m values of the eight Ln₃GaO₆ also increase. The dependence of E_m on the atomic number of Ln is basically not changed when larger computational cells are used (with a little reduction of less than 0.07 eV) as shown in Supplementary Figure S6. Migration path II-1 is neighbored with Ln; therefore, the space for the migration path is related to r_Ln3+. The E_m value on migration path III-1 also increases with the increase in the atomic number of Ln. On the other hand, the E_m values on the other three migration paths slightly change or slightly decrease with increasing atomic number of Ln. However, the migration of V_O²⁺ may be dominated by the lowest E_m; therefore, it can be expected that La₃GaO₆ might have the most active migrations of V_O²⁺ followed by Nd₃GaO₆.

To the best of our knowledge, binary Ln₂O₃ compounds are not known to be applicable as oxygen-ion conductors. To determine the effect of the existence of Ga, we compare the distance between O sites neighbored with only Ln as shown in Figure 9(b). As abovementioned, path II-1 is between the O(II) sites that only have the chemical bonds with only Ln. It is noticeable that the distances between the O(II) sites in the eight Ln₃GaO₆ are shorter than the minimum distances between the O sites in the counterpart Ln₂O₃ by 0.16–0.23 Å (6–8%). In this connection, the E_m values for the migration path between these O sites of the eight Ln₃GaO₆ are much lower than those of the counterpart Ln₂O₃ by 0.12–0.23 eV (17–21%). Considering the smaller volumes per atom of Ln₃GaO₆ than those of Ln₂O₃ (Figure 2(b)), we suggest that the advantage of participation of Ga in Ln₃GaO₆ is to shrink the migration distance for V_O and decrease E_m, which may result in more active diffusion of V_O.

3.6. Oxygen-ion diffusivity (D_O) and conductivity (σ_O)

D_O was calculated to investigate the total movement of O among the eight Ln₃GaO₆ possessing V_O²⁺ by FPMD. The MSDs of O for the computation of D_O are displayed in Supplementary Figure S7. To accelerate the O diffusion, we chose a high temperature of 1873 K. Here, 1873 K is a little lower temperature than the melting temperature reported in the La₂O₃–Ga₂O₃ phase diagram [49–51]. We confirm that MSDs of Ln and Ga, which are more than 1 Ų, are not found during the FPMD for the eight Ln₃GaO₆ (not shown).

To compare the computed σ_O with the experimental values at a lower temperature of 800 °C, we extrapolated σ_O using the D_O value as follows. The temperature dependence of D_O can be expressed by the Arrhenius equation:

(6) $\begin{array}{rl}& D_{\mathrm{O}}\left(T\right)=D_{\mathrm{O},0}exp\left(\frac{-E_a}{k_{B}T}\right)\mathrm{o}\mathrm{r}\\ & \mathrm{l}\mathrm{o}\mathrm{g}{D}_{\mathrm{O}}\left(T\right)=\mathrm{l}\mathrm{o}\mathrm{g}{D}_{\mathrm{O},0}-\frac{E_a}{2.303k_{B}T}\end{array}$(6)

where D_O,0 is a pre-exponential factor, k_B is the Boltzmann constant, and E_a is an activation energy. Here, for the E_a of equation 6(6) $\begin{array}{rl}& D_{\mathrm{O}}\left(T\right)=D_{\mathrm{O},0}exp\left(\frac{-E_a}{k_{B}T}\right)\mathrm{o}\mathrm{r}\\ & \mathrm{l}\mathrm{o}\mathrm{g}{D}_{\mathrm{O}}\left(T\right)=\mathrm{l}\mathrm{o}\mathrm{g}{D}_{\mathrm{O},0}-\frac{E_a}{2.303k_{B}T}\end{array}$(6) , we employed the lowest E_m value on path II-1 obtained by the CI-NEB method, which will be discussed later. The D_O(T) at temperatures lower than 1873 K can be obtained by the following relationship:

(7) $D_{\mathrm{O}}\left(T\right)=D_{\mathrm{O}}\left(1873\mathrm{K}\right)exp\left[\frac{-E_a}{k_B}\left(\frac{1}{T}-\frac{1}{1873}\right)\right]$(7)

The σ_O(T) value can be obtained from D_O(T) by the Nernst–Einstein equation at the dilute level as follows [60]:

(8) ${\sigma }_{\mathrm{O}}\left(T\right)=\frac{D_{\mathrm{O}}\left(T\right){\left(ze\right)}^{2}c}{k_{B}T}=\frac{D_{\mathrm{O},0}{\left(ze\right)}^{2}c}{k_{B}T}exp\left(\frac{-E_a}{k_{B}T}\right)$(8)

where z is the valence of oxygen, e is the unit charge (1.6 × 10⁻¹⁹ C), and c is the ionic concentration of O in Ln₃GaO₆.

The main factors of D_O,0 are the jumping frequency Г, which is approximately 10⁻¹² Hz in oxides, and the squared distances between the sites α² [61]. We expect that the values of log D_O,0 in Equation (6)(6) $\begin{array}{rl}& D_{\mathrm{O}}\left(T\right)=D_{\mathrm{O},0}exp\left(\frac{-E_a}{k_{B}T}\right)\mathrm{o}\mathrm{r}\\ & \mathrm{l}\mathrm{o}\mathrm{g}{D}_{\mathrm{O}}\left(T\right)=\mathrm{l}\mathrm{o}\mathrm{g}{D}_{\mathrm{O},0}-\frac{E_a}{2.303k_{B}T}\end{array}$(6) for these eight Ln₃GaO₆ may be similar to each other considering the same crystal structure and small differences of α² as shown in Figure 9(b). Then, log D_O at the same temperature may depend on E_a, which is assumed to be the lowest E_m. Figure 10(a) shows the relationship of log D_O at 1873 K and the lowest E_m (from Figure 9) of Ln₃GaO₆. A strong linear correlation between these two properties is found with Pearson’s linear correlation coefficient (r_p) of −0.99 [62]. It is noticeable that La₃GaO₆ has the highest D_O at this temperature, followed by Nd₃GaO₆.

Figure 9. (a) E_m for various migration paths of eight Ln₃GaO₆ (Ln = La, Nd, Gd, Tb, Ho, Dy, Er, or Lu) and (b) the lowest E_m of the eight Ln₃GaO₆ and their counterpart Ln₂O₃ (left vertical axis) and the distances between O sites for the migration path in the crystal structure (right vertical axis).

Figure 10. (a) Relationship between D_O obtained by FPMD at 1873 K and the lowest E_m obtained by CI-NEB, and (b) predicted σ_O as a function of T based on D_O at 1873 K for eight Ln₃GaO₆ (Ln = La, Nd, Gd, Tb, Ho, Dy, Er, or Lu). Solid symbols denote computed or experimentally measured values from references. Dashed lines denote predicted values. The concentration of V_O²⁺ of our computation is ~2.1% from the computational cell of Ln₂₄Ga₈O₄₇. The concentration of V_O²⁺ of Nd_2.955Sr_0.045GaO_5.9775 [13] and Gd_2.9Sr_0.1GaO_5.95 [14] are ~0.4 and ~0.8%, respectively. σ_O of 10⁻² S/cm at 1000 K is described to show a reference level of YSZ [2]. The D_O values are averaged from three FPMD trials and error ranges are described between the minimum and maximum values.

Figure 10(b) shows the predicted σ_O of Ln₃GaO₆ as a function of temperature obtained by the above assumption. The experimental σ_O (0.7 × 10⁻² S/cm and 0.3 × 10⁻² S/cm) of Nd₃GaO₆ and Gd₃GaO₆ (at dry Ar condition) are shown as the highly reported value between various dopants at 1073 K (800 °C), which are obtained at the composition of Nd_2.955Sr_0.045GaO_5.9775 [13] and Gd_2.9Sr_0.1GaO_5.95 [14], respectively. Our calculated σ_O values of Nd₃GaO₆ and Gd₃GaO₆ at 800 °C are slightly higher than the experimental values, but lie within the same order of magnitude (approximately 10⁻³–1 × 10⁻² S/cm).

It is worthwhile to discuss that the differences between our prediction and experimental results might be due to the following reasons. The type and amount of dopants affect the concentrations of V_O²⁺ and E_a for the diffusions of oxygen atoms. We used the concentration of V_O²⁺ of ~2.1% from the computational cell of Ln₂₄Ga₈O₄₇, while the concentration of V_O²⁺ of Nd_2.955Sr_0.045GaO_5.9775 [13] and Gd_2.9Sr_0.1GaO_5.95 [14] are ~0.4 and ~0.8%, respectively (assuming that the Sr doping perfectly generates V_O²⁺). In addition, Ref. [14] reported that the Ca-doped Gd₃GaO₆ (Gd_3−xCa_xGaO₆) shows one order lower σ_O of 0.04 × 10⁻² S/cm compared with the counterpart substituted by Sr (Gd_3−xSr_xGaO₆), which shows a σ_O of 0.3 × 10⁻² S/cm. This may be because E_a can be affected by the type of dopant. Previous studies on doped ceria [63,64] suggest that the atomic radius of the dopant is an important factor for changing E_a and finding the optimum type of dopant. We expect that optimum dopants may be different from each other for Ln₃GaO₆ because of different values of r_Ln3+.

The experimental σ_O also depends on the microstructure of the experimental samples (for example, densification and particle size) and the measuring environment (for example, wet or dry). Ref. [14] reported that a proton conductivity of 1 × 10⁻³ S/cm was also detected for doped Gd₃GaO₆ when measured at a slightly lower temperature of 600 °C. When the σ_O of Gd_2.9Sr_0.1GaO_5.95 was measured at wet Ar condition and 800 °C, σ_O reportedly increased to 0.9 × 10⁻² S/cm from 0.3 × 10⁻² S/cm measured at dry Ar condition. This implies the partial existence of proton conductivity in the total conductivity. We also confirm that protons can be stabilized depending on the environmental variables, such as the type and pressure of source gas, and temperature; therefore, the proton conductivity can be found for other Ln₃GaO₆ (see Supplementary Section S.2). When we introduce negatively charged defects by substituting acceptor dopants for Ln³⁺ or Ga³⁺, generation of positively charged defects (proton and V_O²⁺), which may be more dominant than holes considering negative defect formation energies for these defects and large band-gaps of the eight Ln₃GaO₆ compounds, compete with each other to compensate for the charge imbalance. Takahashi et al., suggested that the stability of protons and V_O²⁺ is under a trade-off relationship [65]. Therefore, we need to suppress the incorporation of protons to prevent interruption of V_O²⁺ generation for a suitable application of oxygen-ion conductor. As discussed in Supplementary Section S.2, the increase of generation of V_O²⁺ (decrease of E_v) and decrease of incorporation of protons (increase of protonic defect formation energy) can be achieved by increasing temperature or by decreasing oxygen gas and water vapor pressures.

As with the trend of D_O at 1873 K or the lowest E_m and the atomic number of Ln, the predicted σ_O at 1000 K tends to decrease with increasing atomic number of Ln. La₃GaO₆ and Nd₃GaO₆ show larger σ_O at the same temperature than the other Ln₃GaO₆. The predicted σ_O values of the eight Ln₃GaO₆ at 1000 and 1073 K (800 °C) are in the range of 0.1–1.5 × 10⁻² and 0.3–2.2 × 10⁻² S/cm, respectively. This computational result suggests a possibility that these eight Ln₃GaO₆ can be expected to show a similar order of magnitude of σ_O in the experiment to a common reference value of ~1 × 10⁻² S/cm of YSZ at 1000 K [2]. We expect that La₃GaO₆ may achieve a higher σ_O than the experimentally reported Nd₃GaO₆ and Gd₃GaO₆ if successfully synthesized with the same amount of V_O.

As widely known, E_a is generally obtained using the gradient of Equation (6)(6) $\begin{array}{rl}& D_{\mathrm{O}}\left(T\right)=D_{\mathrm{O},0}exp\left(\frac{-E_a}{k_{B}T}\right)\mathrm{o}\mathrm{r}\\ & \mathrm{l}\mathrm{o}\mathrm{g}{D}_{\mathrm{O}}\left(T\right)=\mathrm{l}\mathrm{o}\mathrm{g}{D}_{\mathrm{O},0}-\frac{E_a}{2.303k_{B}T}\end{array}$(6) with several D_O’s at different temperatures. However, poor statistical data such as insufficient MSD of O are found at temperatures lower than 1873 K. Alternatively, we used the lowest E_m as E_a. In this case, we need to assume that the diffusion mechanism and phase of Ln₃GaO₆ are not changed between the computed and target temperatures, and the diffusion of O is mainly found on the suggested migration path with the lowest E_m. We confirm that these eight Ln₃GaO₆ are dynamically stable so that D_O is exclusively obtained for the investigated crystal structures. In addition, we fixed one V_O in the computational cell; therefore, the diffusion mainly occurs due to V_O (in fact, an O atom migrates in the opposite direction with V_O). We also confirm that the main diffusion trajectory of O is similar to migration path II-1 with the lowest E_m as shown in Supplementary Figures S5 and S8.

Meanwhile, it is worth investigating how the predicted σ_O at 1000 K changes according to the E_a value because the diffusions of O may not always occur only on the dominant migration path with the lowest E_m. Assuming that the D_O at 1873 K is fixed because the FPMD already reflects the features of total movements of O, we investigate the change in predicted σ_O as a function of ΔE_a, which denotes a deviation from the lowest E_m for the E_a of Equations (6(6) $\begin{array}{rl}& D_{\mathrm{O}}\left(T\right)=D_{\mathrm{O},0}exp\left(\frac{-E_a}{k_{B}T}\right)\mathrm{o}\mathrm{r}\\ & \mathrm{l}\mathrm{o}\mathrm{g}{D}_{\mathrm{O}}\left(T\right)=\mathrm{l}\mathrm{o}\mathrm{g}{D}_{\mathrm{O},0}-\frac{E_a}{2.303k_{B}T}\end{array}$(6) )–(8(8) ${\sigma }_{\mathrm{O}}\left(T\right)=\frac{D_{\mathrm{O}}\left(T\right){\left(ze\right)}^{2}c}{k_{B}T}=\frac{D_{\mathrm{O},0}{\left(ze\right)}^{2}c}{k_{B}T}exp\left(\frac{-E_a}{k_{B}T}\right)$(8) ) as shown in Supplementary Figure S9. We confirm that even in the case where ΔE_a is 0.1 eV, most of the predicted σ_O values at 1000 K for these Ln₃GaO₆ keep the same order of magnitude as 10⁻³ S/cm (except for Lu₃GaO₆ with σ_O of 0.8 × 10⁻³ S/cm). It is found that one order of magnitude of σ_O at 1000 K is decreased when a relatively large value of ΔE_a of ~0.4 eV is increased. Although we exceedingly assume that the diffusions of O mainly occur on the migration path of V_O²⁺ with the highest E_m (0.94–1.06 eV) in Figure 9(a), ΔE_a is less than 0.4 eV. In an experiment, Ref. [13] reported that doped Nd₃GaO₆, namely, Nd_2.91Ca_0.09GaO_5.955 and Nd_2.955Sr_0.045GaO_5.9775 showed an E_a value of 0.97 and 1.06 eV, respectively. Ref. [14] reported that doped Gd₃GaO₆, namely, Gd_2.9Ca_0.1GaO_5.95 and Gd_2.9Sr_0.1GaO_5.95 showed an E_a value of 0.98 and 0.92 eV at dry Ar condition, respectively. Although there are some variables such as different concentrations of V_O²⁺ and existence of dopants, the differences between the lowest E_m and the E_a of the more optimum reported value in the experiment of Nd₃GaO₆ and Gd₃GaO₆ are approximately less than ~0.3 eV. Therefore, we expect that the prediction error of σ_O at 1000 K of our model may be less than one order of magnitude.

Before concluding, we suggest another merit of the investigation of Ln₃GaO₆. Recently, there have been many machine learning studies on exploring good functional materials [66–68]. Despite the efficient search by this type of study, one shortcoming is the insufficient number of training data that exclusively have the target variable in a wide range. Although our predicted σ_O values of Ln₃GaO₆ are not the highest among all the reported values, the moderate value of functional property is also useful for an accurate prediction model in a wider search space [68]. We expect that the reported properties of Ln₃GaO₆ in this study and other relevant future works in experiments may be able to contribute to the construction of abundant training data for the machine learning models in the search for oxygen-ion conductors.

4. Conclusion

We systematically investigated the structural and electronic structure of 15 Ln₃GaO₆ (space group Cmc2₁). In general, Ln₃GaO₆ have decreasing tendency of the Ln–O (or O–Ln) bond lengths and volumes with increasing the atomic number of Ln and decreasing r_Ln3+. Then, we distinguished eight Ln₃GaO₆ (Ln = La, Nd, Gd, Tb, Ho, Dy, Er, or Lu) that have E_g(O 2p–Ga 4s) from the others that have the localized f-levels inside E_g(O 2p–Ga 4s) with crosschecks on the previous studies on Ln₂O₃. These eight Ln₃GaO₆ are dynamically stable. Among them, six Ln₃GaO₆ (Ln = Nd, Gd, Tb, Ho, Dy, or Er) are also energetically stable than the competing reference states, which is consistent with experimental reports on successful syntheses of a single phase. La₃GaO₆ and Lu₃GaO₆ are energetically slightly less stable than their competing reference states; however, they might be stabilized within a certain process window considering dynamical stabilities and small energy differences from their competing reference states.

The properties related to V_O were also analyzed. We confirmed that Ln₃GaO₆ may easily generate V_O²⁺ by an extrinsic doping with acceptor dopants considering their low E_v when the Fermi level is near the HOMO region. The CI-NEB results showed that the lowest E_m values of the eight Ln₃GaO₆ compounds among various migration paths of V_O²⁺ ranged from 0.59 to 0.79 eV. The E_m values of these Ln₃GaO₆ compounds are much lower than the counterpart Ln₂O₃, which may be ascribed to the shortened distances between O sites by the participation of Ga³⁺, which has smaller ionic radius. Among the Ln₃GaO₆ compounds, La₃GaO₆ has the lowest E_m, followed by Nd₃GaO₆. The FPMD results showed a strong correlation of the lowest E_m obtained by CI-NEB and the D_O at 1873 K of Ln₃GaO₆. Finally, we predicted the σ_O at 1073 K (800 °C) based on the computed D_O and the E_m obtained by CI-NEB and confirmed that the predicted σ_O and the previously reported experimental values have a similar order of magnitude. In particular, La₃GaO₆ shows a higher predicted σ_O of 2.2 × 10⁻² S/cm at 1073 K (800 °C) than those of previously reported Nd₃GaO₆ and Gd₃GaO₆. In addition, these eight Ln₃GaO₆ compounds are expected to have a similar σ_O value as that of the common reference value of 1 × 10⁻² S/cm of YSZ at 1000 K. Considering that the optimum conditions for YSZ, such as the concentration of dopant, and the processing condition have been developed for decades, we expect that these oxides also have potential as suitable oxygen-ion conductors.

The findings of this study can provide insights into the investigation on Ln₃GaO₆-based oxygen-ion conductors. The σ_O values of these oxides can be further improved by additional cation doping and mixing. In addition, Ln₃GaO₆ compounds have lower melting temperatures than YSZ (~1700 °C [49–51,69] vs. >2700 °C [70]), which can provide an advantage of processing at lower temperatures.

Acknowledgments

The authors thank Drs. Shin Tajima, Seiji Kajita, and Hiroshi Ohno in Toyota Central R&D Lab. for fruitful discussions.

Disclosure statement

No potential conflict of interest was reported by the authors.

Supplementary material

Supplemental data for this article can be accessed here.