Effect of Rare Earth Elements on Stability and Sintering Resistance of Tetragonal Zirconia for Advanced Thermal Barrier Coatings

: The effect of dopant species on the sintering resistance of zirconia-based ceramics remains a huge challenge in terms of both experiment and theory. As one of the most popular materials for high-temperature protective coatings, it is still urgent to obtain rare earth-doped ZrO 2 with high sintering resistance and good phase stability. Here, the sintering resistance and phase stability of rare earth oxides (La 2 O 3 , Nd 2 O 3 , Gd 2 O 3 , and Y 2 O 3 )-stabilized zirconia (ZrO 2 ) were thoroughly studied by theoretical and experimental methods. According to experimental data, ZrO 2 doped with rare earth ions with larger radii (La 3+ , Nd 3+ , and Gd 3+ ) exhibited improved sintering resistance at reduced tetragonal phase stability. Molecular dynamics simulation results revealed that rare earth ions with larger ionic radii are prone to segregation at grain boundaries, which can more effectively reduce the grain boundary energy in the materials under consideration. Therefore, the proposed approach involving doping of NdO 1 . 5 (~1 mol%) and YO 1 . 5 (YbO 1 . 5 , 6 mol%) in ZrO 2 is considered to be a promising route for the effective preparation of sinter-resistant ZrO 2 -based ceramics for refractory and thermal barrier materials.


Introduction
Thermal barrier coatings (TBCs) are advanced coatings applied on metallic surfaces in exhaust heat management systems such as gas turbines or aero-engine parts, which are continuously exposed to high temperatures [1][2][3][4]. For their direct utilization as thermal insulation materials, TBCs with thickness in the range of~150-500 µm are usually deposited onto superalloys to protect them from high-temperature erosion [1]. Traditional TBCs are typically represented by 7-8 wt.% yttria-stabilized zirconia (YSZ), which exhibits excellent thermal and mechanical properties, such as low thermal conductivity and high fracture toughness [2]. Another important characteristic of YSZ-based TBCs subjected to high temperatures is the sintering resistance. When the service temperature of YSZ exceeds 1100 • C, pores inside the YSZ coatings undergo shrinkage, particularly those perpendicular to the heat flow, thus leading to a significant increase in the thermal conductivity of TBCs. Li et al. [5][6][7][8]. systematically investigated the sintering processes of YSZ TBCs, discovering that the closure of two-dimensional (2D) pores noticeably increased the thermal conductivity by nearly 50% after 20 h of calcination at 1200 • C and reduced the mechanical properties of the corresponding coatings. For example, their elastic moduli increased drastically from 70 to 150 GPa, reducing the thermal shock resistance and even leading to premature shedding of the coatings [5].
Sintering temperature and sintering time are among the most essential sintering parameters for YSZ coatings, and the increase of these parameters leads to the increase in the sintering rate of YSZ coatings [7,9,10]. Sintering of YSZ coatings takes place in three stages, i.e., the early, the medium, and the final sintering stages, which directly impact the thermal and mechanical properties, leading to the preliminary failure of YSZ TBCs during operation [7]. Therefore, the sintering activity of YSZ coatings should be suppressed to avoid the rapid shrinkage of pores during the early sintering stage. Furthermore, Pan et al. [11] found that certain impurity contents in YSZ coatings could also accelerate or slow down the sintering rate of the corresponding coatings. For example, a small amount of silica (SiO 2 , 1-2 mol%) made the YSZ grain boundaries cleaner, thus simplifying the atomic diffusion. Moreover, during sintering, a part of SiO 2 was transformed into a liquid phase, which could serve as the medium of mass transfer, thereby accelerating the sintering rate of YSZ TBCs. Jiang et. al. [12] studied the size effect of rare earth ions on the phase stability of tetragonal zirconia (ZrO 2 ) and found that the substitution of Y 3+ with Yb 3+ resulted in the increase in phase stability by decreasing its sintering resistance.
In this context, the development of efficient methods to enhance the anti-sintering ability of YSZ-based TBCs is an urgent task. According to some previous literature reports, two approaches are currently available that aim to improve the anti-sintering ability of YSZ coatings. The first approach involves the introduction of unmolten particles into YSZ coatings, followed by their shrinkage by exposing coatings to high-temperature conditions. This enables the improvement of the sintering resistance, while maintaining low thermal conductivity and elastic moduli. However, the temperatures required for this process are too high (about 1400 • C), inducing the phase transformation of YSZ, and thus reducing the mechanical properties and, consequently, the thermal cycle life of the coatings [13]. The second route involves the variation of the doping elements. For instance, Samth et al. [14] systematically investigated the influence of different dopants on the sintering resistance of YSZ coatings, revealing that dopants with larger ionic radii could improve the sintering resistance and, consequently, the anti-sintering ability of YSZ coatings. Furthermore, even a small amount of dopants (i.e., nearly or slightly less than 4 mol%) was shown to exert a positive influence on the sintering resistance of the coatings by reducing their surface energy via element segregation. Moreover, the total content of rare earth elements (except Sc 3+ ) should not exceed 5 mol% to prevent the formation of a large number of cubic phases and, consequently, maintain a minimum fracture toughness of the coatings [15,16].
Besides, as-sprayed TBCs usually possess a large quantity of 2D/3D pores and some amorphous phases, leading to a relatively high total free energy of the coatings. Therefore, when they are exposed to high temperatures, some phenomena, such as grain growth, pore shrinkage, and pore closure, occur spontaneously owing to the mass transfer [5,7,17]. As a result, the sintering process of TBCs should be fundamentally attributed to a decrease of the systematic free energy when using dopants. The doping elements that are used to stabilize tetragonal ZrO 2 are limited; therefore, rare earth elements are considered as suitable candidates for this purpose. As aforementioned, the substitution of Y with Yb leads to the increase in the phase stability; however, the corresponding sintering resistance is reduced [12]. Therefore, it is still urgent to obtain rare earth-doped ZrO 2 with high sintering resistance and good phase stability. The main objective of this study was to systematically investigate the influence of La 2 O 3 , Nd 2 O 3 , Gd 2 O 3 , and Y 2 O 3 dopants, respectively, on the phase stability and sintering resistance of ZrO 2 ceramics via experimental and theoretical methods. The main contribution and significance of this paper can be further expressed as the following three aspects. Firstly, the sintering resistance of rare earth-doped zirconia was investigated experimentally. The sintering resistance of zirconia-based ceramics could be improved by doping~1 mol% RE 3+ with larger ionic radius to decrease the surface energy. Secondly, according to the molecular dynamics simulations, the dihedral angle between grain boundary and surface was increased after doping RE 3+ with a larger ionic radius, which improved the high-temperature stability inside pores. Thirdly, the decomposition of rare earth-doped zirconia was analyzed experimentally and theoretically. We proposed a theory to explain the decomposition of tetragonal-ZrO 2 in which the formation of defect cluster was the main reason for the phase transformation of RESZ from the tetragonal phase to the cubic and monoclinic phase. According to lattice dynamics simulation, the formation energy of RE 3+ -Vo-RE 3+ increased with the increase of rare earth ionic size, which decreased the phase stability of RESZ.

