A data-driven comparative study of thermomechanical properties in rare-earth zirconate and tantalate oxides for thermal barrier coatings

Computational workflow

A high-throughput computational workflow was employed to establish structure–property relationships in RE₃TaO₇ and RE₂Zr₂O₇ oxides, as shown in Figure 1. The workflow begins with compositional enumeration across 17 RE elements (RE = Sc, Y, La ~ Lu), followed by crystal structure selection based on predicted stable phases. Subsequently, density functional theory (DFT) calculations were performed to extract key thermodynamic and bonding descriptors. These data feed into property prediction models for K_IC and κ_L. Finally, data-driven analyses, including principal component analysis (PCA), Pearson correlation, and SHapley Additive exPlanations (SHAP)-based interpretability methods, were applied to systematically evaluate the influence of each descriptor on the predicted K_IC and κ_L. These analyses identify the most relevant features and provide both statistical and model-driven insight.

Figure 1. Schematic overview of the computational workflow for uncovering the structure–property relationships of the RE₃TaO₇ and RE₂Zr₂O₇ oxides.

Crystal structure

The RE₃TaO₇ and RE₂Zr₂O₇ oxide families exhibit a rich diversity of crystal structures, each governed by the ionic radius of the RE ($$ r_{RE^{3+}} $$) cation and the charge and size of the B-site cation (Ta⁵⁺ or Zr⁴⁺). These structures can be categorized into ordered and disordered fluorite-derived frameworks, with representative models illustrated in Figure 1. In the RE₃TaO₇ system, three primary structure types are observed. The defect fluorite structure (space group Fm$$ \bar{3} $$m) features a disordered arrangement of RE³⁺ and Ta⁵⁺ on the cation sublattice, with one-eighth of oxygen sites vacant for charge compensation^[11]. Larger RE³⁺ cations stabilize the Cmcm phase, whereas medium-sized RE³⁺ elements favor the C222₁ structure, commonly referred to as the Weberite type. Both ordered structures exhibit distinct BO₆ octahedral frameworks and layered arrangements of oxygen vacancies, leading to pronounced structural anisotropy. In contrast, RE₂Zr₂O₇ compounds crystallize in either the defect fluorite (Fm$$ \bar{3} $$m) or ordered pyrochlore (Fd$$ \bar{3} $$m) structures, depending on RE³⁺ size^[45]. The pyrochlore phase features an ordered arrangement of RE and Zr on the A and B sublattices, respectively, with one-eighth of the anion sites vacant in a periodic fashion. Smaller RE³⁺ ions favor the defect fluorite phase, which exhibits full cationic disorder and randomly distributed oxygen vacancies. The transformation from ordered to disordered fluorite structures in both oxide families reflects the increasing ability of the lattice to accommodate strain as ionic radii decrease. Specifically, given the lack of experimental evidence for this compound, the results for Sc₂Zr₂O₇ and Sc₃TaO₇ are included for computational completeness only and should be interpreted with caution.

First-principles calculations

All DFT calculations were carried out using the projector augmented wave (PAW) method with the Perdew–Burke–Ernzerhof (PBE) exchange–correlation functional within the generalized gradient approximation (GGA) as implemented in the Vienna Ab initio Simulation Package (VASP)^[46]. Spin polarization was enabled, and the plane-wave energy cutoff was set to 1.4 times the maximum elemental ENMAX, where ENMAX denotes the recommended kinetic energy cutoff specified in the VASP pseudopotential files. k-point meshes of 4 × 4 × 4 ensured convergence of total energy within 10^-6 eV·atom^-1. Although applying Hubbard U corrections may improve prediction accuracy, the primary focus of this work is on trend-level comparisons across a wide compositional space. Several studies have established the effectiveness of GGA functionals in describing the properties of RE compounds^[47-49]. Therefore, omitting U corrections introduces negligible error in the context of structural, elastic, and thermodynamic analysis. Moreover, this choice avoids the element-specific arbitrariness associated with U parameter selection and ensures methodological consistency throughout the high-throughput workflow^[49]. Accordingly, all calculations in this study were carried out using the PBE functional within the GGA framework without +U corrections^[50], providing a robust and computationally efficient platform for capturing structure–property relationships in RE zirconates and tantalates. The V₀ and equilibrium energy (E₀) were obtained by fitting the total energy–volume data to the four-parameter Birch-Murnaghan equation of state^[51]. Elastic constants were obtained using the energy–strain method, from which mechanical properties including bulk modulus, shear modulus, Young’s modulus, and Poisson’s ratio were derived via the Voigt–Reuss–Hill averaging scheme^[52,53].

