I. INTRODUCTION
Zircon (ZrSiO4) and hafnon (HfSiO4) have attracted considerable interests owing to their attractive properties, such as high chemical stability and ability to accommodate a number of dopant ions, low thermal expansion, and high resistance to thermal shock.Reference Veytizou, Quinson and Jorand1 They are considered as promising candidates for applications in nuclear industry, thermal barrier coating, and solid state laser.Reference Ewing, Lutz and Weber2–Reference Yang, Luff and Townsend4
Inspired by the promising potential applications, many efforts have been devoted to the investigations on their structure, phase stability, and mechanical properties. The crystal of MSiO4 (M = Zr, Hf) is body-centered tetragonal with an I41/amd space group. The conventional unit cell of zircon structure contains four formula units (Z = 4). The crystal structure of MSiO4 is shown in Fig. 1. It consists of alternating (discrete) SiO4 tetrahedra and MO8 dodecahedra, sharing edges to form chains parallel to the c direction. The structure of MO8 dodecahedra is distorted, with four O atoms closer to the M atom than the four other ones. The position of M and Si atoms are located at (0, 3/4, 1/8) and (0, 1/4, 3/8) on the 4a and 4b Wyckoff sites, respectively. The O atoms occupy the 16h Wyckoff sites (0, u, v), where u and v are the internal parameters. The characteristic of crystal of MSiO4 has been fully understood experimentally and theoretically.Reference Crocombette and Ghaleb5–Reference Liu, Liu, Feng, Tian and Zeng16 The high pressure phase transition of zircon is also a hot topic since the zircon–reidite transition is fast and occurs at ambient temperature, which is quite unusual in silicate minerals.Reference van Westrenen, Frank, Hanchar, Fei, Finch and Zha17 Experimentally, the transition pressure ranged from 10 to 23 GPa for room temperature transitions, around 12 for high temperature transitions.Reference van Westrenen, Frank, Hanchar, Fei, Finch and Zha17–Reference Knittle and Williams19 Theoretically, Du et al. estimated a transition pressure of 11.6 GPa.Reference Du, Devanathan, Corrales and Weber11 The mechanical properties of MSiO4, which are important to the applications as refractory materials, are also studied. It was found that the bending strength of zircon ceramics could maintain up to 1000–1400 °C, depending on the purity of sintered body.Reference Shi, Huang and Yan20,Reference Mori, Yamamura, Kobayashi and Mitamura21 Moreover, recently reported experiments revealed that MSiO4 are the main constitutions of oxidation products of Zr- and Hf-based ultrahigh temperature ceramics, especially under low oxygen partial pressure and simulated hypersonic flight conditions.Reference Gao, Zhang, Fu, Xu, Song and Shi22,Reference Parthasarathy, Petry, Cinibulk, Mathur and Gruber23 The presence of high content of zircon phase could improve the oxidation resistance of the matrix.
(a) Crystal structure of MSiO4. (b) Viewed from (001) direction. The yellow polyhedra are SiO4 tetrahedra.
Despite the structure, phase stability and mechanical properties of zircon phase have been reported, the structure–property relationship, anisotropy in mechanical properties and thermal properties of MSiO4 (M = Zr, Hf) are still not fully understood. The thermal properties, especially the thermal conductivity–temperature relationship, are crucial to the applications as thermal barrier coating and solid state laser. In addition, mechanical anisotropy and the ability of withstanding external mechanical damages are the key factors need to be considered when using as structural components. Lacking of this vital knowledge impedes the investigations and applications of MSiO4.
The purpose of this contribution is to investigate the structure–property relationship, anisotropy in mechanical properties and thermal properties of MSiO4 (M = Zr, Hf) by theoretical simulations. The simulation details are given in the next section, which are followed by the results and discussions on the equilibrium structure, mechanical, and thermal properties of MSiO4 from density functional theory (DFT) calculations and chemical bond theory (CBT).
II. SIMULATION DETAILS
The Cambridge Serial Total Energy Package (CASTEP) codeReference Segall, Lindan, Probert, Pickard, Hasnip, Clark and Payne24 was used in the DFT calculations on structural and elastic properties of MSiO4. The Vanderbilt-type ultra-soft pseudopotentialsReference Vanderbilt25 were used to represent the interactions between the ionic cores and the valence electrons. The exchange-correlation energy was treated under the local density approximation (LDA).Reference Perdew, Burke and Ernzerhof26 The plane-wave basis set cut-off energy was fixed at 400 eV in all calculations after convergence test. The special points sampling integration over Brillouin zone was conducted using the Monkhorst–Pack method with a 6 × 6 × 6 special k-points mesh.Reference Monkhorst and Pack27 The population analysis of the bond is performed using a projection of the plane wave states onto a localized basis using a technique described by Sanchez-Portal et al.,Reference Sanchez-Portal, Artacho and Soler28 and the resulting projected states is then performed using the formalism described by Segall et al.Reference Segall, Shah, Pickard and Payne29
The geometry optimization was achieved under the Broyden–Fletcher–Goldfarb–Shanno minimization scheme.Reference Pfrommer, Côté, Louie and Cohen30 The crystal structures, including lattice constants and internal atomic coordinates, were optimized independently to minimize the free enthalpy, interatomic forces and stresses of the unit cell. The tolerances for geometrical optimization were: differences for total energy within 5 × 10−6 eV/atom, maximum ionic Hellmann–Feynman force within 0.01 eV/Å, maximum ionic displacement within 5 × 10−4 Å, and maximum stress within 0.02 GPa.
The theoretical elastic coefficients were determined from the first-principles strain–stress relationships method implemented by Milman and Warren.Reference Milman and Warren31 Both positive and negative strains were applied for each strain component, with a maximum strain value of 0.3%. The criteria for convergence of optimization on atomic internal freedoms were selected as differences in the total energy within 1 × 10−6 eV/atom, Hellmann–Feynman force within 0.002 eV/Å, and maximum ionic displacement within 1 × 10−4 Å. The bulk modulus B and shear modulus G were calculated from the compliance tensor by Voigt–Reuss–Hill approximation.Reference Voigt32–Reference Hill34 The Young's modulus E H and Poisson's ratio ν were calculated using Hill's bulk modulus B H and shear modulus G H using the following equationsReference Green35:
$$E = {{9{B_{\rm{H}}}{G_{\rm{H}}}} \over {3{B_{\rm{H}}} + {G_{\rm{H}}}}}\quad ,$$
$$\nu = {{3{B_{\rm{H}}} - 2{G_{\rm{H}}}} \over {2\left( {3{B_{\rm{H}}} + {G_{\rm{H}}}} \right)}}\quad .$$The minimum thermal conductivity (κmin) of MSiO4 was calculated based on the model proposed by Clarke,Reference Clarke36 who suggested that the theoretical minimum conductivity was calculated after replacing different atoms by an equivalent atom with a mean atomic mass of M/n:
$$\kappa _{{\rm{min}}} \to k_{\rm{B}} \upsilon _{\rm{m}} \left( {{M \over ndN_{{\rm{Av}}} }}} \right)^{ - {2 \mathord{\left/ {\vphantom {2 3}} \right. \kern-\nulldelimiterspace} 3}} \quad ,$$where k B is the Boltzmann's constant, υm the average sound velocity, N Av the Avogadro's number, d the density, M the molecular weight, and n the number of atoms in the molecule. Clarke used a simplification for υm
$\left( {{\upsilon _{\rm{m}}} = 0.87\sqrt {{E \over d}} } \right)$ without considering anisotropy. Liu et al.Reference Liu, Wang, Li and Zhou37 proposed a modification to the Clarke's relation by taking anisotropy of sound velocity into account. This modified model has been successfully applied to estimate the κmin of several complex materials, like Yb3Al5O12,Reference Zhou, Xiang and Feng38 YbAlO3,Reference Xiang, Feng and Zhou39 and Yb2SiO5.Reference Xiang, Feng and Zhou40
The average sound velocity υm of polycrystalline material is defined as the followingReference Anderson41:
$$\upsilon _{\rm{m}} = \left[ {{1 \over 3}\left( {{2 \over {\upsilon _{\rm{s}}^3 }} + {1 \over {\upsilon _{\rm{l}}^3 }}} \right)} \right]^{ - {1 \mathord{\left/ {\vphantom {1 3}} \right. \kern-\nulldelimiterspace} 3}} \quad ,$$where 
$\upsilon _{\rm{l}} = \sqrt {{{{{B + 4G} \mathord{\left/ {\vphantom {{B + 4G} 3}} \right. \kern-\nulldelimiterspace} 3}} \over d}}$ and 
${\upsilon _{\rm{s}}} = \sqrt {{G \over d}}$ are longitudinal and transverse sound velocities, respectively.
The temperature dependent thermal conductivity of MSiO4 was calculated from Slack's modelReference Slack42:
$${\kappa _{\rm{L}}} = A{{\bar M{\Theta ^3}\delta } \over {{\gamma ^2}{n^{{2 \mathord{\left/ {\vphantom {2 3}} \right. \kern-\nulldelimiterspace} 3}}}T}}\quad ,$$where A is a physical constants (A ≈ 3.1 × 10−6 if κL is in W m−1 K−1, and δ in Å), 
$\bar M$ is the average mass of the atoms in the crystal, δ3 is the volume per atom, γ is the high temperature limit of phonon mode Grüneisen parameter, which could be derived from Poisson's ratio νReference Sanditov, Tsydypov and Sanditov43:
$$\gamma = {{9\left( {\upsilon _{\rm{l}}^2 - {4 \over 3}\upsilon _{\rm{s}}^2} \right)} \over {2\left( {\upsilon _{\rm{l}}^2 - 2\upsilon _{\rm{s}}^2} \right)}} = {3 \over 2}\left( {{{1 + \nu } \over {2 - 3\nu }}} \right)\quad ,$$and ΘD is the Debye temperature, which is defined as:
$${\Theta _{\rm{D}}} = {h \over {{k_{\rm{B}}}}}{\left[ {{{3n} \over {4\pi }}\left( {{{{N_{{\rm{Av}}}}d} \over M}} \right)} \right]^{{1 \mathord{\left/ {\vphantom {1 3}} \right. \kern-\nulldelimiterspace} 3}}}{\upsilon _{\rm{m}}}\quad ,$$where h is the Plank's constant, k B is the Boltzmann's constant, and N Av is the Avogadro's number.
III. RESULTS AND DISCUSSION
A. Equilibrium structure
Different electronic exchange-correlation interaction approximation methods [i.e., generalized-gradient approximation (GGA) and LDA] are used in the calculation of equilibrium lattice configuration of MSiO4. Comparing to the experimental values,Reference Robinson, Gibbs and Ribbe7 the lattice parameters calculated by GGA deviate by 3.5% while those by LDA deviate by 1.7%. Therefore, LDA is implemented in the following calculations for its higher accuracy. Table I summarized the theoretical structural parameters from previousReference Robinson, Gibbs and Ribbe7–Reference Liu, Liu, Feng, Tian and Zeng16 and this work. The calculated lattice constants a and c, as well as the internal parameters u and v are found to be in excellent agreement with the previous theoretical and experimental values, the lattice parameters calculated in the present work deviate from experiments by 1.1% for ZrSiO4 and 1.7% for HfSiO4. This ensures the accuracy of the further study on the mechanical and thermal properties of these two silicates.
Theoretical and experimental lattice parameters of ZrSiO4 and HfSiO4.
f Ref. Reference Speer and Cooper12.
The Mulliken population analyses on M–O and Si–O bonds of MSiO4 are conducted in this work. As tabulated in Table I, for both compounds, the length and population of Si–O bonds are shorter and smaller than those of M–O bonds, showing stronger and more covalent bonding properties of Si–O bonds. This feature can be confirmed by calculating the bond covalency and energy from CBT (for detail description of CBT, one may refer to Refs. Reference Zhou, Xiang and Feng38–Reference Xiang, Feng and Zhou40). The main concept of CBT is the decomposition of the chemical formula and estimating the properties of subformulas. For MSiO4 zircon phase, the chemical formula can be decomposed into the following subformulas:
$${\rm{MSi}}{{\rm{O}}_4} = {{\rm{M}}_{{1 \mathord{\left/ {\vphantom {1 2}} \right. \kern-\nulldelimiterspace} 2}}}{{\rm{O}}_{{4 \mathord{\left/ {\vphantom {4 3}} \right. \kern-\nulldelimiterspace} 3}}} + {{\rm{M}}_{{1 \mathord{\left/ {\vphantom {1 2}} \right. \kern-\nulldelimiterspace} 2}}}{\rm{O}}{\prime _{{4 \mathord{\left/ {\vphantom {4 3}} \right. \kern-\nulldelimiterspace} 3}}} + {\rm{Si}}{{\rm{O}}_{{4 \mathord{\left/ {\vphantom {4 3}} \right. \kern-\nulldelimiterspace} 3}}}\quad .$$ The superscript in the formula marks the oxygen atoms with different bond lengths in MO8 dodecahedra. The chemical parameters of each subformula are calculated and listed in Table II. The previously reported values are also included in the table. Our results are quite close to those reported by Li et al.,Reference Li, Zhou and Zhang44 with the largest difference of 1%, as shown in Table II. For all these bonds, the heteropolar part of energies (C μ) are larger than the homopolar part of energies 
$\left( {E_{\rm{h}}^\mu } \right)$, which indicates complex and more ionic bonding nature of these bonds. The calculated covalency 
$\left( {f_{\rm{c}}^\mu } \right)$ of Si–O bonds is larger than those of M–O bonds, which is consistent with the Mulliken population analysis on the bonding characteristic of MSiO4. In addition, the bond energies (U μ) of Si–O bonds are larger than those of M–O bonds, showing strong bonding nature of Si–O bonds compared with M–O bonds. The bond energies of Hf–O bonds are slightly larger than those of Zr–O bonds, since the radius of 8-coordinated Hf ion (71 pm) is slightly smaller than that of Zr (72 pm).Reference Carter and Norton45
Chemical parameters of subformulas for MSiO4 by CBT. C μ and 
$E_{\rm{h}}^\mu$ are in eV, and U μ is in kJ/mol.
In general, MSiO4 shows a heterogeneous bonding nature. SiO4 tetrahedra are strongly constructed and show rigid characteristic while MO8 dodecahedra are relatively weak and soft. The relationship between these structural features and the elastic and thermal properties of MSiO4 will be discussed in the next sections.
B. Elastic properties
To investigate the elastic properties of MSiO4, the full set of second-order elastic coefficients are calculated first. The calculated elastic constants are listed in Table III, accompanied by previously reported onesReference Liu, Liu, Feng, Tian and Zeng16,Reference Özkan and Jamieson46,Reference Chaplot, Mittal and Choudhury47 as comparison. Since the lattice of ZrSiO4 is smaller than the experimental one in our theoretical simulations, the present calculated elastic constants are overestimated by 5%. With regard to HfSiO4, our theoretical values are close to the values reported by Mittal et al.Reference Chaplot, Mittal and Choudhury47 while deviate from the values reported by Liu et al.Reference Liu, Liu, Feng, Tian and Zeng16 The discrepancy may originate from the different electronic exchange-correlation interaction approximation methods. Since the deviation of lattice parameters is smaller, the results calculated from LDA in current work are more reliable. However, since the lattice parameters of HfSiO4 are overestimated in our simulation, the experimental values of HfSiO4 are expected to be lower than our theoretical predictions. The discrepancy between different theoretical data shows the lack of knowledge on the elastic properties of HfSiO4 and illustrates that more delicate experimental and theoretical investigations are needed.
Calculated independent second-order elastic constants c ij (in GPa) for ZrSiO4 and HfSiO4.
For both materials, the stiffness against principal strains is much higher than that against shear deformations, indicating less resistance to the shear deformations. For instance, c 66 yields only 49 and 43 GPa for ZrSiO4 and HfSiO4, respectively. It is not only much smaller than c 11 and c 33, but also smaller than c 44. It implies that it would be easier to perform shear deformation along [010] direction on (100) planes. The anisotropy of Young's modulus and shear modulus will be discussed in the next section.
The mechanical stability of a crystal structure can be justified from the Born–Huang criterion.Reference Born and Huang48 The full expressions of the criterion for a tetragonal structure can be expressed as:
$$\left( {{c_{11}} - {c_{12}}} \right) > 0,\;\left( {{c_{11}} + {c_{33}} - 2{c_{13}}} \right) > 0,\;{c_{11}} > 0,\;{c_{33}} > 0,\;{c_{44}} > 0,\;{c_{66}} > 0\;{\rm{and}}\;\left( {2{c_{11}} + {c_{33}} + 2{c_{12}} + 4{c_{13}}} \right) > 0\quad .$$The criterions are fully satisfied by both silicates, demonstrating their stability under elastic strain perturbation.
The elastic parameters for polycrystalline aggregates, namely bulk modulus B, shear modulus G, Young's modulus E, and Poisson's ratio ν, are estimated by Voigt-Reuss-Hill (VRH) approximation and tabulated in Table IV. The consistence shown in Table IV between our results and the existing theoretical and experimental values shows the reliability of our estimations. In this work, we also use an empirical method based on CBT to estimate the bulk modulus of these materials and validate our calculations. For detail procedure, one may refer to Refs. Reference Zhou, Xiang and Feng38–Reference Xiang, Feng and Zhou40. The chemical parameters of these two silicates are listed in Table V, B μ is the contribution of one particular bond to the bulk modulus B. The bulk moduli of MSiO4 estimated from CBT are quite close to the ones calculated from first-principles. The consistence confirms the first-principle calculation and shows the suitability of CBT on estimation of bulk moduli of materials with the same structure as MSiO4. The moduli of Si–O bonds are obviously larger than those of M–O bonds, showing less compressibility and rigid characteristic of SiO4 tetrahedra.
Calculated bulk modulus B H (GPa), shear modulus G H (GPa), Young's modulus E (GPa), Poisson's ratio (ν) and Pugh's ratio (B/G) of ZrSiO4 and HfSiO4.
Estimation of bond volume νμ (Å3), lattice energy density u μ (GPa), bulk modulus B (GPa) and linear TEC α(×10−6 K−1) of MSiO4 from CBT.
a Ref. Reference Carter and Norton45.
The shear moduli G of both materials are substantially smaller than bulk moduli B, as shown in Table IV. Pugh introduced the quotient of bulk to shear modulus of polycrystalline phases (B/G), which reflects the ductility of a material.Reference Pugh49 A high (low) B/G value is associated with ductility (brittleness), and the critical value is about 1.75. The theoretical B/G value for ZrSiO4 and HfSiO4 from calculations is 2.13 and 2.28, respectively. The B/G values for MSiO4 are significantly higher than 1.75, suggesting that MSiO4 are “quasi-ductile” materials. Another parameter to evaluate the plasticity is the ratio of H V3/E*2, where H V is the microhardness and E* = E/(1 − ν2) is the effective Young's modulus. As described by Musil et al.,Reference Musil, Kunc, Zeman and Polakova50 high (low) value of H V3/E*2 is associated with high (low) resistance to plastic deformation. For ZrSiO4, the value of H V3/E*2 is 0.010, while for a typical brittle ceramic, α-Al2O3, it is 0.0441. This value is 0.008 for the well-known mechinable ceramic LaPO4.Reference Wang, Zhou and Lin51 Therefore, the resistance to plastic deformation of MSiO4 is quite low and is considered as “quasi-ductile” materials.
The Poisson's ratio is associated with the volume change during uniaxial deformation.Reference Ravindran, Fast, Korzhavyi, Johansson, Wills and Eriksson52 If ν = 0.5, no volume change occurs during elastic deformation. Low Poisson's ratio for MSiO4 means that a large volume change is associated with deformation. Moreover, Poisson's ratio provides information about the characteristics of the bonding forces. For a typical metal, the value is supposed to be 0.33; for the ionic-covalent crystal, the value is situated between 0.2 and 0.3; the strong covalent crystal has even smaller Poisson's ratio. The calculated Poisson's ratio is 0.30 and 0.31 for ZrSiO4 and HfSiO4, respectively, indicating their ionic-covalent bonding nature.
The intrinsic hardness of MSiO4 is estimated by the model proposed by Chen et al.Reference Chen, Niu, Li and Li53 The empirical model proposed by Chen et al. correlates not only with shear modulus G, but also with bulk modulus B. In this model, the Vickers hardness is estimated by the following formula:
$${H_{\rm{V}}} = 2{\left( {{k^2}G} \right)^{0.585}} - 3\quad ,$$where k is the shear modulus to bulk modulus ratio. The estimated intrinsic hardness is 10 and 8.8 GPa for ZrSiO4 and HfSiO4, respectively. Experimentally, the hardness of ZrSiO4 is closely related to the sintering techniques. Shi et al. reported a hardness of 10 GPa for hot-pressed zircon ceramic.Reference Shi, Huang and Yan20 Rendtorff et al. fabricated dense zircon ceramic by high energy ball milling and spark plasma sintering, the hardness they given was 11–13 GPa.Reference Rendtorff, Grasso, Hu, Suarez, Aglietti and Sakka54 The hardness of ZrSiO4 estimated from this model is close to the experimental ones, implying suitability of this model in predicting the hardness of MSiO4. Though the hardness of HfSiO4 is smaller than that of ZrSiO4, relatively high hardness is adverse to the machinability of these two materials.
C. Elastic anisotropy
It is widely accepted that microcracks are induced in materials owing to the anisotropy of the coefficient of thermal expansion as well as elastic anisotropy. Thus, understanding the elastic anisotropy is essential to the application of a material. As listed in Table IV, the theoretical Young's moduli of MSiO4 are slightly different along different principle directions, implying MSiO4 are not isotropic in elastic properties. To describe the elastic anisotropy of a crystal, a more practical method is the spatial surface construction of mechanical moduli as a function of direction. The direction dependent Young's modulus (E) for tetragonal crystal is defined as:Reference Nye55
$${1 \over E} = {s_{11}}\left( {l_1^4 + l_2^4} \right) + \left( {2{s_{12}} + {s_{66}}} \right)l_1^2l_2^2 + {s_{33}}l_3^4 + \left( {2{s_{13}} + {s_{44}}} \right)l_3^2\left( {1 - l_{3}^2} \right)\quad ,$$where s ij is the elastic compliance constant, and l 1, l 2, and l 3 are the directional cosines of angles with the three principal directions, respectively.
The overall surface contours of Young's modulus of MSiO4 are shown in Fig. 2, and the planar projections on (100), (001), and (110) crystallographic planes are also shown. For different crystallographic planes, A and B represent different crystallographic directions. A represents [100], and B represents [010] for (001) plane; in the case of (100) plane, they are [100] and [001]; while for (110) plane, they are 
$\left[ {1\bar 10} \right]$ and [001] directions. The anisotropy of Young's moduli of MSiO4 is quite noticeable. The anisotropy on the (001) plane is much stronger than the other two planes. The minimum Young's moduli are parallel to [110] direction for both silicates, which are 160 and 144 GPa for ZrSiO4 and HfSiO4, respectively. Low Young's modulus is associated with low fracture energy, thus [110] direction is the weakest direction under uniaxial tension for both zircon phases. The maximum Young's modulus of ZrSiO4 is 399 GPa, which is 2.49 times larger than the minimum. For HfSiO4, the largest Young's modulus is 433 GPa. It is 3 times larger than the smallest one. The anisotropy of Young's modulus in HfSiO4 is severer than that in ZrSiO4.
(a, c) Surface contour of direction dependent Young's modulus of ZrSiO4 and HfSiO4, respectively. (b, d) Planar projections on (100), (001), and (110) crystallographic planes for ZrSiO4 and HfSiO4, respectively.
To get a clear view on the shear behavior of MSiO4 on different shear planes, it is practical to discuss the directionality of shear modulus. However, it is impossible for shear modulus to construct a spatial surface since there are two determining directions while a surface presents only one radius vector. It can only be presented by the polar plots superposed on a particular shear plane or direction. In this work, we choose three most likely shear planes [(100), (001), and (110)] to examine the shear anisotropy of MSiO4. The shear modulus is related to the compliances as follow:
$$\mathop {{G_{ij}}}\nolimits^\prime = {1 \over {4\mathop {{s_{ijij}}}\nolimits^\prime }}\quad ,$$where s ijij′ is the compliances under certain coordination system. The shear modulus on the considered planes can be expressed as:
$$\mathop {{G_{\left( {100} \right)}}}\nolimits^\prime = {1 \over {{s_{66}} + \left( {{s_{44}} - {s_{66}}} \right){{\sin }^2}\theta }}\quad ,$$
$$\mathop {{G_{\left( {001} \right)}}}\nolimits^\prime = {1 \over {{s_{44}}}}\quad ,$$
$$\mathop {{G_{\left( {110} \right)}}}\nolimits^\prime = {1 \over {{s_{44}} + \left( {2{s_{11}} - 2{s_{12}} - {s_{44}}} \right){{\cos }^2}\theta }}\quad .$$The deduction procedure has been discussed thoroughly in previous reported work.Reference Turley and Sines56 The shear modulus on different planes is plotted in Fig. 3. The anisotropy of shear modulus of ZrSiO4 is quite similar to that of HfSiO4, they both exhibit obvious shear anisotropy. The shear anisotropy on different planes is quite different from each other. To provide a measure of the degree of anisotropy in bonding between atoms on different planes, we define a shear anisotropic factor:
$$A = {{{G_{{\rm{max}}}}} \over {{G_{{\rm{min}}}}}}\quad .$$Directionality of shear modulus on different shear planes for (a) ZrSiO4 and (b) HfSiO4.
For (100) shear plane, 
${A_{(100)}} = {{{S_{66}}} \over {{S_{44}}}}$, and for (001) shear plane, A (001) = 1, while for (110) shear plane, 
${A_{\left( {110} \right)}} = {{{S_{44}}} \over {2\left( {{S_{11}} - {S_{12}}} \right)}}$.
For an isotropic crystal, the factors A (100), A (001), and A (110) must be unity, while any value that unequal to unity is a measure of the degree of elastic anisotropy. The anisotropic factors for (100), (001), (110) shear planes are 2.30, 1, 1.66 and 2.49, 1, 1.96 for ZrSiO4 and HfSiO4, respectively. We can find out that the shear anisotropy on (100) plane is much severer than the other two planes for both silicates, and the shear anisotropy of HfSiO4 is larger than that of ZrSiO4. The lowest shear moduli, which are only 49.4 and 42.8 GPa for ZrSiO4 and HfSiO4, respectively, are associated with [010](100) slip system. This suggests that the resistance to the [010](100) shear deformation for MSiO4 is quite low.
The linear compressibility of a crystal is the relative decrease in length of a line when the crystal is subjected to unit hydrostatic pressure. In general, it varies with direction. For a tetragonal crystal, the directionality of linear compressibility is expressed as:
$$\beta = \left( {{s_{11}} + {s_{12}} + {s_{13}}} \right) - \left( {{s_{11}} + {s_{12}} - {s_{13}} - {s_{33}}} \right)l_3^2\quad .$$The compressibility of MSiO4 is rotationally symmetrical about the four fold axis (c axis), thus it is sufficient to describe the characteristic of the directionality of linear compressibility by projecting the surface contour on (100) and (001) planes. The compressibility of MSiO4 is shown in Fig. 4. The linear compressibility of a (1.64 × 10−3 GPa−1 for ZrSiO4, 1.53 × 10−3 GPa−1 for HfSiO4) is much larger that of c (1.00 × 10−3 GPa−1 for ZrSiO4, 0.98 × 10−3 GPa−1 for HfSiO4) for both silicates. The linear compressibility of ZrSiO4 measured by Hazen et al. was 1.6 × 10−3 GPa−1 for a, 0.95 × 10−3 GPa−1 for c.Reference Segall, Shah, Pickard and Payne29 Our prediction reproduces the experimental values quite well.
Planar projections of linear compressibility on (100) and (001) crystallographic planes for (a) ZrSiO4 and (b) HfSiO4, respectively.
D. Thermal expansion coefficient
Thermal expansion coefficient (TEC) is a fundamental parameter of a material. In this work, the TECs of MSiO4 are evaluated by a simply but efficient method involved CBT. The detail procedure can be found in Refs. Reference Sanditov, Tsydypov and Sanditov43–Reference Carter and Norton45. All the chemical parameters of subformulas are calculated and tabulated in Table V αμ is the contribution of one particular bond to the average linear TEC of MSiO4. With strong chemical bonding and rigid characteristic, the response of Si–O bonds to temperature is much smaller than that of M–O bonds. The linear expansion coefficient of MSiO4 is 5.1 × 10−6 and 4.4 × 10−6 K−1 for ZrSiO4 and HfSiO4, respectively. The experimental TEC of ZrSiO4 is quite scattered, extends from 4.1 × 10−6 to 5.5 × 10−6 K−1.Reference Subbarao, Agrawal, McKinstry, Sallese and Roy57 Our prediction locates in this range, which means CBT is an effective method to estimate the TEC of MSiO4.
E. Thermal conductivity
The calculated sound velocities (υl, υs, υm), Debye temperature (ΘD) and Grüneisen parameter (γ) that are required for estimation of temperature-dependent thermal conductivity of MSiO4 are listed in Table VI. Since lower elastic moduli and higher density result in lower sound velocities, the sound velocities of HfSiO4 are much lower than those of ZrSiO4. Relatively low Debye temperature of HfSiO4 implies relatively weak atomic vibrations in the lattice. The Grüneisen constants γ are 1.72 and 1.83 for ZrSiO4 and HfSiO4, respectively. The higher Grüneisen constant of HfSiO4 indicates severer anharmonic phonon scattering in HfSiO4, compare to ZrSiO4.
Theoretical density d (in g cm−3), sound velocities (υl, υs, υm in km s−1), Debye temperature ΘD (in K), Grüneisen constant γ, and minimum thermal conductivity (in W m−1 K−1) of Polycrystalline MSiO4 ceramics.
The temperature-dependent thermal conductivity of MSiO4 is estimated from Slack's model, as shown in Fig. 5. In the temperature range where Slack's model is suitable, the thermal conductivity of MSiO4 is smaller than 5.0 W m−1 K−1. With increasing temperature, the thermal conductivity of MSiO4 decreases over the studied temperature range. As temperature increases further, the thermal conductivity would approach a minimum when phonons are completely uncoupled and energy is transferred between neighboring atoms, i.e., the phonon mean-free path decreases to the average atomic distance.Reference Slack42 The minimum thermal conductivity κmin can be evaluated by modified Clarke's model, as illustrated by Eqs. (3) and (4). The minimum thermal conductivity was predicted to be 1.54 and 1.24 W m−1 K−1 for ZrSiO4 and HfSiO4, respectively. Combining these two models, the behavior of thermal transportation of MSiO4 can be described as follow: with increasing temperature, the thermal conductivities decline as κ = 3133.89/T and κ = 2419.13/T for ZrSiO4 and HfSiO4, respectively. At around 1900 K, they reach the minimum value and keep unchanged with temperature.
Temperature dependence of thermal conductivity of MSiO4. Minimum thermal conductivities (dash line) are also shown.
At high temperatures, the minimum thermal conductivities are 1.54 and 1.24 W m−1 K−1 for ZrSiO4 and HfSiO4, respectively. Low thermal conductivity of MSiO4 indicates intense phonon scattering process, which could be traced back to the characteristics of crystal structure. The origin of low thermal conductivity of MSiO4 is mainly attributed to the heterogeneous bonding. The rigid Si–O tetrahedra are efficient in transporting heat, while softer M–O dodecahedra are equivalent to the thermal rattle structures, and provide “weak zones” that scatter phonons. In addition, as suggested by Clarke,Reference Clarke36 the strongest influence on the minimum thermal conductivity is the ratio of the atomic weight to the number of atoms per molecule, i.e., mean atomic weight M/n in Eq. (3). The mean atomic weight of HfSiO4 (45.1 amu/atom) is approximately 47% larger than that of ZrSiO4 (30.6 amu/atom), which results in lower thermal conductivity of HfSiO4.
IV. CONCLUSIONS
The structural, elastic, and thermal properties of MSiO4 were investigated by using the combination of DFT calculations and CBT. The calculated ground state properties including lattice constants, internal structure parameters and elastic constants are in good agreement with experimental values and previously reported theoretical values. Mulliken population analyses show that SiO4 tetrahedra are strongly constructed while MO8 dodecahedra are relatively weak and soft. This is confirmed by CBT. The polycrystalline mechanical properties, such as bulk moduli, shear moduli, Young's moduli, and hardness are estimated from theoretical elastic constants. Relatively large B/G and small H V3/E*2 suggest that MSiO4 are “quasi-ductile” ceramics. The discussion on the anisotropy in Young's modulus shows that both silicates show strong anisotropy in Young's modulus, and the smallest Young's modulus is in [110] directions for both silicates, which means [110] directions have the lowest fracture energy. Shear anisotropic factors show that the largest anisotropy in shear modulus is on (100) shear planes, and the lowest shear moduli, which are only 49.4 and 42.8 GPa for ZrSiO4 and HfSiO4, respectively, are associated with [010](100) slip system. The TECs estimated from CBT for these two compounds are 5.1 × 10−6 and 4.4 × 10−6 K−1, respectively. The intrinsic thermal conductivity of ZrSiO4 and HfSiO4 decreases as κ = 3133.89/T and κ = 2419.13/T, respectively. At high temperatures, it approaches a minimum value (1.54 and 1.24 W m−1 K−1 for ZrSiO4 and HfSiO4, respectively).
ACKNOWLEDGMENTS
This work was supported by the National Outstanding Young Scientist Foundation for Y.C. Zhou under Grant No. 59925208, and the National Natural Science Foundation of China under Grant No. U1435206.
