Calculation of Thermal Expansion Coefficient of Rare Earth Zirconate System at High Temperature by First Principles

Compounds of rare earth zirconates with pyrochlore structure are candidates for the application of thermal barrier coatings of next generation. In order to modify the mechanic properties and maintain the low thermal conductivity, other trivalent rare-earth element substitution is commonly used. Presently, investigation on the evaluation of the property of thermal expansion is attracting more attention. In this paper, a feature parameter of thermal expansion coefficient at high temperature (α∞) was proposed by combining Grüneisen’s equation and the Debye heat capacity model. Using α∞ model, the thermal expansion property of different compounds can be easily figured out by first principles. Firstly, α∞ of ZrO2, HfO2, were calculated, and results are in good agreement with the experimental data from the literature. Moreover, α∞ of La2Zr2O7, Pr2Zr2O7, Gd2Zr2O7, and Dy2Zr2O7 were calculated, and results demonstrated that the model of α∞ is a useful tool to predict the thermal expansion coefficient at high temperature. Finally, Gd2Zr2O7 with 4 different Yb dopant concentrations (Gd1-xYbx)2Zr2O7 (x = 0, 0.125, 0.3125, 0.5) were calculated. Comparing with the experimental data from the literature, the calculation results showed the same tendency with the increasing of Yb concentration.