Lattice thermal conductivity

The thermodynamic combinatorial model is employed to predict the Debye temperature Θ, Grüneisen parameters γ to the temperature-dependent κ_L in oxides^[27]. Specifically,

have been selected as the best descriptors for calculating κ_L in the Slack equation:

where A is a constant, M_av is average atomic mass, n denotes the number of atoms in the unit cell, V₀ represents the atomic volume. This formulation was applied to RE₂Zr₂O₇ and RE₃TaO₇ oxides with temperatures ranging from 200 to 1,800 K, as reported in Ref^[27].

The effective thermal conductivity κ_eff was normalized as

where the porosity ϕ is set to 5% in accordance with experiment^[54].

Fracture toughness

The Griffith Criterion is a foundational theory in fracture mechanics for brittle materials. It postulates that crack propagation occurs when the elastic strain energy released by crack extension is sufficient to overcome the energy required to create the new crack surfaces (surface energy, γ). K_IC is related to the Young’s modulus (E) and fracture energy (G_c = 2γ for mode I fracture in brittle materials) through an expression of the form^[55,56]

where ν is Poisson’s ratio. Calculation details for the surface energy (γ) are provided in Ref^[57]. It is important to note that Griffith theory generally predicts the fracture strength of materials containing pre-existing flaws, rather than the ideal theoretical strength of a perfect crystal. In practice, semi-empirical predictive models, derived from fitting experimental data, are often employed to provide correction factors for different material systems^[58,59].

Structural stability and thermodynamic trends

The thermodynamic stability, V₀, and E₀ of RE₃TaO₇ and RE₂Zr₂O₇ oxides are systematically evaluated across 17 RE elements (RE = Sc, Y, La ~ Lu) using first-principles calculations in the relevant crystal structures as illustrated in Figure 2. For the RE₃TaO₇ oxides, V₀ exhibits a consistent decreasing trend corresponding to the well-documented lanthanide contraction as the ionic radius of RE ions decreases. Among the investigated configurations, the orthorhombic C222₁ phase is thermodynamically most favorable for RE³⁺ cations of intermediate size (e.g., Sm ~ Dy, Y), while the Cmcm structure is stabilized by the largest cations (La, Pr). The disordered fluorite structure becomes energetically competitive for smaller RE ions due to its ability to accommodate local strain through enhanced structural flexibility^[60,61]. Similarly, the RE₂Zr₂O₇ oxides also demonstrate volume contraction consistent with decreasing $$ r_{RE^{3+}} $$. However, these zirconates exhibit notably lower E_form compared to their tantalate counterparts, indicating superior intrinsic thermodynamic stability. In RE₂Zr₂O₇ oxides, an ordered pyrochlore structure (Fd$$ \bar{3} $$m) is the most stable for larger RE ions, while a disordered fluorite structure is increasingly stabilized as the $$ r_{RE^{3+}} $$ decreases, driven primarily by structural flexibility. Notably, the transition to disorder occurs at a larger ionic radius in RE₂Zr₂O₇ than in RE₃TaO₇, underscoring the greater rigidity and lower defect tolerance of the Zr–O framework compared to Ta–O^[19].