Experimental
Ln 2 O 3 (Ln = La, Nd, Gd, and Y)-stabilized ZrO 2 (LnSZ) powders were prepared via a sol-gel method [18,19]. First, the appropriate amounts of Ln(NO 3 ) 3 ·6H 2 O and ZrOCl·8H 2 O precursors were taken according to a stoichiometric ratio (Ln 0 . 09 Zr 0 . 91 O 1 . 955 ) and then dissolved in distilled water in a beaker. The solution was stirred for 2 h on a magnetic stirrer at a temperature of 80 • C, then a suitable amount of citric acid was added (so that the ratio of metal ions to citric acid was 1:1), and finally stirred again until it became saline. Further, the obtained salt was subjected to heating at 80 • C for 12 h to achieve the dry gel consistency. The former was ground in a mortar and exposed to 6 h of heat treatment at 1100 • C, allowing the production of the corresponding LnSZ powders. The samples were finally sintered for 90 h at 1300 • C. Moreover, Gd 2 O 3 (2.5%) and Y 2 O 3 (2%) co-doped ZrO 2 was also prepared to study the co-doping effect on the sintering behavior of corresponding materials.
The phase characterization of the obtained LnSZs samples was performed via X-ray diffraction (XRD). The X-ray profiles were recorded at room temperature using a Bruker D8 Advance (Bruker, Massachusetts, United States) diffractometer within a 2θ range of 20-80 • at a scanning rate of 10 • /min. The morphology of the samples was characterized by scanning electron microscopy (SEM) using a JSM-7000F (JEOL, Japan) microscopic system. The average crystalline size of prepared powders was calculated by Scherrer's equation [7]: where λ is the X-ray wavelength and K is a constant that is set to 0.89. θ is the diffraction angle and β is the correct full-width at half maximum (FWHM). To adjust the instrument error, the Gaussian-Gaussian relationship is employed: where B is FWHM and b is the standard width of a reference Al 2 O 3 sample.

Numerical Approach
Molecular dynamics simulation was employed to investigate the segregation behavior of doping elements in ZrO 2 . The segregation ability of rare earth elements was calculated using the LAMMPS open source code [20] and the GULP code [21] combined with a classical short-range Born-Mayer-Buckingham (BMB) potential and a long-range Coulomb potential, and the expressions are represented as follows [22]: where V(r) and r ij are the potential energy of ions i and j, respectively; and A, C, and ρ are the potential parameters listed in Table 1 [22,23]. Table 1. Potential parameters of A, C, and ρ in Equation (3).