Introduction
It is well known that yttria-stabilized zirconia with a mass fraction of 8%(8YSZ) is widely used as top coating of thermal barrier coatings (TBCs) for aero engines [1,2]. However, with the continuously increasing demand of the thrust-to-weight ratio, 8YSZ is not workable because of its phase change and sintering since the temperature of the turbine front inlet is much greater than 1200 • C [3,4]. Gd 2 Zr 2 O 7 , as the representative of rare earth zirconates, is evidently one of the most promising candidates for the application of next generation TBCs due to its lower thermal conductivity and higher phase stability [5,6]. However, it still suffers from the problem that mechanical properties are not high enough and thermal cycling performance is poor [7]. The thermal expansion properties play the key role. Comparatively, the thermal expansion coefficient of rare earth zirconates is about 9-10 × 10 −6 K −1 (1073 K) [8], which is much lower than that of 8YSZ, about 11 × 10 −6 K −1 (1073 K) [3]. The thermal expansion coefficient of NiCoCrAlY bonding layer, is about 17.5 × 10 −6 K −1 (1273 K) [5]. The mismatch of the thermal expansion between the ceramic top layer and the bonding layer causes thermal stresses during thermal cycling, which can lead to cracks and failure of the TBCs system [9,10].
In order to improve the thermal expansion property of rare earth zirconates, doping with other trivalent rare earth elements is typically used. Excellent thermophysical properties such as high thermal stability, lower thermal conductivity, and high thermal expansion have been demonstrated in doped rare earth zirconates, such as Sm 2 (Zr 1-x Ce x ) 2 O 7 , (Nd 1-x Gd x ) 2 Zr 2 O 7 , Gd 2 (Zr 1-x Ti x ) 2 O 7 and (Sm 1-x Gd x )Zr 2 O 7 [6,[11][12][13], etc. Nevertheless, the research described above was carried out mainly by experiments, which are time costly and not fully adequate. For instance, heat conduction and thermal expansion of materials are closely related to lattice vibrations. However, to investigate the lattice vibrations of materials experimentally, neutron scattering or Raman spectroscopy is required [14]. Correspondingly, computational material simulations are more efficiently implemented to accelerate the design of materials. As we know, the first principles exhibit powerful capabilities in material design. Zhao [15] investigated the structure, mechanical properties, minimum thermal conductivity, and electronic properties of a series of Gdsite and Zr-site substituted Gd 2 Zr 2 O 7 pyrochlores by first principles. Atsushi Togo [16] calculated the thermal expansion properties of Ti 3 SiC 2 , Ti 3 AlC 2 , and Ti 3 GeC 2 by the first principles combing with quasi-harmonic approximation (QHA). Feng [17] investigated the thermal expansion properties of rare earth zirconates (Ln 2 Zr 2 O 7 , Ln = La, Nd, Sm and Gd) of pyrochlore structures by using first principles. Meanwhile, the calculation of properties of doped compounds by first principles is very difficult especially for the calculation of the thermal expansion property. For instance, the conventional cell of Gd 2 Zr 2 O 7 contains 88 atoms, and cell expansion must be performed to obtain solid solution structures with different doping concentrations, which makes the calculation of phonons extremely difficult.
Classically, Grüneisen proposed a thermal expansion equation [18], in which the thermal expansion coefficient(α) is determined by elastic properties and heat capacity (C V ). Further, C V is the function of Debye temperature and temperature. When the temperature is much greater than Debye temperature, C V can be considered as a constant. In this case, the calculation of α is simplified to the calculation of elastic properties, which makes the calculation much easier and faster. Coincidentally, the TBCs for aero engines operate at a high temperature, which is much higher than Debye temperature. Moreover, investigations showed that the coefficient of thermal expansion gradually increases with temperature increasing at high temperature, but the increasing rate gradually decreases. Reasonably, the thermal expansion coefficient at super high temperature(α ∞ ) can be used as a comparison factor to characterize the thermal expansion property of different dopant/concentration compound materials. Therefore, in this paper, based on the Grüneisen's equation, a computational model of α ∞ was developed, by which the calculation of a series of rare earth zirconates were implemented by the first principles.

Methodology
The supercell containing 88 atoms of Gd 2 Zr 2 O 7 was used for calculation. The structures of Yb doped Gd 2 Zr 2 O 7 were formed by replacing Gd atoms with different amounts of Yb atoms. The structural models for (Gd 1-x Yb x ) 2 Zr 2 O 7 were built using the cluster expansion approach by calculating the lowest forming energy [19][20][21]. Further, the structures were optimized by the Birch-Murnaghan equation of state [22]. The elastic constants of the material were calculated by the stress-strain method [23].
The first principles calculations were based on density functional theory using the Vienna Ab initio Simulation Package (VASP) [24] with the generalized gradient approximation (GGA) for exchange-correlation energy, in the form of Perdew-Burke-Ernzerhof (PBE) [25]. The kinetic cut-off energy for the plane wave expansion was taken to be 600 eV in the Brillouin zone integrations using 2 × 2 × 2 k points. The average force acting on ions was reduced to 0.05 eV/Ang. Valence electrons included for distinct atoms were O 2s 2 2p 4 , Zr 5s 1 4d 3 , Hf 5d 3 6s 1 , Gd 6s 2 5p 6 5d 1 , Yb 6s 2 5p 6 .

Yb Doped Gd 2 Zr 2 O 7 Structure
The supercell of Gd 2 Zr 2 O 7 contains 88 atoms, including 16 Gd atoms, 16 Zr atoms and 56 oxygen atoms. A total of 2, 5 and 8 Yb atoms were used to replace the Gd atoms to obtain three different concentrations: (Gd 0.875 Yb 0.125 ) 2 Zr 2 O 7 , (Gd 0.875 Yb 0.3125 ) 2 Zr 2 O 7 , and (Gd 0.5 Yb 0.5 ) 2 Zr 2 O 7 , respectively. Correspondingly, the possible numbers of (Gd 1-x Yb x ) 2 Zr 2 O 7 doped structures were C 2 16 , C 5 16 , and C 8 16 . Excluding the equivalent structures, the number of unequal possible doped structures are 3, 35, and 97, separately. The forming energy E of each structure was calculated using the cluster expansion approach, according to Equation (1) [21], and the final doped structure was identified by the structure with the lowest forming energy.
wherein E 0 is the energy of the doped structure; E 1 and E 2 are the energy of the single cell of Gd 2 Zr 2 O 7 and Yb 2 Zr 2 O 7 , and n 1 and n 2 are the number of dopant atoms, respectively. The calculation results were shown in Figure 1. According to the lowest forming energy, the geometrical configurations of three different doping structures were elaborated in Figure 2. ions was reduced to 0.05 eV/Ang. Valence electrons included for distinct atoms were O 2s 2 2p 4 , Zr 5s 1 4d 3 , Hf 5d 3 6s 1 , Gd 6s 2 5p 6 5d 1 , Yb 6s 2 5p 6 .

Yb Doped Gd2Zr2O7 Structure
The supercell of Gd2Zr2O7 contains 88 atoms, including 16 Gd atoms, 16 Zr atoms and 56 oxygen atoms. A total of 2, 5 and 8 Yb atoms were used to replace the Gd atoms to obtain three different concentrations: (Gd0.875Yb0.125)2Zr2O7, (Gd0.875Yb0.3125)2Zr2O7, and (Gd0.5Yb0.5)2Zr2O7, respectively. Correspondingly, the possible numbers of (Gd1-xYbx)2Zr2O7 doped structures were 2 16 C , 5 16 C , and 8 16 C . Excluding the equivalent structures, the number of unequal possible doped structures are 3, 35, and 97, separately. The forming energy E of each structure was calculated using the cluster expansion approach, according to Equation (1) [21], and the final doped structure was identified by the structure with the lowest forming energy.
wherein E0 is the energy of the doped structure; E1 and E2 are the energy of the single cell of Gd2Zr2O7 and Yb2Zr2O7, and n1 and n2 are the number of dopant atoms, respectively. The calculation results were shown in Figure

Lattice Constant and Elastic Modulus
Based on the geometrical configurations shown in Figure 2, the lattice constants and elastic modulus were calculated and listed in Table 1. Table 1. Lattice constant, the elastic constants (C 11 , C 12 , and C 44 ), bulk modulus (B), shear modulus (G), and Poisson's ratio (µ) of rare earth zirconates.

Lattice Constant and Elastic Modulus
Based on the geometrical configurations shown in Figure 2, the lattice constants and elastic modulus were calculated and listed in Table 1. It can be seen that the lattice constants decrease with the increasing of Yb dopant. The calculated bulk and shear modulus of Gd2Zr2O7 are 176.6 Gpa and 91.9 GPa, respectively, which meets agreement with the data measured by experiments [26]. Generally, the bulk modulus and shear modulus of (Gd1-xYbx)2Zr2O7 decrease with the increasing of Yb content, which are possibly caused by structure change. Subramanian M A [27] pointed out that doping of Yb atoms reduce the average cation radius ratio r(A 3+ )/r(B 4+ ) and change the crystal structure of Gd2Zr2O7 from pyrochlore to disorders in the structures. For the cubic phase, there are three independent elastic constants, C11, C12, and C44 [28]. The calculated elastic constants are all positive, satisfying the generalized elastic stability criterion, namely, C11 + 2C12 > 0; C44 > 0; C11 − C12 > 0, indicating that all studied structures are mechanically stable [29]. According to Pugh's theory [30], when G/B < 0.5, the material is ductile; otherwise, the material is brittle. The G/B value of all materials calculated is greater than 0.5, indicating that they are brittle materials. It can be seen that the lattice constants decrease with the increasing of Yb dopant. The calculated bulk and shear modulus of Gd 2 Zr 2 O 7 are 176.6 Gpa and 91.9 GPa, respectively, which meets agreement with the data measured by experiments [26]. Generally, the bulk modulus and shear modulus of (Gd 1-x Yb x ) 2 Zr 2 O 7 decrease with the increasing of Yb content, which are possibly caused by structure change. Subramanian M A [27] pointed out that doping of Yb atoms reduce the average cation radius ratio r(A 3+ )/r(B 4+ ) and change the crystal structure of Gd 2 Zr 2 O 7 from pyrochlore to disorders in the structures. For the cubic phase, there are three independent elastic constants, C 11 , C 12 , and C 44 [28]. The calculated elastic constants are all positive, satisfying the generalized elastic stability criterion, namely, C 11 + 2C 12 > 0; C 44 > 0; C 11 − C 12 > 0, indicating that all studied structures are mechanically stable [29]. According to Pugh's theory [30], when G/B < 0.5, the material is ductile; otherwise, the material is brittle. The G/B value of all materials calculated is greater than 0.5, indicating that they are brittle materials.

Thermal Expansion of Rare Earth Zirconates System
According to Grüneisen's equation [18], the volumetric coefficient of thermal expansion (β) can be expressed as Equation (2).
Here, γ is Grüneisen's constant, C V is the heat capacity, B is the bulk modulus, and V is the molar volume.
As we know, the rare earth zirconate is in cubic phase. For cubic phase, γ is defined as function of the Poisson's ratio (µ) as the following Equation (3) [31].
wherein N = nN A , N A is Avogadro's constant, and n is the number of atoms in the molecular formula; k B is the Boltzmann constant; T D is the Debye temperature. Meanwhile, T D can be calculated by the following Equation (5) [33].
υ L , υ S is the longitudinal and transverse sound velocity and can be expressed as Equations (7) and (8), separately.
According to the calculation result shown in Table 1, B, µ, and V of different (Gd 1-x Yb x ) 2 Zr 2 O 7 vary within 4% difference, which can be considered as constant. Thus, it is indicated from Equation (10) that α is proportional to T D −3 at the same temperature, which is compliant with Ruffa's equation [34].
Furthermore, supposing that the temperature is much greater than the Debye temperature, the integral term in Equation (9) can be mathematically simplified and finally Equation (11) can be obtained.
Here, α ∞ representatives the linear thermal expansion coefficient at super high temperature. Actually, the TBCs are working under the temperature (e.g., the temperature of combustion chamber in F135 turbine engine can be up to 2253 K [35]) much higher than T D of rare earth zirconate (about 500 K) [36]. Another one, α of rare earth zirconate, increases with the increase in temperature, meanwhile the increasing rate slows down more and more [37][38][39]. Therefore, α ∞ can be likely used to compare the difference of thermal expansion property of different dopant/concentration for the same compound.

The Validity of α ∞ Model
Cubic ZrO 2 and HfO 2 are of typical fluorite structure which is the same as that of rare earth zirconates [40]. Firstly, the lattice parameter and elastic properties were calculated by first principles. Secondly utilizing the α ∞ model, the thermal expansion properties of cubic ZrO 2 and HfO 2 were calculated, and both results were listed in Table 2. Comparing to the data from the material project database, calculation results of lattice parameter and elastic properties are very close to the same level. Hong [41] and Irshad, K.A. [42] measured the linear thermal expansion coefficient (α) of cubic ZrO 2 and HfO 2 by in situ high-temperature X-ray diffraction, and it is (12 ± 3) × 10 −6 K −1 and 8.80 × 10 −6 K −1 , respectively. Comparably, the calculated results of α ∞ are 9.72 × 10 −6 K −1 and 9.05 × 10 −6 K −1 . It is revealed that both are a good match. In order to further verify the validity of the α ∞ model, α ∞ of serial rare earth zirconates, including La 2 Zr 2 O 7 , Pr 2 Zr 2 O 7 , Gd 2 Zr 2 O 7 , and Dy 2 Zr 2 O 7 were calculated. The data of elastic property was cited from the literature [36], and both were shown in Table 3. α ∞ and α were plotted in Figure 3. Compared to α measured by experiment [8], α ∞ are clearly higher because α were measured at 800 • C which is not too much higher than Debye temperature. Henry Lehmann [37] measured the thermal expansion coefficients of Gd 2 Zr 2 O 7 and La 2 Zr 2 O 7 , which were 10.652 × 10 −6 K −1 (1473 K) and 9.09 × 10 −6 K −1 (1373 K), respectively. The results are very close to α ∞ calculated. It is further demonstrated that α ∞ can be a useful tool to predict the thermal expansion coefficient at high temperature. Table 3. Lattice constant a 0 , bulk modulus (B), shear modulus (G), Poisson's ratio (µ) of rare earth zirconates system [36].

The Effect of Yb Doping of Gd2Zr2O7 on α∞
α∞ of Gd2Zr2O7 with 4 different Yb doping contents Gd2Zr2O7, (Gd0.875Yb0.125)2Zr2O7, (Gd0.875Yb0.3125)2Zr2O7, and (Gd0.5Yb0.5)2Zr2O7 were calculated. Table 4 shows the calculation result of theoretical density, sound velocity, and Debye temperature. Table 4. Density (ρ), longitudinal wave velocity (υL), shear wave velocity of sound (υS), average velocity of sound (υm) and Debye temperature (TD) of rare earth zirconates. It can be seen that TD of Gd2Zr2O7 was calculated to be 511 K, which is in good agreement with that measured by Toshiaki Kawano [8]. TD is decreases with the increase of Yb dopant, which is dominated by the decrease of average velocity of sound (υm). Figure 4 shows the difference between calculated α∞ and α measured by experiment at 1073 K [43].   It can be seen that T D of Gd 2 Zr 2 O 7 was calculated to be 511 K, which is in good agreement with that measured by Toshiaki Kawano [8]. T D is decreases with the increase of Yb dopant, which is dominated by the decrease of average velocity of sound (υ m ). Figure 4 shows the difference between calculated α ∞ and α measured by experiment at 1073 K [43]. The content of Yb exp. [43] cal. It is clear that Yb doping can increase the thermal coefficient greatly because theoretical density increases meanwhile TD decreases with Yb doping, as shown in Table 4. However, α∞ is a bit lower than α. Feng [17] also discussed and explained the problem. Actually, the first principles were developed based on the material being from an ideally perfect crystal. However, for real bulk materials, various defects (e.g., vacancies and dislocations) and pores ineluctably existed. In general, the total energy of a crystal with defects is higher than that of an ideally perfect crystal, and the anharmonic effect may be affected It is clear that Yb doping can increase the thermal coefficient greatly because theoretical density increases meanwhile T D decreases with Yb doping, as shown in Table 4. However, α ∞ is a bit lower than α. Feng [17] also discussed and explained the problem. Actually, the first principles were developed based on the material being from an ideally perfect crystal. However, for real bulk materials, various defects (e.g., vacancies and dislocations) and pores ineluctably existed. In general, the total energy of a crystal with defects is higher than that of an ideally perfect crystal, and the anharmonic effect may be affected by various defects and pores in the structure. In addition, the density of the tested ceramic coupon is certainly lower than that of theoretical density. That is the reason why the thermal expansion coefficient measured by experiment is higher than that calculated. Figure 4 also reveals that both α ∞ and α increase with an increase in Yb content. The increase of thermal coefficient is more remarkable with lower Yb doping concentration. With higher Yb doping concentration, the growth rate slows down. Both the measured and calculated results have the same tendency.

ρ/(kg·m −3 ) υL/(m·s −1 ) υS/(m·s −1 ) υm/(m·s −1 ) TD/(K)
Thermal expansion is related to the crystal structure and the electronic structure [44]. The PDOS of (Gd 1-x Yb x ) 2 Zr 2 O 7 crystals are shown in Figure 5. In terms of bonding, Gd/Yb-5d, Zr-4d, and O-2p overlap, which means that the electronic state around the Fermi level is primarily determined by the relatively weak p-d bond between O-2p and Zr-4d (or Gd/Yb-5d). The main change from doping is the relative position changes in Gd/Yb-5d, which may be caused by the difference in the valence electron layers between Gd and Yb. With the increase of Yb content, the total state density curve moves slightly to the lower energy level. Low crystal energy means high coefficient of thermal expansion [45], and the rare earth element Yb influences the O-2p states due to hybridization [44]. The p-d bond strength decreased with the increase in Yb content. According to the value of the PDOS ordinate, the PDOS of Zr-4d is the largest, so the p-d bond strength of Zr-O is higher than Yb-O and Gd-O.

Conclusions
Rare earth zirconates are candidates for next generation TBCs, and it is important to develop a method to characterize the thermal expansion property. Combining the Grüneisen's equation and the Debye heat capacity model, an efficient model of α∞ to charac-

Conclusions
Rare earth zirconates are candidates for next generation TBCs, and it is important to develop a method to characterize the thermal expansion property. Combining the Grüneisen's equation and the Debye heat capacity model, an efficient model of α ∞ to characterize the thermal expansion coefficient at super high temperature was established. Firstly, using the α ∞ model, the high temperature thermal expansion coefficients of cubic ZrO 2 and cubic HfO 2 were calculated to be 9.72 × 10 −6 K −1 and 9.05 × 10 −6 K −1 , respectively, which are in agreement with those shown in the literature. Secondly, α ∞ of serial rare earth zirconates, including La 2 Zr 2 O 7 , Pr 2 Zr 2 O 7 , Gd 2 Zr 2 O 7 , and Dy 2 Zr 2 O 7 were calculated, and results demonstrated that α ∞ can be a useful tool to predict the thermal expansion coefficient at high temperature. Lastly, α ∞ of (Gd 1-x Yb x ) 2 Zr 2 O 7 with four different doping contents were calculated, and results showed the same tendency as that measured by experiments. Generally, by characterizing the thermal expansion coefficient at high temperature through the elastic properties and Debye temperature of the material, the complicated calculation of phonon spectrum can be avoided. Thus, the model of α ∞ has the broad application prospect to predict the thermal expansion property at high temperature for other rare earth zirconates.

Data Availability Statement:
The data presented in this study are available upon request from the corresponding author.