Figure 2. Thermodynamic analysis for (A-C) RE₃TaO₇ and (D-F) RE₂Zr₂O₇ comparisons (RE = Sc, Y, La ~ Lu), including (A and D) Equilibrium energy, (B and E) equilibrium volume and (C and F) formation energies of different phases. Symbols denote space groups: C222₁, Cmcm and Fm$$ \bar{3} $$m for RE₃TaO₇; Fd$$ \bar{3} $$m and Fm$$ \bar{3} $$m for RE₂Zr₂O₇ oxides. RE: Rare-earth.

To identify underlying correlations among structural and thermodynamic parameters, PCA was conducted using four standardized descriptors: $$ r_{RE^{3+}} $$, E₀, V₀, and E_form. Given the heterogeneous physical dimensions of these features, z-score normalization was applied, and the correlation matrix was selected as the basis for decomposition. For both RE₃TaO₇ and RE₂Zr₂O₇ systems, the first two principal components account for over 95% of the total variance, indicating that the dimensionality reduction retains the vast majority of structural–thermodynamic information. The loading patterns suggest that E_form and V₀ contribute most strongly to the first principal component (PC1), whereas $$ r_{RE^{3+}} $$ and E₀ dominate the second component (PC2). This orthogonal separation implies that atomic packing (reflected by V₀) and chemical stability (reflected by E_form) capture complementary axes of variability across the RE series.

Thermal conductivity and fracture toughness

The high-temperature performance of RE₂Zr₂O₇ and RE₃TaO₇ oxides as TBCs is strongly governed by their κ_L and K_IC. In this study, κ_L and K_IC were systematically evaluated for the pyrochlore (Fd$$ \bar{3} $$m) RE₂Zr₂O₇ and Weberite-type (C222₁) RE₃TaO₇ phases using a combined model for κ_L and the Griffith-based energetic model for K_IC. The adopted models for predicting κ_L and intrinsic K_IC have been previously validated by prior studies on HEOs^[27,62]. Specifically, the Debye-Slack formalism and its extensions have shown predictive accuracy within 15% when benchmarked against experimental κ_L for RE oxides, particularly those with well-defined phonon spectra and moderate anharmonicity^[27,44]. For K_IC, the Griffith-based energetic models employed here, especially the variant incorporating surface energy and Young’s modulus, are recognized as robust estimators for brittle systems, capturing relative trends and chemical effects despite neglecting extrinsic toughening mechanisms^[62]. These models provide a reliable basis for high-throughput screening, balancing physical realism with computational feasibility.

To ensure the validity of thermal conductivity predictions across the studied temperature range (200-1,800 K), potential structural phase transitions were carefully considered. For RE₂Zr₂O₇, a well-established order–disorder transition from the pyrochlore to defect-fluorite phase occurs at elevated temperatures, typically between 1,530 and 2,400 °C (1,803-2,673 K), depending on the $$ r_{RE^{3+}} $$^[18] Similarly, RE₃TaO₇ oxides are known to undergo a transition from the orthorhombic weberite-type (C222₁) structure to a disordered fluorite phase at temperatures generally above 1,973 K^[10,63]. These transition points lie beyond the temperature range considered in this work. Accordingly, all κ_L calculations were conducted using the most stable DFT-predicted phase at 0 K for each composition. The quasi-harmonic Debye–Grüneisen model employed here assumes a single, temperature-invariant phase and is therefore valid up to the onset of these high-temperature transitions. As no phase changes are expected within the modeled range (200-1,800 K), the resulting predictions reliably capture the intrinsic thermal transport behavior of the respective ground-state structures.

