I. INTRODUCTION
Yttrium disilicate (Y2Si2O7) is well known for its polymorphism and exhibits up to six structure types, known as у, α, β, γ, δ, and probably z.Reference Felsche1 Experimental trials and theoretical calculations have found some interesting properties of Y2Si2O7 such as high melting point (1775 °C), good machinability, stability in oxidizing and hot corrosive environments, and low thermal conductivity, which guarantee its potential application as a candidate for high-temperature structural materials or environmental barrier coatings (EBCs).Reference Sun, Zhou, Wang and Li2–Reference Sun, Li and Zhou5 For example, it was reported that Y2Si2O7 is an important grain boundary phase in Si3N4 sintered with Y2O3 additive and endows good high-temperature properties.Reference Hong, Yoshida, Ikuhara, Sakuma, Nishimura and Mitomo6 Sun's research revealed that Y2Si2O7/ZrO2 composite with improved mechanical properties and wear resistance is a candidate for advanced EBCs against thermal load, chemical attack, and humid hot-gas corrosion.Reference Sun, Wu, Li and Zhou7 However, a crucial challenge rising from the polymorphism of Y2Si2O7 is that phase transitions with temperature variations might result in undesirable changes in crystal structure and properties, which will further lead to catastrophic degradation of materials in real application environment.Reference Fernandez-Carrion, Allix, Florian, Suchomel and Becerro8 Therefore, it is necessary to comprehensively investigate and understand the relationship between crystal structure and elastic/thermal properties of different Y2Si2O7 phases.
Polymorphism is a common phenomenon among rare-earth disilicates.Reference Fleet and Liu9,Reference Maier, Rixecker and Nickel10 In the present work, we choose Y2Si2O7 as a model system to explore how elastic stiffness and thermal conductivity react to the changes of chemical bonding and crystal structure in different polymorphs. Among all the six polytypes of Y2Si2O7, four high temperature phases, α-, β-, γ-, and δ-Y2Si2O7, attract the most interests. The phase transition path is described as: 
$\alpha \buildrel {1250 \pm 10{}^{\rm{o}}{\rm{C}}} \over \longrightarrow \beta \buildrel {1445 \pm 10{}^{\rm{o}}{\rm{C}}} \over \longrightarrow \gamma \buildrel {1535 \pm 10{}^{\rm{o}}{\rm{C}}} \over \longrightarrow \delta$.Reference Ito and Johnson11 The β, γ, and δ phases contain similar Si2O7 pyrosilicate units; therefore, they constitute a suitable polymorphic family for the study of relationship between crystal structure and elastic/thermal properties of Y2Si2O7 disilicates. Whilst, α-Y2Si2O7 has isolated SiO4 tetrahedra and Si3O10 units,Reference Parmentier, Bodart, Audoin, Massouras, Thompson, Harris, Goursat and Besson12 which are obviously different with the crystal structures of other polymorphic members. Therefore, α-Y2Si2O7 is not included in the present study.
Our previous investigation has revealed the mechanisms dominating the low thermal conductivity of γ-Y2Si2O7.Reference Luo, Wang, Wang, Li and Hu13 The key factors include significant bonding heterogeneity, high level of anharmonic phonon scattering, and low group velocity of acoustic phonon transport. In this study, we take a step further and focus on comparing the intrinsic elastic/thermal properties among β-, γ-, and δ-Y2Si2O7 from the bases of the complexity of crystal structure, phonon group velocity, and anharmonicity of phonon scattering. We first calculate the equilibrium crystal structure, second-order elastic coefficients, bulk modulus, shear modulus, Young's modulus, and Poisson's ratio of β-, γ-, and δ-Y2Si2O7. Thereafter, mean sound velocity, Debye temperature, and theoretical minimum thermal conductivity are obtained from the calculated structural and mechanical parameters using the Debye–Clarke model. The temperature dependences of the lattice thermal conductivities are also predicted from the Debye–Slack model. The results demonstrate how the structural and elastic/thermal parameters significantly affect the thermal conductivities of polymorphic Y2Si2O7.
II. CALCULATION METHOD
A. Structural properties
First-principles calculations are performed based on the density functional theory using the CASTEP code, wherein the Vanderbilt-type ultrasoft pseudopotential and the localized density approximation (LDA) are incorporated. After convergence test, we fix the plane-wave basis set cutoff energy at 450 eV in all calculations. The special points sampling integration over Brillouin zone are realized by using the Monkhorst–Pack method with special k-points meshes of 4 × 3 × 5, 5 × 2 × 5, and 2 × 5 × 3 for β-, γ-, and δ-Y2Si2O7, respectively. Equilibrium crystal structures of studied materials are obtained from adequate geometry optimization. During this process, lattice parameters and internal atomic coordinates are independently modified to minimize the total energy, interatomic forces, and unit-cell stresses. In this work, the Broyden–Fletcher–Goldfarb–Shanno (BFGS) minimization scheme is used and the convergence tolerances are set as follows: The difference in total energy within 5.0 × 10−6 eV per atom, the maximum ionic Hellmann–Feynman force within 0.01 eV/Å, the maximum stress within 0.02 GPa, and the maximum ionic displacement within 5.0 × 10−4 Å. The reliability of present calculations had been verified by previous works.Reference Wang, Zhou and Lin4,Reference Luo, Wang, Wang, Li and Hu13
Degree of distortion (Δd) of a structural unit could be estimated by the deviations of individual bond lengths from their mean value, expressed in the form:Reference Cong, Zhang, Wang, Yu, Fan, Cheng, Sun, Zhang, Lu, Jiang and Boughton14
$${\Delta _{\rm{d}}} = {1 \over t}{\sum\limits_i {\left[ {{{d{{\left( {{\rm{M}} - {\rm{O}}} \right)}_i} - \overline {d\left( {{\rm{M}} - {\rm{O}}} \right)} } \over {\overline {d\left( {{\rm{M}} - {\rm{O}}} \right)} }}} \right]} ^2}\quad {\rm{,}}$$where t is the total number of Si–O bonds in a SiO4 tetrahedron or Y–O bonds in a YO6 or YO7 polyhedron; d(M–O)i and 
$\overline {d\left( {M - O} \right)}$ are the i th and average length of Si–O or Y–O bonds, respectively. A higher value of ∆d represents a more distorted structural unit.
B. Elastic properties
As the first step to investigate the elastic stiffness of β-, γ-, and δ-Y2Si2O7, second-order elastic coefficients are calculated from the stress–strain relationship. We apply four homogeneous strain patterns and introduce six amplitudes (three positive and three negative) for each pattern with the maximum strain value |ε| of 0.5%. The corresponding stresses are calculated with respect to optimizing the internal atomic degrees of freedom, as is proposed by Milman.Reference Milman and Warren15 The elastic constants are determined from linear fitting of the stress–strain dependences. Criteria of convergence are set as follows: The difference in total energy within 1.0 × 10−6 eV per atom, the maximum ionic Hellmann–Feynman force within 0.002 eV/Å, and the maximum ionic displacement within 1.0 × 10−4 Å.
The compliance tensor matrix [S], inverse of stiffness tensor matrix [C], is calculated from the second order elastic coefficients. Thereafter, the bulk modulus B and shear modulus G could be derived from Voigt–Reuss–Hill (VRH) approximations.Reference Sun, Liu, Wang, Wang, Zhou and Hu16,Reference Hill17
According to Voigt's approximation, the bulk modulus B V and shear modulus G V are expressed as:
$$\eqalign{ & {B_{\rm{V}}} = {1 \over 9}\left( {{c_{11}} + {c_{22}} + {c_{33}}} \right) + {2 \over 9}\left( {{c_{12}} + {c_{13}} + {c_{23}}} \right) \cr & {G_{\rm{V}}} = {1 \over {15}}\left( {{c_{11}} + {c_{22}} + {c_{33}} - {c_{12}} - {c_{13}} - {c_{23}}} \right) <$$><$$>+ {1 \over 5}\left( {{c_{44}} + {c_{55}} + {c_{66}}} \right)\quad {\rm{.}} \cr}$$According to Reuss's approximation, B R and G R are expressed as:
$$\eqalign{ & {B_{\rm{R}}} = {1 \over {\left( {{s_{11}} + {s_{22}} + {s_{33}}} \right) + 2\left( {{s_{12}} + {s_{13}} + {s_{23}}} \right)}} \cr & {G_{\rm{R}}} = {{15} \over {4\left( {{s_{11}} + {s_{22}} + {s_{33}}} \right) - 4\left( {{s_{12}} + {s_{13}} + {s_{23}}} \right) + 3\left( {{s_{44}} + {s_{55}} + {s_{66}}} \right)}}\quad {\rm{.}} \cr}$$Voigt's and Reuss's schemes set the upper and lower limits of the elastic modulus for polycrystalline, respectively. Hill proposed to use the mean arithmetic values B H and G H to present the elastic parameters of polycrystalline:Reference Hill17
$$B_{\rm{H}} = {1 \over 2}\left( {B_{\rm{V}} + B_{\rm{R}} } \right)\quad G_{\rm{H}} = {1 \over 2}\left( {G_{\rm{V}} + G_{\rm{R}} } \right)\quad {\rm{.}}$$$$The Young's modulus E and Poisson's ratio ν are determined by:Reference Hill17
$$E = {{9{B_{\rm{H}}}{G_{\rm{H}}}} \over {3{B_{\rm{H}}} + {G_{\rm{H}}}}}\quad {\rm{,}}$$
$$\nu = {{3{B_{\rm{H}}} - 2{G_{\rm{H}}}} \over {2\left( {3{B_{\rm{H}}} + {G_{\rm{H}}}} \right)}}\quad {\rm{.}}$$C. Thermal properties
Lattice thermal conductivity determines the heat transport in a material with temperature gradient.Reference Clarke18 In the successful model of describing lattice thermal conductivity, Debye introduced a kinetic theory, by which a collection of phonons with volumetric heat capacity C v and group velocity νm undergo random walks with mean free paths Λ in a temperature field,Reference Kittel19
$$\kappa = {1 \over 3}{C_{\rm{v}}}{\nu _{\rm{m}}}\Lambda \quad {\rm{.}}$$At high temperatures, C v approaches to the constant value of 3k B per atom.
As is suggested by Clarke,Reference Clarke18 Λ approaches to the average interatomic distance at the high temperature limit when phonons are completely uncoupled, and hence the thermal conductivity reaches to the minimum value. The minimum thermal conductivity (κmin) is a qualitative criterion for identifying novel candidate with low thermal conductivity. Clarke introduced mass approximation by replacing different atoms in a molecule with an equivalent atom having a mean atomic mass of M/m and expresses the theoretical minimum thermal conductivity as:Reference Clarke18
$$\kappa \to {k_{\rm{B}}}{\nu _{\rm{m}}}\Lambda \to {k_{\rm{B}}}{\nu _{\rm{m}}}{\Lambda _{\min }} \to {k_{\rm{B}}}{\nu _{\rm{m}}}{\left( {{M \over {m\rho {N_{\rm{A}}}}}} \right)^{ - {2 \mathord{\left/ {\vphantom {2 3}} \right. \kern-\nulldelimiterspace} 3}}}\quad {\rm{,}}$$where k B is the Boltzmann constant, N A is the Avogadro number, ρ is the density, M is the molecular mass, and m is the number of atoms per molecule.
The average sound velocity v m, almost independent of temperature change, is expressed by:Reference Kittel19
$${\nu _{\rm{m}}} = {\left[ {{1 \over 3}\left( {{2 \over {\mathop {{\nu _{\rm{s}}}}\nolimits^3 }} + {1 \over {\mathop {{\nu _1}}\nolimits^3 }}} \right)} \right]^{ - {1 \mathord{\left/ {\vphantom {1 3}} \right. \kern-\nulldelimiterspace} 3}}}\quad {\rm{,}}$$where v s and v l are the velocities of transverse and longitudinal acoustic waves, respectively, which can be derived from the bulk modulus, shear modulus, and mass density,
$${\nu _{\rm{s}}} = \sqrt {\left( {{{{G_{\rm{H}}}} \over \rho }} \right)} \quad {\rm{and}}\quad {\nu _{\rm{l}}} = \sqrt {{{\left( {{B_{\rm{H}}} + {4 \over 3}{G_{\rm{H}}}} \right)} \mathord{\left/ {\vphantom {{\left( {{B_{\rm{H}}} + {4 \over 3}{G_{\rm{H}}}} \right)} \rho }} \right. \kern-\nulldelimiterspace} \rho }} \quad {\rm{.}}$$Debye temperature is also an important parameter to evaluate thermal properties of a material and could be expressed as:
$${\theta _{\rm{D}}} = {h \over {{k_{\rm{B}}}}}{\left[ {{{3m} \over {4\pi }}\left( {{{\rho {N_{\rm{A}}}} \over M}} \right)} \right]^{{1 \mathord{\left/ {\vphantom {1 3}} \right. \kern-\nulldelimiterspace} 3}}}{\nu _{\rm{m}}}\quad {\rm{,}}$$where h is the Plank constant.
The minimum thermal conductivity serves as a qualitative guideline for the selection of potential thermal insulator at high temperature, but it alone tells too little information about heat transport over a wide temperature range. In the present calculation, the temperature dependence of lattice thermal conductivity is calculated from the well established Slack's equation:Reference Morelli, Slack, Shindé and Goela20
$$\kappa = \matrix{ {A{{\bar M\mathop {\theta _{\rm{D}} }\nolimits^3 \delta } \over {\gamma ^2 n^{{2 \mathord{\left/ {\vphantom {2 3}} \right. \kern-\nulldelimiterspace} 3}} T}}} & {\left( {T \,\ge\, \theta _{\rm{D}} } \right)} \cr } \quad {\rm{,}}$$where 
$\bar M$ [kg/mol] is the mean atomic mass, δ3 [m3] is the average volume of one atom in the primitive unit cell, θD [K] is the Debye temperature, T [K] is the absolute temperature, n is the number of atoms per primitive unit cell, γ is the acoustic Grüneisen constant derived from Poisson's ratio ν,Reference Sanditov, Tsydypov and Sanditov21,Reference Badehian, Salehi and Ghoohestani22 and A [W mol kg−1 m−2 K−3] is a coefficient determined by γ:Reference Morelli, Slack, Shindé and Goela20,Reference Bruls23
$$\gamma = {{9\left( {\mathop {{\nu _1}}\nolimits^2 - {\textstyle{4 \over 3}}\mathop {{\nu _{\rm{s}}}}\nolimits^2 } \right)} \over {2\left( {\mathop {{\nu _1}}\nolimits^2 - 2\mathop {{\nu _{\rm{s}}}}\nolimits^2 } \right)}} = {3 \over 2}\left( {{{1 + \nu } \over {2 - 3\nu }}} \right)\quad {\rm{,}}$$
$$A\left( \gamma \right) = {{5.720 \times {{10}^7} \times 0.849} \over {2 \times \left( {1 - {{0.514} \over \gamma } + {{0.228} \over {{\gamma ^2}}}} \right)}}\quad {\rm{.}}$$There are some approximate formulations for the Grüneisen constant, such as the thermodynamic γ determined from thermal expansion data,Reference Morelli, Slack, Shindé and Goela20 microscopic γ calculated from mode-specific Grüneisen constant,Reference Morelli, Slack, Shindé and Goela20 and acoustic γ obtained from the elastic parameters near equilibrium structure.Reference Sanditov, Tsydypov and Sanditov21 Our previous investigation has shown that acoustic Grüneisen constant deduced from Poisson's ratio yields a reliable prediction of the anharmonicity of yttrium silicates.Reference Luo, Wang, Wang, Li and Hu13 Furthermore, parameters θD, δ, and γ are reasonably approximated as independent of temperature.
III. RESULTS AND DISCUSSION
A. Crystal structure
The calculated lattice parameters of β-, γ-, and δ-Y2Si2O7 are listed in Table I, together with reported experimental data from x-ray diffraction for comparison. Both β- and γ-Y2Si2O7 have the monoclinic lattice but in different space groups, C2/m for β-Y2Si2O7 and P21/C for γ-Y2Si2O7; δ-Y2Si2O7 has the orthorhombic lattice with Pna21 space group. The theoretical lattice parameters of β-Y2Si2O7 deviate from the experimental values by 1.43% for a, 1.33% for b, 1.14% for c, and 0.10% for β; deviations for γ-Y2Si2O7 are 1.30% for a, 1.47% for b, 1.40% for c, and 0.11% for β; and differences for δ-Y2Si2O7 are 1.18% for a, 1.16% for b, and 1.33% for c. Calculated lattice parameters are in good agreement with referenced data.Reference Liddell and Thompson24–Reference Dinger, Rai and Thomas26
Calculated lattice parameters of β-, γ-, and δ-Y2Si2O7, together with experimental data for comparison.
Crystal structures of β-, γ-, and δ-Y2Si2O7 are presented in Fig. 1. The overall characteristics of the three Y2Si2O7 polymorphs could be described as Si2O7 units stacking in the framework of Y–O polyhedra. However, difference arises from detailed investigations. The unit cell of β-Y2Si2O7 contains 22 atoms (two formulas). As is shown in Fig. 1(a), there is one Y site, one Si site, and three O sites. The Y atom bonds to six O atoms and the YO6 octahedron consists of Y–O bonds with lengths in the range of 2.214–2.274 Å. Each YO6 octahedron shares three edges with neighboring octahedra and builds up the distorted honeycomb structure that is parallel to plane ab. Si is bonded with four O atoms via one shorter bond (1.601 Å) and three longer bonds (two of 1.612 Å and one of 1.613 Å). Two adjacent SiO4 tetrahedra share one bridging OB atom and are linked in the Si2O7 pyrosilicate units. The Si–OB–Si bridge is linear (with a bond angle of 180°) and is perpendicular to axis b. The triangular facets of the two SiO4 tetrahedra are oriented in opposite directions in the Si2O7 unit. The unit cell of γ-Y2Si2O7, as is shown in Fig. 1(b), contains 22 atoms (two formulas) with one Si site, one Y site, and four different O sites. The sandwich-like stacking configuration of YO6 octahedra remains, but the honeycomb arrangement breaks. Instead, YO6 polyhedra share edges with each other and form ribbon-like arrangement along plane bc. Apart from these differences, the crystal structure of γ-Y2Si2O7 resembles that of β-Y2Si2O7. Linear Si–OB–Si bridge commonly exists, and the bridging OB atom is located right at the center of symmetry. As is shown in Fig. 1(c), the unit cell of δ-Y2Si2O7 contains 44 atoms (four formulas) with two Y sites, two Si sites, and seven O sites. Y occupies two nonequivalent sites, forming two different types of YO7 polyhedra that alternately stack up in a wave-like arrangement parallel to (010) atomic plane. Bond lengths of Y–O bonds differ slightly in the two sets of YO7 polyhedra (2.193–2.441 Å and 2.194–2.443 Å, respectively). Both Si1 and Si2 atoms are coordinated by four nonequivalent O atoms with Si–O bond lengths ranging from 1.587 to 1.647 Å and 1.595 to 1.655 Å, respectively. The Si2O7 unit consists of two nonequivalent SiO4 tetrahedra, and the Si–OB–Si bond bends with the angle of 158.1°. The bond lengths of two nonequivalent Si–OB bonds are 1.647 and 1.655 Å. What worthy of emphasis is that in β and γ phases, the Si–OB–Si bridge is linear and the bridging OB atom is not bonded with Y atoms; whereas in the δ phase, the Si–OB–Si bridge is nonlinear, and the bridging OB atom is bonded with two neighboring Y atoms.
Crystal structures of (a) β-, (b) γ-, and (c) δ-Y2Si2O7.
General trend of the stability of Si2O7 unit is related to the morphology of Si–OB–Si bridge and deserves to be discussed further. On the one hand, the strong electrostatic force between adjacent SiO4 tetrahedra in [Si2O7]6− groups tends to repulse each other apart and opens the Si–OB–Si angle toward 180°. Simultaneously, a π-type bonding forms between the bridging OB atom and Si atom.Reference Cruickshank27 The π-type interaction within Si2O7 unit will also stabilize the linear [O3Si]–OB–[SiO3] structure by modulating the bridging Si–OB bond length close to the values of the outer Si–O bonds. As a result, the linear [O3Si]–OB–[SiO3] configuration is energetically preferred in β and γ phases. In δ phase, however, stretching of Si–OB bond at high temperature results in the nonlinear configuration of Si–OB–Si bridge, and the inner Si–OB bonds are obviously longer than the outer Si–O bonds. Therefore, it is expected that the bended [O3Si]–OB–[SiO3] structure has lower stability and strength than the linear one.
Degrees of distortion ∆d of Si–O tetrahedra and Y–O polyhedra in β-, γ-, and δ-Y2Si2O7 are calculated and summarized in Table II. The average lengths of Si–O bond and Y–O bond in β-Y2Si2O7 are 1.610 and 2.237 Å, respectively, which are quite close to those in γ-Y2Si2O7 (1.608 Å for Si–O bond and 2.240 Å for Y–O bond). The two phases also have close values of ∆d for the SiO4 tetrahedra or YO6 polyhedra. For δ-Y2Si2O7, the average lengths of the two types of Si–O bonds differ only 0.3% and both are quite close to those in β- and γ-Y2Si2O7. The noticeable feature is that the values of ∆d of the two SiO4 tetrahedra in δ-Y2Si2O7 are much larger than those in β- and γ-Y2Si2O7. This result corresponds to the obviously different bond lengths for Si–OB and outer Si–O bonds. The average bond length of Y–O bond in YO7 polyhedra has higher values in δ-Y2Si2O7. Not surprisingly, the two types of YO7 polyhedra in δ-Y2Si2O7 have notably high values of ∆d, mainly because both the bridging OB atoms and the terminal O atoms participate in the first coordination shell of Y atoms and hence enhance the distortion of YO7 polyhedra.
Calculated bond lengths (in Å) and degree of distortion of SiO4 tetrahedra and YO6 or YO7 polyhedra in β-, γ-, and δ-Y2Si2O7 (CN denotes the corresponding coordination number).
Mulliken population, based on calculations of overlap populations for the nearest neighbors in the crystal, is a vivid parameter to display bonding strength. Larger Mulliken population suggests a higher strength of the covalent bonding.Reference Segall, Shah, Pickard and Payne28Figure 2 shows the Mulliken populations and bond lengths for various Y–O and Si–O bonds in the polymorphs. It is clear that Y–O bonds are considerably weaker than Si–O bonds in all of β-, γ-, and δ-Y2Si2O7. This result is also found in previous works by studying bond length and energy level of bonding states.Reference Wang, Zhou and Lin4,Reference Ching, Quyang and Xu29–Reference Sun, Zhou, Wang and Li31 In addition, both bond lengths and Mulliken populations of Y–O bonds are more dispersive than those of Si–O bonds, which suggests the higher levels of bonding heterogeneity in Y–O polyhedra. We also noted that the corresponding values of Y–O or Si–O bonds are quite close in β- and γ-Y2Si2O7, indicating similar bonding features in the two phases. It is interesting to find that either Si–O or Y–O bond contains particularly weaker one (Si–OB or Y–OB bond) in δ-Y2Si2O7 that drags down the average bonding strength of SiO4 or YO7 unit. The bridging Si–OB bond has higher bond length and lower Mulliken population than the outer Si–O bonds in δ-Y2Si2O7. Therefore, compared with the linear [O3Si]–OB–[SiO3] structure in β and γ phases, the nonlinear one in δ-Y2Si2O7 is weaker in bonding strength, which may obviously affect elastic and thermal properties of this phase. We will discuss in the following sections how the structural characteristics relate to the elastic properties and thermal conductivities of Y2Si2O7 polymorphs.
Mulliken populations of Y–O and Si–O bonds in β-, γ-, and δ-Y2Si2O7.
B. Elastic properties
Table III lists the full set of second-order elastic coefficients of β-, γ-, and δ-Y2Si2O7. Material with monoclinic crystal structure has thirteen independent elastic coefficients (c ij), and orthorhombic lattice has nine. From the second-order elastic coefficients, elastic parameters including bulk modulus B, shear modulus G, Poisson's ratio ν, and Young's modulus E are derived and listed in Table IV. All calculated elastic coefficients and elastic parameters are very close for β- and γ-Y2Si2O7, with differences less than 5%. The similar elastic stiffness could reasonably be attributed to their close crystal structures. In addition, the maximum value of elastic coefficients is 303 GPa (c 11) for β-Y2Si2O7 and 288 GPa (c 22) for γ-Y2Si2O7, both are nearly along the crystal directions of the [O3Si]–OB–[SiO3] units. The bulk modulus of δ-Y2Si2O7 is 17% lower than that of γ-Y2Si2O7, which is originated from the lower average bonding strength in the δ phase. Similar trend has also been reported for amorphous SiO2 for which the elastic modulus abnormally increases at elevated temperature by opening of the Si–OB–Si bond angle.Reference Liu, Wang, Zhou and Li32 The shear moduli of the three phases, on the other hand, do not vary significantly. In fact, the shear deformation resistance of Y2Si2O7 is dominated by the distortion of soft Y–O polyhedral, whereas Si2O7 pyrosilicate unit is rigid and remains stable under shearing.Reference Wang, Zhou and Lin4
Calculated second-order elastic coefficients (in GPa) of β-, γ-, and δ-Y2Si2O7.
Calculated bulk modulus B, shear modulus G, G/B ratio, Poisson's ratio ν, and Young's modulus E of β-, γ-, and δ-Y2Si2O7.
It is widely accepted that damage tolerance or quasi-plasticity of a ceramic is characterized by its capability to achieve shear deformation. According to the well-known Pugh's criterion, ductile and brittle materials are typically screened by the shear to bulk modulus ratio G/B of 0.57.Reference Yang, Guo, Linghu and Yang33 Calculated G/B ratios of β- and γ-Y2Si2O7 are 0.43 and 0.45, respectively; but the value of δ-Y2Si2O7 is 0.56, quite approaching to the critical Pugh's criterion. Besides, experimentally measured Vicker's hardness is 6.4 GPaReference Boakye, Keller, Mogilevsky, Parthasarathy, Ahrens, Hay and Cinibulk34 for β-Y2Si2O7 and 6.2 GPaReference Sun, Zhou, Wang and Li31 for γ-Y2Si2O7, which are comparable to the typical damage tolerant ceramic LaPO4 (H V = 4.68 GPaReference Sun, Zhou, Wang and Li31,Reference Wang, Pan, Chen, Fang, Cao and Luo35 ). Another parameter to evaluate quasi-plasticity is the ratio of H V/E,Reference Lawn36 which are 0.037 and 0.035 for β and γ phases, respectively. Their magnitudes are comparable to 0.036 for LaPO4.Reference Sun, Zhou, Wang and Li31 Therefore, β- and γ-Y2Si2O7 will show low shear deformation resistance, damage tolerance, and quasi-plasticity; while δ-Y2Si2O7 should be classified as a more brittle phase.
From the characteristics of crystal structures and second-order elastic coefficients, it is expected that β-, γ-, and δ-Y2Si2O7 might have anisotropic elastic modulus. We further estimate the direction dependent Young's modulus of the polymorphs by:Reference Nye37
$$\eqalign{ & {1 \over E} = \mathop {{l_1}}\nolimits^4 {s_{11}} + 2\mathop {{l_1}}\nolimits^2 \mathop {{l_2}}\nolimits^2 {s_{12}} + 2\mathop {{l_1}}\nolimits^2 \mathop {{l_3}}\nolimits^2 {s_{13}} + 2\mathop {{l_1}}\nolimits^3 {l_3}{s_{15}} + \mathop {{l_2}}\nolimits^4 {s_{22}} + 2\mathop {{l_2}}\nolimits^2 \mathop {{l_3}}\nolimits^2 {s_{23}} + 2{l_1}\mathop {{l_2}}\nolimits^2 {l_3}{s_{25}} \cr & \quad\quad\quad\quad\quad\!\!\! + \mathop {{l_3}}\nolimits^4 {s_{33}} + 2{l_1}\mathop {{l_3}}\nolimits^3 {s_{35}} + \mathop {{l_2}}\nolimits^2 \mathop {{l_3}}\nolimits^2 {s_{44}} + 2{l_1}\mathop {{l_2}}\nolimits^2 {l_3}{s_{46}} + \mathop {{l_1}}\nolimits^2 \mathop {{l_3}}\nolimits^2 {s_{55}} + \mathop {{l_1}}\nolimits^2 \mathop {{l_2}}\nolimits^2 {s_{66}}\quad {\rm{,}} \cr}$$where s ij is the elastic compliance constant, l 1, l 2, and l 3 are the direction cosines in any arbitrary direction. For orthorhombic lattice, s 15, s 25, s 35, and s 46 should be omitted. Figures 3(a)–3(c) present the topologies of the Young's modulus of β-, γ-, and δ-Y2Si2O7, respectively, along various directions. Calculated maximum and minimum Young's moduli are 309 and 105 GPa, respectively, for β-Y2Si2O7; 276 and 113 GPa for γ-Y2Si2O7; and 243 and 109 GPa for δ-Y2Si2O7. The ratios of E max/E min, which are 2.94 for β-Y2Si2O7, 2.44 for γ-Y2Si2O7, and 2.23 for δ-Y2Si2O7, are taken to estimate the anisotropic characters of Young's moduli along different crystallographic orientations. Anisotropic elastic stiffness of a material typically relates to arrangements of chemical bonding.Reference Liu, Wang, Zhou, Liao and Li38 For β-Y2Si2O7, calculated maximum Young's modulus lies in the direction very close to the orientation of linear [O3Si]–OB–[SiO3] units (the directional divergence ∆φ ≤ 9°). And the case is similar for γ-Y2Si2O7, although there are two equivalent orientations due to crystallographic symmetry. Detailed analysis of the anisotropic elastic modulus of γ-Y2Si2O7 could be found in our previous work.Reference Luo, Wang, Wang, Li and Hu13 As discussed in the previous section, Si2O7 unit is more rigid than YO6 polyhedron, so β- and γ-Y2Si2O7 exhibit higher resistance along the linear [O3Si]–OB–[SiO3] structure. For δ-Y2Si2O7, the maximum value of Young's modulus lies in two equivalent directions according to crystallographic symmetry, but none of them is relevant to the nonlinear [O3Si]–OB–[SiO3] units. Therefore, the bended Si–OB–Si bridge is far less effective than the linear Si–OB–Si bridge to provide high elastic stiffness.
Illustrations of directional dependent Young's moduli of (a) β-, (b) γ-, and (c) δ-Y2Si2O7.
C. Thermal conductivity
Table V presents the calculated parameters including sound velocity, Debye temperature, theoretical minimum thermal conductivity, and Grüneisen constant of β-, γ-, and δ-Y2Si2O7. The average sound velocities of the three polymorphs are close in magnitudes, with the maximum difference within 1.11%. As is indicated in Eqs. (8) and (11), similar values of average sound velocity νm and average interatomic distance δ yield close values of Debye temperature and minimum thermal conductivity. The calculated κmin values are 1.14, 1.16, and 1.18 W·m−1·K−1 for β-, γ-, and δ-Y2Si2O7, respectively, which suggest them as the candidates with low thermal conductivity.
Calculated sound velocity (v s, v l, and v m), Debye temperature θD, theoretical minimum thermal conductivity κmin, and Grüneisen constant γ of β-, γ-, and δ-Y2Si2O7.
Anharmonic phonon interactions are the origin of thermal resistance and its intensity could be measured by the Grüneisen constant γ. As is shown in Table V, the calculated acoustic Grüneisen constant is 1.84, 1.77, and 1.60 for β-, γ-, and δ-Y2Si2O7, respectively. Specifically, the Grüneisen constant γ of the δ phase is obviously smaller than the values of other two phases, which could be attributed to its distinct crystal structure. According to Eqs. (6) and (13), it is interesting to link the G/B ratio with Grüneisen constant γ as:
$$\gamma = {3 \over 2}{{\left( {1 + {{\left( {3 - 2{G \over B}} \right)} \mathord{\left/ {\vphantom {{\left( {3 - 2{G \over B}} \right)} {2\left( {3 + {G \over B}} \right)}}} \right. \kern-\nulldelimiterspace} {2\left( {3 + {G \over B}} \right)}}} \right)} \mathord{\left/ {\vphantom {{\left( {1 + {{\left( {3 - 2{G \over B}} \right)} \mathord{\left/ {\vphantom {{\left( {3 - 2{G \over B}} \right)} {2\left( {3 + {G \over B}} \right)}}} \right. \kern-\nulldelimiterspace} {2\left( {3 + {G \over B}} \right)}}} \right)} {\left( {{{2 - 3\left( {3 - 2{G \over B}} \right)} \mathord{\left/ {\vphantom {{2 - 3\left( {3 - 2{G \over B}} \right)} {2\left( {3 + {G \over B}} \right)}}} \right. \kern-\nulldelimiterspace} {2\left( {3 + {G \over B}} \right)}}} \right)}}} \right. \kern-\nulldelimiterspace} {\left( {{{2 - 3\left( {3 - 2{G \over B}} \right)} \mathord{\left/ {\vphantom {{2 - 3\left( {3 - 2{G \over B}} \right)} {2\left( {3 + {G \over B}} \right)}}} \right. \kern-\nulldelimiterspace} {2\left( {3 + {G \over B}} \right)}}} \right)}}\quad {\rm{.}}$$This equation builds a correlation between quasi-plasticity and the level of anharmonicity of phonon interaction. Increment of quasi-plasticity (lower G/B ratio) will lead to the higher level of anharmonicity (larger γ). Since lattice thermal conductivity is inversely proportional to the value of γ2, Y2Si2O7 polymorphs with damage tolerance and quasi-plasticity will also have low lattice thermal conductivity.
In our previous study, Slack's equation had been used to predict the temperature dependence of lattice thermal conductivities of γ-Y2Si2O7 and Y2SiO5.Reference Luo, Wang, Wang, Li and Hu13 In this work, we follow the same method to compare the κ–T relationships of β-, γ-, and δ-Y2Si2O7 over a wide range of temperature. Results are depicted in Fig. 4, together with the theoretical minimum thermal conductivity at high temperature limit (marked with dashed-dotted line). We only plot the line at κmin = 1.16 W·m−1·K−1 for the three phases, because their deviations are too small to be distinguished in the figure. Above the phase transition temperature of each phase, the theoretical lattice thermal conductivity is presented by dashed line and is a guide for the eyes. Generally, the thermal conductivities of β-, γ-, and δ-Y2Si2O7 decrease as ∼1/T from 500 K (close to Debye temperatures) to high temperatures when Umklapp phonon scattering is the dominant mechanism.Reference Clarke18 The κ–T relationship yields κ = 2485/T for β-Y2Si2O7, κ = 1761/T for γ-Y2Si2O7, and κ = 1406/T for δ-Y2Si2O7. Values of κ(T) below κ min have no physical meaning since κ min bounds the lowest value for lattice thermal conductivity.
Temperature dependences of thermal conductivities of β-, γ-, and δ-Y2Si2O7.
To better elucidate the origin of the different κ–T relationships of Y2Si2O7 polymorphs, we factorize the Slack's equation into three parts: Structural factor 
${{\overline M \delta } \mathord{\left/ {\vphantom {{\overline M \delta } n}} \right. \kern-\nulldelimiterspace} n}^{{2 \over 3}}$ characterizing crystal structure characteristics, thermal parameter θDReference Lee, Fox and Bansal3 derived from elastic constants, and anharmonic coefficient A(γ)/γ2 representing anharmonicity of phonon scattering. The lattice thermal conductivity is determined by the balance of the three parts. Since the Debye temperatures of the three phases have close values, we mainly discuss the contributions of the other two factors. For β- and γ-Y2Si2O7, the difference of A(γ)/γ2 is only 8.12%, but that of 
${{\overline M \delta } \mathord{\left/ {\vphantom {{\overline M \delta } n}} \right. \kern-\nulldelimiterspace} n}^{{2 \over 3}}$ is as high as 59.32%. As analyzed above, β- and γ-Y2Si2O7 have similar crystal structures and bonding characteristics. Both phases have the same values of mean atomic mass 
$\bar M$ and similar values of interatomic distance δ. Thus, the discrepancy (41.11%) between thermal conductivities of β- and γ-Y2Si2O7 is mainly caused by parameter n, the number of atoms per primitive unit cell, which is 11 for β-Y2Si2O7 and 22 for γ-Y2Si2O7. For δ-Y2Si2O7, 
${{\overline M \delta } \mathord{\left/ {\vphantom {{\overline M \delta } n}} \right. \kern-\nulldelimiterspace} n}^{{2 \over 3}}$ is 37.29% lower than that of γ-Y2Si2O7 because the δ phase has 44 atoms per primitive unit cell; but coefficient A(γ)/γ2 is 24.70% higher than that of γ-Y2Si2O7 owing to the lower Grüneisen constant γ. Hence, the lattice thermal conductivity of δ-Y2Si2O7 is lower than that of γ-Y2Si2O7 when the contribution of structural difference is partly mediated by that of anharmonicity.
It is worthy of discussing the physical picture how the number of atoms in primitive cell modifies the lattice thermal conductivities of the polymorphs. Number of atoms per primitive cell, n, corresponds to the complexity of crystal structure and plays an important role in determining the maximum acoustic mode frequency. With larger n, the lattice parameters increase in real space and the position of Brillouin zone boundary in reciprocal space moves inward and cuts the acoustic phonon branch with lower frequency. In addition, the phonon dispersion folds in on itself as n increases. These two physical mechanisms will decrease the group velocity of acoustic mode phonons, as well as the acoustic Debye temperature θa (approximated by θa = θDn −1/3, where θD is the traditional Debye temperature derived from elastic modulus). Therefore, the thermal conductivity will be obviously affected when n changes in the primitive cell.
As for methodology of evaluating thermal conductivity among polymorphs, Slack's equation by itself provides important information for materials selection. For a group of polymorphic ceramics, wherein, the mean atomic mass is mostly the same and the average interatomic distance is quite similar, the lattice thermal conductivity predominantly depends on anharmonicity of phonon scattering, Debye temperature, and the number of atoms per primitive unit cell. That is to say, higher levels of structural complexity and anharmonicity, which are characterized by higher values of n and γ, will lead to the lower thermal conductivity. Based on the Debye–Slack model, one can conveniently investigate what specific factors govern overall lattice thermal conductivity and select candidates with desired thermal conductivity among polymorphs.
IV. CONCLUSION
In summary, the crystal structure, elastic stiffness, and temperature dependence of lattice thermal conductivities of β-, γ-, and δ-Y2Si2O7 are predicted and compared based on first-principle pseudopotential total energy calculations. We comprehensively investigate the full set of second-order elastic coefficients, bulk modulus, shear modulus, Young's modulus, Poisson's ratio, sound velocity, Debye temperature, theoretical minimum thermal conductivity, and Grüneisen constant. It is found that β-, γ-, and δ-Y2Si2O7 display anisotropic elastic moduli and the trends are explained in terms of structural and bonding characteristics. The results further show that the nonlinear [O3Si]–OB–[SiO3] unit (Si–OB–Si bond angle is 158.1°) in δ-Y2Si2O7 lowers the strength and stability but enhance the heterogeneity of chemical bonding. Besides, calculated parameters are used to predict the minimum lattice thermal conductivity at high temperature limit and the κ–T relationship using the Debye–Slack model. Temperature dependences of lattice thermal conductivities are: κ = 2485/T for β-Y2Si2O7, κ = 1761/T for γ-Y2Si2O7, and κ = 1406/T for δ-Y2Si2O7. At high temperature limit, the lattice thermal conductivities are bounded by minimum values around 1.14–1.18 W·m−1·K−1. We also propose that the divergences of lattice thermal conductivities of Y2Si2O7 polymorphs mainly stem from the different levels of anharmonicity of phonon scattering and complexity of crystal structure. Finally, this work uses integrated method to investigate the specific factors dominating lattice thermal conductivity, which could promisingly be generalized to search novel candidates with extremely low lattice thermal conductivity.
ACKNOWLEDGMENTS
This work was supported by the Natural Sciences Foundation of China under Grant Nos. 51032006 and 51372252.