Species Pair
A ij (eV) ρ ij ( lomb potential, and the expressions are represented as follows [22]: where V(r) and rij are the potential energy of ions i and j, respectively; and A, C, and ρ are the potential parameters listed in Table 1 [22,23]. 25 During simulation, the ions were assumed to be rigid spheres. Figure 1 depicts two cubic ZrO2 models with two grain boundaries and two crystal surfaces, respectively, simulated in obedience to the periodical boundary conditions. The grain boundary model consisted of 5663 O2 ions, 2830 Zr4+ ions, and 2 Ln3+ (Ln = La, Nd, Sm, and Gd) ions. The two rare earth ions were first inserted far away from the grain boundary (about 2.5 nm) to isolate the interaction between Ln ions and the grain boundary (the cutoff distance of the BMB potential was 10 Å). Comparatively, the two Ln ions were inserted in the grain boundary to create dopant segregation in the GB model. Similarly, the surface model consisted of 4607 O2 ions, 2302 Zr4+, and 2 Ln3+ ions. The modeling of the grain boundary was previously described in detail [22]. The conjugate gradient (CG) approach was used to calculate the cohesive energies of the simulated cubic ZrO2 structures. The energies of the grain boundary and the crystal surface, induced by the segregation of rare earth dopants, could be obtained by subtracting the energy of a perfect cubic ZrO2 structure with the same number of ions, as follows:  C ij (eV lomb potential, and the expressions are represented as follows [22]: where V(r) and rij are the potential energy of ions i and j, respectively; the potential parameters listed in Table 1 [22,23]. During simulation, the ions were assumed to be rigid spheres. F cubic ZrO2 models with two grain boundaries and two crystal surfaces ulated in obedience to the periodical boundary conditions. The grai consisted of 5663 O2 ions, 2830 Zr4+ ions, and 2 Ln3+ (Ln = La, Nd, Sm two rare earth ions were first inserted far away from the grain bound to isolate the interaction between Ln ions and the grain boundary (th the BMB potential was 10 Å). Comparatively, the two Ln ions were in boundary to create dopant segregation in the GB model. Similarly, the sisted of 4607 O2 ions, 2302 Zr4+, and 2 Ln3+ ions. The modeling of t was previously described in detail [22]. The conjugate gradient (CG) a to calculate the cohesive energies of the simulated cubic ZrO2 structur the grain boundary and the crystal surface, induced by the segregatio pants, could be obtained by subtracting the energy of a perfect cubic Z the same number of ions, as follows:  During simulation, the ions were assumed to be rigid spheres. Figure 1 depicts two cubic ZrO 2 models with two grain boundaries and two crystal surfaces, respectively, simulated in obedience to the periodical boundary conditions. The grain boundary model consisted of 5663 O 2 ions, 2830 Zr 4+ ions, and 2 Ln 3+ (Ln = La, Nd, Sm, and Gd) ions. The two rare earth ions were first inserted far away from the grain boundary (about 2.5 nm) to isolate the interaction between Ln ions and the grain boundary (the cutoff distance of the BMB potential was 10 Å). Comparatively, the two Ln ions were inserted in the grain boundary to create dopant segregation in the GB model. Similarly, the surface model consisted of 4607 O 2 ions, 2302 Zr 4+ , and 2 Ln 3+ ions. The modeling of the grain boundary was previously described in detail [22]. The conjugate gradient (CG) approach was used to calculate the cohesive energies of the simulated cubic ZrO 2 structures. The energies of the grain boundary and the crystal surface, induced by the segregation of rare earth dopants, could be obtained by subtracting the energy of a perfect cubic ZrO 2 structure with the same number of ions, as follows: where γ surf and γ GB are the crystal surface energy and the grain boundary energy, respectively.
Crystals 2021, 11, x FOR PEER REVIEW 4 of 1 a classical short-range Born-Mayer-Buckingham (BMB) potential and a long-range Cou lomb potential, and the expressions are represented as follows [22]: where V(r) and rij are the potential energy of ions i and j, respectively; and A, C, and ρ ar the potential parameters listed in Table 1 [22,23]. Table 1. Potential parameters of A, C, and ρ in Equation (3). During simulation, the ions were assumed to be rigid spheres. Figure 1 depicts two cubic ZrO2 models with two grain boundaries and two crystal surfaces, respectively, sim ulated in obedience to the periodical boundary conditions. The grain boundary mode consisted of 5663 O2 ions, 2830 Zr4+ ions, and 2 Ln3+ (Ln = La, Nd, Sm, and Gd) ions. Th two rare earth ions were first inserted far away from the grain boundary (about 2.5 nm to isolate the interaction between Ln ions and the grain boundary (the cutoff distance o the BMB potential was 10 Å). Comparatively, the two Ln ions were inserted in the grain boundary to create dopant segregation in the GB model. Similarly, the surface model con sisted of 4607 O2 ions, 2302 Zr4+, and 2 Ln3+ ions. The modeling of the grain boundary was previously described in detail [22]. The conjugate gradient (CG) approach was used to calculate the cohesive energies of the simulated cubic ZrO2 structures. The energies o the grain boundary and the crystal surface, induced by the segregation of rare earth do pants, could be obtained by subtracting the energy of a perfect cubic ZrO2 structure with the same number of ions, as follows:  Before calculating the energy differences, the stable configuration of dopants and oxygen vacancies was studied. The binding energies of different configurations were calculated and compared. The model used in the present study contained 5 × 5 × 5 ZrO 2 unit cells (1500 atoms) in the initial configuration, two rare earth ions (Y 3+ ) were first placed in two different sites at a distance of >15 Å, and the oxygen vacancy was also created in the O 2 lattice site at a distance of >15 Å from the two rare earth ions. The lattice energy was subsequently calculated through the CG approach. To obtain the optimal configuration, the rare earth ion and oxygen vacancy were made to gradually approach each other to create a series of configurations. In the final step, the binding energy of the Y 3+ -Y 3+ -V O defect cluster was calculated to determine the most stable configuration. The effect of dopant size on the phase stability of rare earth ion-stabilized zirconia (RESZ) was also evaluated by calculating the decomposition energy to form the RE 3+ -V O -RE 3+ defect clusters.

Results and Discussion
3.1. The Anti-Sintering Ability of LnSZ Figure 2 shows the XRD patterns of LnSZ powders calcined for 6 h at 1100 • C, revealing the presence of pyrochlore LZ and tetragonal Ln stabilized ZrO 2 phases. This indicates that Ln elements precipitated and segregated to form the LZ phases after the heat treatment. When the Ln 2 O 3 concentration in LnSZ was below 1.5 mol%, the LnSZ powder transformed from the tetragonal (T) to the monoclinic (M) phase during cooling [24], which seems to disagree with the experimental results shown in Figure 2. In this study, the existence of the tetragonal phase in LaSZ powders is due to the following reasons: (1) the concentration of the tetragonal La 2 O 3 in the LaSZ structure is still higher than 1.5 mol%; and (2) the nanostructure of LaSZ powder results in a more stable tetragonal phase than the cubic phase when the grain size is less than 100 nm [25]. Analysis of the other three powders reveals the presence of only the tetragonal phase, and neither the pyrochlore nor the monoclinic phase is observed in Figure 2. This is a testimony to the lack of segregation of the Ln 2 O 3 particles in the LnSZ powders, making their tetragonal structure more stable.  Figure 3 shows the evolution of the average grain size in LnSZ powders with varying composition. With the increasing radius of rare earth ion, the average grain size of LnSZ powder first decreases and then increases. As a result, NdSZ powders possess the smallest grains. Consequently, the average grain size of YSZ can be effectively reduced by the full or partial substitution of Y2O3 phases in YSZ using rare earth elements with large ionic radii. Figure 3 exhibits a noticeable reduction in the growth rate of grains in the case of the Nd2O3 and La2O3 phases. Therefore, both are expected to improve the sintering resistance of ZrO2. However, owing to the precipitation of the second LZ phase, the La2O3 content in LaSZ decreases, leading to the reduction in the inhibition effect of La2O3 on the ZrO2 grain growth to a certain extent. These experimental results are consistent with the results of some previous reports. For example, Matsumoto et al. found that, compared with Er2O3, Al2O3, or SrO oxides, the addition of 1 mol% La2O3 to a 4YSZ phase could significantly improve the sintering resistance of the YSZ system [14].  Figure 3 shows the evolution of the average grain size in LnSZ powders with varying composition. With the increasing radius of rare earth ion, the average grain size of LnSZ powder first decreases and then increases. As a result, NdSZ powders possess the smallest grains. Consequently, the average grain size of YSZ can be effectively reduced by the full or partial substitution of Y 2 O 3 phases in YSZ using rare earth elements with large ionic radii. Figure 3 exhibits a noticeable reduction in the growth rate of grains in the case of the Nd 2 O 3 and La 2 O 3 phases. Therefore, both are expected to improve the sintering resistance of ZrO 2 . However, owing to the precipitation of the second LZ phase, the La 2 O 3 content in LaSZ decreases, leading to the reduction in the inhibition effect of La 2 O 3 on the ZrO 2 grain growth to a certain extent. These experimental results are consistent with the results of  To investigate the influence of co-doping on the sintering resistance of ZrO2-based ceramics, Gd2O3 and Y2O3 co-doped ZrO2 was also prepared. The most pronounced decrease in the average grain size (by approximately 13%) was achieved by replacing 2 mol% Y2O3 with Gd2O3. This coincides with the results of some previous studies, reflecting the ability to substantially improve the sintering resistance of corresponding materials by totally or partially replacing Y2O3 with Gd2O3. Figure 4 displays the SEM images of GdYSZ samples with different Gd2O3 concentrations, revealing that the average grain sizes in the powders are below 100 nm, which is in good agreement with the above-mentioned XRD results. The smallest average grain size is observed in the GdSZ powder, whereas the largest one refers to the YSZ powder. This result indicates that the YSZ powder is much denser than the Gd2O3-doped YSZ powder, with an obviously smaller number of micro-pores in its structure. Therefore, the addition of Gd2O3 effectively inhibits the mass transfer of ZrO2 and improves the stability of the internal holes, as shown in Figure 4. To investigate the influence of co-doping on the sintering resistance of ZrO 2 -based ceramics, Gd 2 O 3 and Y 2 O 3 co-doped ZrO 2 was also prepared. The most pronounced decrease in the average grain size (by approximately 13%) was achieved by replacing 2 mol% Y 2 O 3 with Gd 2 O 3 . This coincides with the results of some previous studies, reflecting the ability to substantially improve the sintering resistance of corresponding materials by totally or partially replacing Y 2 O 3 with Gd 2 O 3 . Figure 4 displays the SEM images of GdYSZ samples with different Gd 2 O 3 concentrations, revealing that the average grain sizes in the powders are below 100 nm, which is in good agreement with the above-mentioned XRD results. The smallest average grain size is observed in the GdSZ powder, whereas the largest one refers to the YSZ powder. This result indicates that the YSZ powder is much denser than the Gd 2 O 3 -doped YSZ powder, with an obviously smaller number of micro-pores in its structure. Therefore, the addition of Gd 2 O 3 effectively inhibits the mass transfer of ZrO 2 and improves the stability of the internal holes, as shown in Figure 4.

High-Temperature Phase Stability of LnSZ
An improvement in the sintering resistance should guarantee the high-temperature stability of YSZ. According to data reported in the literature [26,27], a 4YSZ system (4 mol% Y 2 O 3 partially stabilized zirconia) possesses a metastable tetragonal form, wherein the inner Y 2 O 3 phase is prone to segregation and precipitation during the heat treatment, resulting in the inhomogeneous distribution of Y 2 O 3 across the YSZ structure [2,28,29]. In turn, the Y 2 O 3 -enriched regions undergo transformation from the tetragonal to the cubic phase, whereas those with Y 2 O 3 deficiency maintain the tetragonal structure. Moreover, further cooling leads to the transformation of the Y 2 O 3 -depleted regions from the tetragonal to the monoclinic phase, which is also accompanied by the increase in volume (about 3.5%) [29]. The phase stability of YSZ is also influenced by external stress, i.e., the application of a sufficiently high tensile stress (~0.5 GPa) induces the transformation of the corresponding phase into a monoclinic state [29]. In a word, the premature phase transformation of LnSZ should be prevented to suppress the formation of a large number of internal cracks. Crystals 2021, 11, x FOR PEER REVIEW 8 of 15

High-Temperature Phase Stability of LnSZ
An improvement in the sintering resistance should guarantee the high-temperature stability of YSZ. According to data reported in the literature [26,27], a 4YSZ system (4 mol% Y2O3 partially stabilized zirconia) possesses a metastable tetragonal form, wherein the inner Y2O3 phase is prone to segregation and precipitation during the heat treatment, resulting in the inhomogeneous distribution of Y2O3 across the YSZ structure [2,28,29]. In turn, the Y2O3-enriched regions undergo transformation from the tetragonal to the cubic phase, whereas those with Y2O3 deficiency maintain the tetragonal structure. Moreover, further cooling leads to the transformation of the Y2O3-depleted regions from the tetragonal to the monoclinic phase, which is also accompanied by the increase in volume (about  Figure 5 shows the XRD patterns of LnSZ powders exposed to 90 h heat treatment at 1300 • C. Figure 5 demonstrates that the full width at a half-maximum (FWHM) of the diffraction peak is significantly decreased, indicating the noticeable grain growth after high-temperature annealing. Moreover, a large fraction of monoclinic ZrO 2 was detected in LaSZ and NdSZ. Furthermore, no traces of the monoclinic phase were highlighted in GdSZ, GdYSZ, and YSZ powders after the long-term exposure at 1300 • C. The cubic phase along with some amount of the tetragonal phase was present in all the powders, except for LaSZ. This indicates that the tetragonal LaSZ decomposed into the cubic, monoclinic, and pyrochlore phases, respectively, which could be attributed to the precipitation of La 2 O 3 in LaSZ. As a result, M-ZrO 2 was found to be more stable than T-ZrO 2 when the content of La 2 O 3 was less than 1.5 mol% with a large grain size (large than 1 µm). This is also applicable to NdSZ. However, with the increasing Nd content, Nd 2 O 3 -stabilized ZrO 2 with a cubic structure is expected to be formed, rather than Nd 2 Zr 2 O 7 (NZ) with a pyrochlore structure. Furthermore, a serious segregation of Gd also occurred in GdSZ, which was caused by the decomposition of GdSZ into the cubic and tetragonal Gd 2 O 3 -depleted phases. Figure 5b presents that GbSZ possesses two pronounced peaks corresponding to (111) and (110) planes, attributed to the cubic and tetragonal phases, respectively. A similar decomposition of the tetragonal phase also emerged in GdYSZ and YSZ, leading to the formation of the cubic phase. As a result, the phase stability of LnSZ increased with the decrease of the ionic radius of the dopant.
in LaSZ and NdSZ. Furthermore, no traces of the monoclinic phase were highlighted in GdSZ, GdYSZ, and YSZ powders after the long-term exposure at 1300 °C. The cubic phase along with some amount of the tetragonal phase was present in all the powders, except for LaSZ. This indicates that the tetragonal LaSZ decomposed into the cubic, monoclinic, and pyrochlore phases, respectively, which could be attributed to the precipitation of La2O3 in LaSZ. As a result, M-ZrO2 was found to be more stable than T-ZrO2 when the content of La2O3 was less than 1.5 mol% with a large grain size (large than 1 μm). This is also applicable to NdSZ. However, with the increasing Nd content, Nd2O3-stabilized ZrO2 with a cubic structure is expected to be formed, rather than Nd2Zr2O7 (NZ) with a pyrochlore structure. Furthermore, a serious segregation of Gd also occurred in GdSZ, which was caused by the decomposition of GdSZ into the cubic and tetragonal Gd2O3-depleted phases. Figure 5(b) presents that GbSZ possesses two pronounced peaks corresponding to (111) and (110) planes, attributed to the cubic and tetragonal phases, respectively. A similar decomposition of the tetragonal phase also emerged in GdYSZ and YSZ, leading to the formation of the cubic phase. As a result, the phase stability of LnSZ increased with the decrease of the ionic radius of the dopant. The aforementioned results indicate that, even though the stabilizers with larger ionic radii allow for the effective improvement of the sintering resistance of ZrO2, these rare earth ions, while precipitated, could also be formed more easily than the cubic or pyrochlore phase, which significantly reduced the phase stability of the corresponding materials. Therefore, the excellent phase stability of TBCs was ensured by using ZrO2 with large-sized ions under the operating conditions at moderate and/or low temperatures (e.g., 500-700 ℃ in the case of a diesel engine), while the amount of rare earth elements with a large radii should be less than 1 mol% at higher temperatures (about 1000 ℃). Thus, The aforementioned results indicate that, even though the stabilizers with larger ionic radii allow for the effective improvement of the sintering resistance of ZrO 2 , these rare earth ions, while precipitated, could also be formed more easily than the cubic or pyrochlore phase, which significantly reduced the phase stability of the corresponding materials. Therefore, the excellent phase stability of TBCs was ensured by using ZrO 2 with large-sized ions under the operating conditions at moderate and/or low temperatures (e.g., 500-700 • C in the case of a diesel engine), while the amount of rare earth elements with a large radii should be less than 1 mol% at higher temperatures (about 1000 • C). Thus, the optimum strategy involved the utilization of small-sized rare earth ions (such as Y 3+ , Er 3+ , or Yb 3+ ) at small doping concentrations (0.5-1.5 mol%). As a result, the improvement of sintering resistance of the corresponding materials was achieved through the reduction of the surface energy and the grain boundary energy of ZrO 2 along with a decrease in the driving force for the sintering reaction. Figure 6 displays the backscattered electron image of the LaSZ sample exposed to 90 h heat treatment at 1300 • C. The white particles highlighted in this image are attributed to La-enriched grains, indicating the precipitation of the LZ phase at a high (above 1.5 mol%) concentration of La 2 O 3 in LaSZ. Besides, numerous micro-cracks between the neighboring grains are observed in Figure 6, whose emergence is due to the volume change of ZrO 2 transforming from the tetragonal to the monoclinic phase during cooling. Figure 6 displays the backscattered electron image of the LaSZ sample exposed to 90 h heat treatment at 1300 ℃. The white particles highlighted in this image are attributed to La-enriched grains, indicating the precipitation of the LZ phase at a high (above 1.5 mol%) concentration of La2O3 in LaSZ. Besides, numerous micro-cracks between the neighboring grains are observed in Figure 6, whose emergence is due to the volume change of ZrO2 transforming from the tetragonal to the monoclinic phase during cooling.

Simulation and Discussion
As described in the previous study [7], the sintering process can be controlled by the mass transfer. The concentrations of rare earth elements used in sintering were small; therefore, the ion diffusion coefficient and the grain boundary mobility could be assumed to be the same for different rare earth element-doped ZrO2. Figure 7 depicts the pore shrinking during the high temperature treatment. During high temperature annealing, pores will shrink such that the surface area of pores (red line in Figure 7) generally decreases via mass transport through the grain boundary (black line in Figure 7). As a result, the dihedral angle will change to a certain value at which the surface energy and grain boundary energy are in a balanced state. Obviously, the elimination of two surfaces is at the expense of the creation of a grain boundary. As a result, the grain growth and pore shrinkage were determined by the grain boundary energy and surface energy. Moreover, the grain growth rate and pore stability in the powders can be described using the following equations [7,30]:

Simulation and Discussion
As described in the previous study [7], the sintering process can be controlled by the mass transfer. The concentrations of rare earth elements used in sintering were small; therefore, the ion diffusion coefficient and the grain boundary mobility could be assumed to be the same for different rare earth element-doped ZrO 2 . Figure 7 depicts the pore shrinking during the high temperature treatment. During high temperature annealing, pores will shrink such that the surface area of pores (red line in Figure 7) generally decreases via mass transport through the grain boundary (black line in Figure 7). As a result, the dihedral angle will change to a certain value at which the surface energy and grain boundary energy are in a balanced state. Obviously, the elimination of two surfaces is at the expense of the creation of a grain boundary. As a result, the grain growth and pore shrinkage were determined by the grain boundary energy and surface energy. Moreover, the grain growth rate and pore stability in the powders can be described using the following equations [7,30]: where r 0 is the initial average crystal size before heat treatment; r t is the average grain size after heat treatment at high temperatures for time 't'; M GB , γ GB , and γ surf are the grain boundary mobility, grain boundary energy, and surface energy, respectively; and ϕ is the dihedral angle. where r0 is the initial average crystal size before heat treatment; rt is the average grain size after heat treatment at high temperatures for time 't'; MGB, γGB, and γsurf are the grain boundary mobility, grain boundary energy, and surface energy, respectively; and φ is the dihedral angle. To study the grain boundary segregation of rare earth ions and oxygen vacancies via molecular dynamics simulation, the most stable configuration was selected. The binding energy between Y 3+ and VO (oxygen vacancy) is displayed in Figure 8. The binding energy To study the grain boundary segregation of rare earth ions and oxygen vacancies via molecular dynamics simulation, the most stable configuration was selected. The binding energy between Y 3+ and V O (oxygen vacancy) is displayed in Figure 8. The binding energy increases with the decrease of the Y 3+ -V O distance until the Vo occupied the second nearest neighbor site of the Y 3+ ion. The smallest energy difference was −0.474 eV, which shows a good agreement with the references' reports [31]. The Y 3+ ion with a large ionic size (compared with Zr 4+ ) preferred to have eight neighboring oxygen ions, whereas the Vo tended to occupy the nearest neighbor site of the Y 3+ ion. As a result, the optimal site of V O in YSZ is the second nearest neighbor site to accommodate atomic strain and Coulombic interaction. Moreover, the binding energy of the Y 3+ -Y 3+ defect cluster is −0.23 eV, verifying that the rare earth ions also prefer to form a cluster to decrease the atomic strain induced by the ionic radial mismatch. The most stable configuration of YSZ is the configuration with the Y 3+ -V O -Y 3+ defect cluster. To study the grain boundary segregation of rare earth ions and oxygen vacancies via molecular dynamics simulation, the most stable configuration was selected. The binding energy between Y 3+ and VO (oxygen vacancy) is displayed in Figure 8. The binding energy increases with the decrease of the Y 3+ -VO distance until the Vo occupied the second nearest neighbor site of the Y 3+ ion. The smallest energy difference was −0.474 eV, which shows a good agreement with the references' reports [31]. The Y 3+ ion with a large ionic size (compared with Zr 4+ ) preferred to have eight neighboring oxygen ions, whereas the Vo tended to occupy the nearest neighbor site of the Y 3+ ion. As a result, the optimal site of VO in YSZ is the second nearest neighbor site to accommodate atomic strain and Coulombic interaction. Moreover, the binding energy of the Y 3+ -Y 3+ defect cluster is −0.23 eV, verifying that the rare earth ions also prefer to form a cluster to decrease the atomic strain induced by the ionic radial mismatch. The most stable configuration of YSZ is the configuration with the Y 3+ -VO-Y 3+ defect cluster. It can be deduced from Equation (4) that the suppressed grain growth in LnSZ powders originates from the reduction of the grain boundary energy due to the dopant segregation. This coincides with the experimental results, revealing the grain growth suppression by rare earth dopants with large ionic radii. Figure 9 demonstrates that, among the rare earth element-based dopants, La 2 O 3 is found to be the most effective to reduce the surface energy without considering the precipitation of pyrochlore LZ.
Equation (5) reveals the equilibrium relationship between the grain boundary and the crystal surface. However, the dihedral angle between grains is influenced by the grain boundary energy and the surface energy, and the segregation of rare earth elements can reduce the grain boundary energy and surface energy. As a result, the ultimate stability of the dihedral angle around the pores increases owing to a rapidly reducing surface energy, which further increases the pore stability after heat treatment. Figure 9 displays the variances in energies associated with grain boundaries and surfaces undergoing the segregation of rare earth dopants. The segregation of elements with large ionic radii at the grain boundary or at the surface was found to effectively reduce their energies, and the larger the ionic radius, the lower the grain boundary (surface) energy. Furthermore, Figure 9 illustrates that the rare earth dopants with larger ionic radii prefer to segregate on the surface rather than at the grain boundary, because, the larger the rare earth element, the stronger the difference between the grain boundary energy and the surface energy. As a result, the dihedral angle around the pores increases with the ionic radius of the rare earth element, resulting in the higher porosity in Gd 2 O 3 -doped ZrO 2 compared with YSZ. Figure 8. The binding energy between Y 3+ and VO with respect to the ionic distance.
It can be deduced from Equation (4) that the suppressed grain growth in LnSZ powders originates from the reduction of the grain boundary energy due to the dopant segregation. This coincides with the experimental results, revealing the grain growth suppression by rare earth dopants with large ionic radii. Figure 9 demonstrates that, among the rare earth element-based dopants, La2O3 is found to be the most effective to reduce the surface energy without considering the precipitation of pyrochlore LZ. Figure 9. The differences between the grain boundary energy and surface energy in rare earthdoped ZrO2. Equation (5) reveals the equilibrium relationship between the grain boundary and the crystal surface. However, the dihedral angle between grains is influenced by the grain boundary energy and the surface energy, and the segregation of rare earth elements can reduce the grain boundary energy and surface energy. As a result, the ultimate stability of the dihedral angle around the pores increases owing to a rapidly reducing surface energy, which further increases the pore stability after heat treatment. Figure 9 displays the variances in energies associated with grain boundaries and surfaces undergoing the segregation of rare earth dopants. The segregation of elements with large ionic radii at the grain boundary or at the surface was found to effectively reduce their energies, and the larger the ionic radius, the lower the grain boundary (surface) energy. Furthermore, Figure 9 illustrates that the rare earth dopants with larger ionic radii prefer to segregate on the surface rather than at the grain boundary, because, the larger the rare earth element, the stronger the difference between the grain boundary energy and the surface energy. As a result, the dihedral angle around the pores increases with the ionic radius of the rare earth element, resulting in the higher porosity in Gd2O3-doped ZrO2 compared with YSZ.
The aforementioned discussion indicates that the most stable configuration of YSZ is the configuration with the Y 3+ -VO-Y 3+ defect cluster, and the decomposition of rare earthdoped ZrO2 could be ascribed to the formation of RE 3+ -Vo-RE 3+ clusters. From this perspective, the stability of rare earth-doped ZrO2 could be evaluated by using the molecular dynamics simulations. The decomposition energy is defined in terms of the following relationship: The aforementioned discussion indicates that the most stable configuration of YSZ is the configuration with the Y 3+ -V O -Y 3+ defect cluster, and the decomposition of rare earth-doped ZrO 2 could be ascribed to the formation of RE 3+ -Vo-RE 3+ clusters. From this perspective, the stability of rare earth-doped ZrO 2 could be evaluated by using the molecular dynamics simulations. The decomposition energy is defined in terms of the following relationship: where E decom is the decomposition energy and E Cluster and E RSZ are the energy of the decomposed RESZ and the initial configuration of RESZ, respectively. Figure 10 exhibits that the decomposition energy of RESZ shows a nearly linear relationship with the atomic size of RE 3+ ions; that is, it increases with the increase of the atomic size of rare earth elements. In other words, the phase stability of RESZ decreases with the increase of dopant size, which agrees well with the above-mentioned experimental results.
where Edecom is the decomposition energy and ECluster and ERSZ are the energy of the decomposed RESZ and the initial configuration of RESZ, respectively. Figure 10 exhibits that the decomposition energy of RESZ shows a nearly linear relationship with the atomic size of RE 3+ ions; that is, it increases with the increase of the atomic size of rare earth elements. In other words, the phase stability of RESZ decreases with the increase of dopant size, which agrees well with the above-mentioned experimental results. Figure 10. The decomposition energy of rare earth-stabilized ZrO2 calculated using lattice dynamics simulation.

Conclusions
The sintering behavior and phase stability of tetragonal ZrO2 doped with rare earth elements with large ionic radii were thoroughly investigated, and the key results of this study are as follows: Figure 10. The decomposition energy of rare earth-stabilized ZrO 2 calculated using lattice dynamics simulation.