The calculated temperature-dependent κ_L from 200 to 1,800 K, incorporating 5% porosity corrections for representative oxides in both families is shown conceptually in Figure 3A and B. At lower temperatures, κ_L rapidly decreases following a typical inverse temperature dependence (1/T), governed by intrinsic phonon–phonon Umklapp scattering. At higher temperatures, κ_L approaches a minimum limit (κ_min) due to phonon mean free paths converging toward interatomic distances, indicative of pronounced phonon scattering. As shown in Figure 3C, RE₃TaO₇ consistently exhibits lower κ_L (0.76-1.52 W·m^-1·K^-1 at 300 K) compared to RE₂Zr₂O₇ (1.7-2.8 W·m^-1·K^-1 at 300 K), a trend maintained across the entire temperature range. The lower κ_L of RE₃TaO₇ oxides primarily results from their inherently lower symmetry (C222₁), greater structural complexity and consequently increased phonon scattering. Additional factors, including larger atomic mass differences and oxygen sublattice disorder, further reduce κ_L by intensifying phonon scattering channels.

Figure 3. Temperature-dependent κ_L for (A) C222₁ RETa₃O₇ and (B) RE₂Zr₂O₇ pyrochlore oxides by the combinatorial approach (κ₍₂₎(Θ₍₂₎(γ⁽³⁾))) with 5% porosity correction; (C) κ_L values at 300 and 1,473 K. Calculated K_IC by K₀⁽ⁿ⁾ where n = 1-4 of for (D) RE₃TaO₇ and (E) RE₂Zr₂O₇ oxides, together with the reported result^[69-80]; (F) Comparison of the minimum and maximum predicted K_IC values by K₀⁽¹⁾.

For the identification of the optimal model for predicting K_IC, Figure 3D-F compares the results obtained from various models (K₀⁽¹⁾ and K₀⁽²⁾ for energetic, K₀⁽³⁾ and K₀⁽⁴⁾ for empirical) for oxides, alongside the reported results. Among these models, K₀⁽¹⁾ = $$ \sqrt{\frac{2E\gamma }{1-\nu ^2}} $$^[55] and K₀⁽²⁾ = $$ \sqrt{2E\gamma} $$^[56] represent toughness estimates particularly suitable for brittle fracture systems. K₀⁽³⁾ = (1 + α)V₀^1/6G(B/G)^1/2[58] (α = 0) and K₀⁽⁴⁾ = α₀^-1/2V₀^1/6(ζ(v)E)^3/2[59] (α₀ = 8,840 GPa), are empirical models that consider modulus ratios and bond characteristics and provide complementary insight into ionic and covalent contributions. The shaded regions represent K_IC limits calculated using the highest and lowest surface energy values for each structure. Through plotting surface energy (γ) as a function of relative orientation along planes [i.e., from (001) → (012) → (011) → (021) → (010)], the C222₁ and pyrochlore structures have the lowest γ with the most energetic favorable (031) and (110) surfaces, respectively. As shown in Figure 3D and E, K_IC increases with decreasing $$ r_{RE^{3+}} $$ in both systems. Among the tested models, K₀⁽¹⁾ was selected for subsequent analysis due to its consistent performance and theoretical robustness. Intriguingly, the calculations predict that RE₂Zr₂O₇ generally exhibits slightly higher intrinsic K_IC (~1.80 MPa·m^1/2) than RE₃TaO₇ (~1.6 MPa·m^1/2), particularly for mid-sized RE ions. This is attributed to the highly ordered pyrochlore structure in RE₂Zr₂O₇, where periodic oxygen vacancy arrangements may enable crack deflection or microcrack formation, where mechanisms favorable for intrinsic toughening. However, this theoretical prediction contrasts with experimental observations, where RE₃TaO₇ often demonstrates higher or comparable K_IC (1.0-2.6 MPa·m^1/2)^[11] than RE₂Zr₂O₇. This discrepancy arises from the difference between intrinsic (defect-free) and extrinsic (microstructure-sensitive) K_IC. Theoretical models inherently exclude microstructural effects such as grain boundaries, porosity, domain switching and residual stress, which play a crucial role in experimental measurements.

Figure 3F summarizes these trends, highlighting that while RE₂Zr₂O₇ oxides demonstrate uniformly high intrinsic K_IC across the RE series, the C222₁-type RE₃TaO₇ oxides exhibit greater variability, with maximum toughness predicted for Tm₃TaO₇ and Lu₃TaO₇. The superior intrinsic performance of RE₂Zr₂O₇ is closely linked to its symmetric bonding network and ordered defect structure, which provide effective energy dissipation near crack tips^[64]. In contrast, the more complex and disordered lattice of RE₃TaO₇ limits such toughening routes under idealized conditions. Moreover, the combined plot of K_IC vs. κ_L at room temperature directly illustrates the inherent trade-off between mechanical robustness and thermal insulation, as shown in Figure 4. RE₂Zr₂O₇ oxides cluster within a region of high K_IC and moderate thermal conductivity, ideal for applications prioritizing mechanical integrity. In contrast, the orthorhombic C222₁-type RE₃TaO₇ oxides occupy the low-κ_L, moderate-K_IC region, which is advantageous for maximizing thermal insulation, albeit potentially at the expense of mechanical robustness. These results, validated against available experimental data, underscore the effectiveness of the data-driven predictive framework employed.

Figure 4. The relationship between lattice thermal conductivity at 300 K by and fracture toughness for RE₃TaO₇ and RE₂Zr₂O₇ with reported data^{[10,76,81-83]}.

Atomic and electronic attributes of lattice distortion

To clarify the atomistic origins governing the thermal conductivity and K_IC of RE-based RE₃TaO₇ and RE₂Zr₂O₇ oxides, detailed analyses of bonding environments were conducted. These include bonding charge density distributions and quantitative descriptors related to bond length and charge heterogeneity. Such atomic-scale characterizations directly inform the macroscopic property trends observed, offering a unified descriptor-based framework justified by the structural similarities between the two oxide families - both composed of interconnected BO₆ octahedra (B = Ta or Zr) and distorted RE–O polyhedra, despite differences in nominal RE:O ratios.

As shown in Figure 5, bonding charge density (Δρ)^[65] isosurfaces reveal marked differences in electronic distribution between the orthorhombic C222₁-type RE₃TaO₇ and pyrochlore RE₂Zr₂O₇ oxides. In both cases, yellow and blue regions around RE and O atoms indicate strong ionic character. However, the RE₃TaO₇ oxides exhibit pronounced anisotropy in charge distribution, consistent with their lower crystal symmetry and indicative of stronger local lattice distortions. In contrast, the zirconates show more isotropic and symmetric Δρ patterns, aligned with their higher lattice symmetry and more ordered vacancy configuration.

Figure 5. The bonding charge density isosurface of (A) C222₁ RETa₃O₇ and (B) RE₂Zr₂O₇ pyrochlore (RE = Sc, Y, La ~ Lu). The blue and yellow isosurfaces respectively correspond to Δρ ± 0.02 e^-Å^-3. RE: Rare-earth.

To quantitatively assess local structural disorder, two descriptors were defined: bond-length variance (D_b) and Bader charge variance (D_e):

where b_i and e_i represent individual bond lengths and Bader charges within RE–O or Zr/Ta–O coordination environments, and $$ \bar{b} $$ and $$ \bar{e} $$ are their respective averages. N is the number of bonds or charges considered. These disorder metrics quantify fluctuations in bond distance and electronic environment, which have been linked to enhanced phonon scattering and reduced κ_L in complex oxides^[14,30].

As shown in Figure 6A, the relative bond lengths across different RE–O and Zr/Ta–O coordination units (REO₇, REO₈, and mixed REO_7/8 polyhedra) decrease with decreasing $$ r_{RE^{3+}} $$, reflecting the well-known lanthanide contraction. The insets illustrate representative local coordination geometries for these polyhedral units. Figure 6B presents the corresponding bond-length disorder degree D_b, which varies notably across both RE species and structural types. Materials with broader bond-length distributions, such as C222₁-type RE₃TaO₇ oxides, exhibit stronger phonon scattering and hence lower κ_L, especially at elevated temperatures^[18]. Within RE₃TaO₇ oxides, smaller RE ions tend to form shorter, stronger RE–O bonds, but the simultaneous increase in D_b moderates thermal transport by amplifying disorder-induced scattering. However, the distribution of bond lengths around the average provides a crucial descriptor related to local structural heterogeneity and strain fields within the lattice. A measure of this heterogeneity, such as the distortion degree of RE–O bond lengths, also varies across the RE series and between structures, as shown in Figure 6B. This bond length heterogeneity emerges as a significant factor influencing thermal conductivity. Materials exhibiting broader bond length distributions (higher distortion) tend to have lower κ_L particularly at high temperatures, due to enhanced phonon scattering from local disorder^[18]. For instance, the generally lower κ_L of RE₃TaO₇ oxides compared to zirconates can be partly attributed to the greater structural complexity and potential for bond length variation in the lower-symmetry C222₁ structure. Within the RE₃TaO₇ oxides, smaller RE ions lead to shorter, stronger bonds, but the interplay with bond length heterogeneity dictates the overall κ_L. Additionally, K_IC is also intrinsically linked to the bonding environment. The trend of increasing K_IC with decreasing $$ r_{RE^{3+}} $$ in RE₃TaO₇ oxides correlates with the formation of shorter, stronger average RE–O bonds, which increases the intrinsic resistance to bond breaking at the crack tip^[66,67], as shown in Figure 4. The generally higher toughness of pyrochlore zirconates suggests their specific ordered structure provides more effective toughening mechanisms compared to the C222₁ RE₃TaO₇ oxides.

Figure 6. Bonding and charge disorder analysis in RE₃TaO₇ and RE₂Zr₂O₇ oxides. (A) Relative bond lengths for RE–O (left) and Zr/Ta–O (right) polyhedra, where “REO₇”, “REO₈” and “REO_7/8” denote the 7- and 8-fold and mixed-coordination RE–O environments. Insets illustrate representative polyhedral units; (B) Bond-length disorder degree D_b of RE–O (left) and Zr/Ta–O (right) polyhedra, using the same symbol scheme as in (A); (C) Valence charges of cations RE³⁺ (left) and anions O^2- (right), derived from Bader charge analysis; (D) Valence charge disorder degree D_e for RE³⁺ (D_e^RE) and O^2- (D_e^O), quantifying charge inhomogeneity within coordination polyhedra. RE: Rare-earth.

In Figure 6C and D, electronic structure insights from Bader charge analysis further clarify structure–property relationships. Most RE ions exhibit Bader charges consistent with a nominal +3 oxidation state. However, Eu and Yb in the RE₃TaO₇ system show significantly lower calculated charges, indicative of a preference for the divalent (+2) state. As shown in Figure 6D, this deviation introduces local charge imbalance and strain fields, broadening the bond-length distribution and increasing D_e, which in turn destabilizes the ordered C222₁ structure (reflected in higher relative E_form, as shown in Figure 2). Interestingly, such pronounced effects are not observed in the corresponding zirconates, suggesting the more flexible fluorite-derived framework of the zirconates is better able to accommodate such charge and strain fluctuations. These mechanisms provide vital guidelines for strategically balancing K_IC and thermal insulation in RE-based oxide ceramics for advanced TBC applications.

Quantifying structure–property relationships

To establish direct links between atomic-scale descriptors and macroscopic thermomechanical behavior, quantitative correlation and interpretability analyses were conducted. Figure 7 presents Pearson correlation coefficient matrices and SHAP feature importance rankings for both RE₃TaO₇ and RE₂Zr₂O₇ oxides. The SHAP package was utilized to interpret the ML models effectively^[68]. These tools collectively identify the most influential parameters controlling κ_L and K_IC, including the $$ r_{RE^{3+}} $$, V₀, E₀, E_form, bond energy (E_b) for RE–O, Ta–O and Zr–O bonds, and disorder metrics such as bond-length heterogeneity (D_b) and Bader charge disorder (D_e). Hierarchical clustering within the correlation matrices further highlights groups of interrelated descriptors, as visualized through dendrograms and colored sidebars.

Figure 7. Correlation analysis across RE₃TaO₇ and RE₂Zr₂O₇ between target properties and fundamental descriptors, including are-earth ionic radius ($$ r_{RE^{3+}} $$), equilibrium volume (V₀), total energy (E₀), formation energy (E_form), bond energy (E_b) for RE–O, Ta–O, and Zr–O bonds, bond-length disorder (D_b) and Bader charge disorder (D_e). Pearson correlation heatmap for (A) RE₃TaO₇ and (B) RE₂Zr₂O₇. Colored sidebars indicate feature clusters derived from hierarchical clustering; (C) SHAP analysis showing the top five most influential features (ranked by mean SHAP value) for predicting K_IC and κ_L in each system. RE: Rare-earth; SHAP: SHapley Additive exPlanations.

In RE₃TaO₇ oxides, as shown in Figure 7A, strong positive correlations are observed between $$ r_{RE^{3+}} $$, E₀ and E_b^RE–O with both K_IC and κ_L. This indicates that larger RE ions and stronger RE–O bonding contribute to enhanced structural rigidity and thermal transport capacity. Conversely, descriptors capturing local lattice disorder, such as bond-length heterogeneity (D_b^REO7/8) and RE-site Bader charge disorder (D_e^RE), also show positive correlation with K_IC, suggesting that structural distortion aids crack resistance by disrupting crack propagation pathways. Interestingly, Ta–O E_b (E_b^Ta–O) and its associated disorder (D_b^TaO6) display inverse correlations with κ_L, indicating that stronger, more heterogeneous Ta–O interactions exacerbate phonon scattering, thereby suppressing thermal conductivity. These findings reinforce the dual role of bonding strength and disorder in tailoring functional performance. As shown in Figure 7B, the RE₂Zr₂O₇ pyrochlore oxides exhibit similar, though somewhat attenuated, correlation patterns. This attenuation is attributed to their higher symmetry and reduced chemical complexity. Nevertheless, key descriptors such as Zr–O E_b (E_b^Zr–O), RE Bader charge disorder (D_e^RE) and RE–O bond-length disorder (D_b^REO8) remain strongly correlated with both K_IC and κ_L, confirming their general importance across systems.

To further quantify the relative importance of these descriptors, SHAP analysis was applied to trained regression models for each oxide family, as shown in Figure 7C. For RE₃TaO₇, dominant predictive features include RE–O E_b, bond-length disorder (D_b), and E_form - consistent with the physical picture of configurational complexity and bond robustness as key regulators of performance. In RE₂Zr₂O₇, the most influential features shift slightly to Zr–O E_b, Bader charge disorder, and V₀, emphasizing the role of atomic packing and local charge environment. Taken together, both correlation-based and model-based analyses converge on three dominant descriptors: E_b, D_b and Bader charge heterogeneity (D_e), as the primary determinants of K_IC and κ_L across the RE series. These interpretable, physically grounded parameters elucidate the mechanisms of phonon scattering, lattice stiffening and crack resistance. Moreover, they serve as a robust foundation for data-driven materials design. Importantly, these insights explain why RE₃TaO₇ oxides, despite exhibiting intrinsically lower thermal conductivity, maintain moderate toughness. The reinforcement of RE–O bonds and local disorder-induced energy dissipation mechanisms helps compensate for structural complexity. In contrast, RE₂Zr₂O₇ oxides achieve a superior balance of high fracture resistance and moderate κ_L, attributable to their chemically ordered and structurally accommodating lattice framework.