Relationships between elastic anisotropy and thermal expansion in A2Mo3O12 materials

We report calculated elastic tensors, axial Grüneisen parameters, and thermal stress distributions in Al2Mo3O12, ZrMgMo3O12, Sc2Mo3O12, and Y2Mo3O12, a series of isomorphic materials for which the coefficients of thermal expansion range from low-positive to negative. Thermal stress in polycrystalline materials arises from interactions between thermal expansion and mechanical properties, and both can be highly anisotropic. Thermal expansion anisotropy was found to be correlated with elastic anisotropy: axes with negative thermal expansion were less compliant. Calculations of axial Grüneisen parameters revealed that the thermal expansion anisotropy in these materials is in part due to the Poisson effect. Models of thermal stress due to thermal expansion anisotropy in polycrystals following cooling showed thermal stresses of sufficient magnitude to cause microcracking in all cases. The thermal expansion anisotropy was found to couple to elastic anisotropy, decreasing the bulk coefficient of thermal expansion and leading to lognormal extremes of the thermal stress distributions.


Introduction
Materials with the general formula A 2 M 3 O 12 , where A is a trivalent cation and M is W or Mo, have the unusual property of linear coefficients of thermal expansion spanning the range from low positive (2.4 Â 10 À6 K À1 ) to negative (À9 Â 10 À6 K À1 ), including zero. [1][2][3][4][5] Their structures have considerable chemical flexibility, as solid solutions can be made incorporating cations with large differences in size (e.g., Al 3+ and Y 3+ ). 3 Aliovalent substitutions of Zr 4+ or Hf 4+ with Mg 2+ for A 3+ can be performed, offering additional possibilities for the tailoring of physical properties within this family. [5][6][7][8] The structures of these materials are composed of vertex-linked networks of AO 6 octahedra and MO 4 tetrahedra; 9-11 negative contributions to thermal expansion arise from transverse motions of bridging oxygen atoms accompanied by distortion of the AO 6 octahedra. 9,12 Negative thermal expansion (NTE) has been correlated with the compliance of the A-O ionic bond, and the subsequent distortability of the AO 6 octahedra. 5,9,12 The open-framework structure allows for unusual thermal expansion and chemical flexibility, 4,12,13 and it typically leads to low stiffness in comparison to other oxide ceramics. 5,[14][15][16][17] The thermal expansion of materials from the A 2 M 3 O 12 family is known to be highly anisotropic. 15,18 Such anisotropy can cause large thermal stresses and subsequent microcracks in polycrystalline materials with randomly oriented grains. 19,20 However, the elastic anisotropy of these materials is generally unknown; furthermore, it is more difficult to determine experimentally than the anisotropy of the coefficient of thermal expansion (CTE). Critically, the mechanical properties determine the ability of thermomiotic (negative thermal expansion) 15 materials to counteract positive thermal expansion 15,21 and affect the stress distributions in polycrystalline A 2 M 3 O 12 materials. 15 Therefore, the elastic tensors of Al 2 Mo 3 O 12 , Sc 2 Mo 3 O 12 , and ZrMgMo 3 O 12 have been calculated and reported herein. Together with that of Y 2 Mo 3 O 12 (reported previously), 16 we now have elastic constants of an isomorphic family of materials ranging from low positive to large negative thermal expansion.
The study of the anisotropy of the elastic properties of these materials generally has previously been limited to determination of bulk moduli by variable-pressure XRD; the exception is Y 2 Mo 3 O 12 . 16 The variable-pressure XRD experiment can yield information about directional compressibilities; for example in Sc 2 W 3 O 12 17 and Sc 2 Mo 3 O 12 14 the compressibility along an axis was shown to be positively correlated with the thermal expansion along that axis. However, this technique has limitations; for example the compressibilities in orthorhombic Al 2 W 3 O 12 could not be accurately determined due to a pressure-induced phase transition, 14 and no information about shear elasticity was obtained. The use of ultrasonic measurement to obtain shear moduli of bulk polycrystalline samples 5,16,21 contributes some additional information, but experimental determination of the full elastic tensor would require single crystals for Brillouin spectroscopy or resonant ultrasound spectroscopy. Therefore, for practical reasons, computational methods are the best choice to study the full elastic anisotropy of A 2 M 3 O 12 materials, and this is the approach taken herein. There is an associated drawback: the DFT calculations are performed for T = 0 K. Since A 2 Mo 3 O 12 materials (except Y 2 Mo 3 O 12 ) transition to a monoclinic phase upon cooling, the calculations were performed on a phase which is not stable at absolute zero. Therefore the phononic structures contain vibrational modes with negative energy that destabilize the orthorhombic structure. [22][23][24] However, since the orthorhombic-monoclinic phase transition involves rotation of the coordination polyhedra, 2 the unstable mode is located away from the G point of the Brillouin zone 25,26 and the instability is not caused by negative terms in the elastic stiffness tensor. 27 Additionally, because many thermomiotic materials have highly temperature-dependent bulk moduli, 18,28-30 DFT results will not necessarily match room-temperature elastic constants, but our approach is validated by comparison with experimental behaviour at room temperature.
The energies of the optic phonons at the G point were also calculated using DFT. Calculation of the full dispersion relationship would be computationally expensive, however with the large unit cell of the orthorhombic A 2 M 3 O 12 materials (68 atoms) 1,5 there is little dispersion. 16,18 Low-energy optic phonon modes in A 2 M 3 O 12 materials corresponding to polyhedral distortions can have large negative Grüneisen parameters, 31,32 and therefore their energies play an important role in the thermal expansion behaviour. Determination of optic phonon frequencies also allows further validation of the calculations by comparison to experimental Raman spectra.

Methods
DFT calculations were carried out using the ABINIT 33 code in order to determine the elastic tensors and G-point phonon energy of A 2 M 3 O 12 and AMgM 3 O 12 materials. Computation of elastic tensors and G-point phonon energies were performed using the ACENET 34 and WestGrid 35 clusters, taking advantage of the massive parallelization functionality of ABINIT. 36 ABINIT composes wavefunctions as linear combinations of a plane wave basis set and performs calculations in reciprocal space. 37,38 Calculations of elastic tensors and phonon frequencies were performed using the response-function capability of ABINIT. [39][40][41] Phonon frequencies were corrected for errors due to finite sampling on the effective charges by imposition of the acoustic sum rule (i.e., requiring that the frequencies of the acoustic branches are zero at G). 39 The CIF2Cell program 42,43 was used to create input files from published crystal structure data. 5,44,45 PBE GGA exchange-correlation functionals were used in all cases. 46 In the cases of Al 2 Mo 3 O 12 and Sc 2 Mo 3 O 12 , norm-conserving 47 pseudopotentials generated using the OPIUM code 47-50 from the Bennet and Rappe pseudopotential library 51 were used. For ZrMgMo 3 O 12 , two-projector optimized norm-conserving Vanderbilt pseudopotentials 52 were used. These pseudopotentials were tested by comparison of calculated bulk moduli to experimental values for MgO (À7.2% deviation), 53 ZrO 2 (À1.6% deviation), 54 and Mo (À3.4% deviation), ‡ 55 yielding reasonable results. 56 The plane wave cut-off energy and k-point grid spacing were determined by convergence studies (Table 1); the criterion used was convergence of the internal pressure to within 1%. A k-point grid spacing smaller than 0.04 Å À1 was used in all cases. The structures were relaxed until the remaining pressure was less than 1 MPa. Changes in the unit cell axis lengths following relaxation are shown in Table 1 16 were used to create models of thermal stress in polycrystals of random texture, following a drop in temperature of 700 K, using the teFFT algorithm described in ref. 57. The teFFT algorithm simulates the elastic response to a change in temperature using a fast Fourier transform (FFT) approach to iteratively solve the constituent equations of thermoelasticity. 57,58 Since this method is only computationally limited by the speed of the FFT, simulations should scale with order n log n, where n is the number of nodes used in the model. Additionally, the spectral approach to solving the constituent equations of thermoelasticity only requires a structured rectilinear grid or image as an input, called the Fourier grid, forgoing the need for meshing. 57 These aspects resulted in considerable computational time savings by comparison to finite-element methods, allowing much larger models to be studied. The microstructure models used were generated in DREAM.3D 59 and consisted of 1190 randomly oriented, equiaxial grains on a 256 3 Fourier grid. A section of a modeled microstructure is shown in Fig. 1. Since the spectral approach requires periodic boundary conditions for the model, in the thermoelastic models polycrystalline structures were embedded in a compliant buffer layer. The buffer layer represents a region of infinite elastic compliance, which is created by  35 2 Â 3 Â 3 3.0 2.1 2.1 ‡ Mo was chosen rather than an oxide as the test material because MoO 3 has a layered structure and therefore is an unsuitable point of comparison, 56 and the bulk modulus of MoO 2 has not been reported in the literature. setting the stress values to zero, and is meant to approximate material surfaces. 55 The usage of buffer layers ensures that the effect of elastic anisotropy on the volume thermal expansion could be determined.
Relationships between the software packages and computational methods used in this work are shown graphically as a flowchart in the ESI. †

Elasticity
The calculated elastic stiffness tensor (c) of orthorhombic Al 2 Mo 3 O 12 (Pbcn setting) is given by: (1) and the corresponding directional elastic moduli are shown in Table 2. The thermal expansion tensor (a) from the literature 60 is shown for comparison of its anisotropy to that of c: (2) Al 2 Mo 3 O 12 showed elastic anisotropy similar to that of its thermal expansion, with only a small difference between the b-and c-axes. The overall elastic anisotropy of Al 2 Mo 3 O 12 , expressed as the universal anisotropy parameter, 61 A U , is equal to 1.9. By comparison. Y 2 Mo 3 O 12 has A U = 3.9. 16 Unlike in the case of Y 2 Mo 3 O 12 , 16 the axial Young's moduli and directional compressibilities followed the expected trends and were inversely proportional to each other. The Voight-Reuss-Hill (VRH) values 62 for the bulk elastic constants of Al 2 Mo 3 O 12 were determined to be K = 35.3 GPa and G = 36.2 GPa. These stiffnesses are ca. double those determined experimentally for monoclinic Al 2 Mo 3 O 12 (K = 13.5 AE 1.3 GPa, G = 18.2 AE 0.7 GPa). 63 A similar discrepancy was seen in the case of Y 2 Mo 3 O 12 , 16 obscuring whether the difference is due to temperature dependence of the elastic constants or the difference between the monoclinic and orthorhombic phases.
The calculated elastic stiffness tensor of orthorhombic ZrMgMo 3 O 12 (P2 1 nb setting) is given by: and the corresponding directional elastic moduli are shown in Table 3. The experimental thermal expansion tensor is shown for comparison of its anisotropy to that of c: 5 (4) As for Al 2 Mo 3 O 12 , ZrMgMo 3 O 12 shows increased stiffness along its thermomiotic axes both in terms of Young's modulus and compressibility. The stiffnesses of the b-and c-axes were similar even though their CTEs are quite different in magnitude. The overall elastic anisotropy of ZrMgMo 3 O 12 was very similar to that of Al 2 Mo 3 O 12 ; both have A U = 1.9. The VRH bulk and shear moduli were 26.5 GPa and 23.8 GPa, respectively. These values are close to those obtained experimentally (K = 31 AE 3 GPa, G = 22 AE 1 GPa). 5 The calculated elastic stiffness tensor of orthorhombic Sc 2 Mo 3 O 12 (Pbcn setting) is given by:   and the corresponding directional elastic moduli are shown in Table 4. The experimental thermal expansion tensor 2 is shown for comparison of its anisotropy to that of c: The elastic anisotropy of Sc 2 Mo 3 O 12 (A U = 2.9) was larger than that of Al 2 Mo 3 O 12 and ZrMgMo 3 O 12 , which could be expected given its larger CTE anisotropy. The VRH values for the bulk and shear moduli were K = 24.0 GPa and G = 22.3 GPa; unlike in the cases of Al 2 Mo 3 O 12 and Y 2 Mo 3 O 12 16 the calculated bulk modulus of Sc 2 Mo 3 O 12 is significantly less than the experimental value (K = 32 AE 2 GPa from variable-pressure XRD). 14 Correlations of axial elastic properties with axial thermal expansion, axial Young's moduli (Y ii ), axial compressibilities (b ii ), and directional shear moduli (G jk ) are visualized in Fig. 2 as functions of the axial CTEs (a ii ) 2,5,60 for the three materials described above and Y 2 Mo 3 O 12 . 3,16 The elastic behaviours of the thermomiotic b-and c-axes are correlated with their axial CTEs, while such a correlation is not seen for the a-axis. The b-and c-axes are stiffer than the a-axis except in the case of Y 2 Mo 3 O 12 , where the a-axis has a higher Young's modulus but also higher compressibility. Stiffer thermomiotic axes could decrease the bulk CTE of a polycrystal, as the increased stiffness would yield a larger contribution to the bulk expansion. As described below, this effect can be predicted in some cases. The trends in the VRH-averaged isotropic moduli (Fig. 2(d)) are more pronounced than those in the axial stiffnesses. A decrease in overall stiffness with increasingly negative thermal expansion can be seen clearly. When a thermomiotic material is used to counteract positive thermal expansion, its stiffness can be as important as its CTE in determining the amount of CTE reduction achieved, 15,20,21,63 and therefore the negative correlation between stiffness and CTE might act to discourage the use of materials such as Y 2 Mo 3 O 12 in applications. The trend relating the elastic moduli to the CTE is very similar to the previously reported correlation between the A site ionic force and the CTE, 5 suggesting that polyhedral distortability is a common factor between thermal expansion and elastic properties in this group of materials.
The materials studied show significantly different behaviour when subjected to different stress conditions. The b-and c-axes  show a decrease in uniaxial stiffness ( Fig. 2(a)) with increasingly negative thermal expansion. However, when subjected to isotropic pressure (Fig. 2(b)), less variation in the stiffness is seen. These differences illustrate the importance of computational studies of the full elastic tensor, as only b ii can be obtained from experimental variable-pressure XRD measurements on polycrystalline samples. The directional shear moduli generally show increased stiffness for G 23 by comparison with G 13 and G 12 , which is consistent with the increased stiffness of the b-and c-axes compared to the a-axis. This trend leads to the shear modulus in a plane (G jk ) increasing with the thermal expansion of the axis perpendicular to the plane (a ii  Fig. 3.

Phonon frequencies
The G-point optic phonon frequencies calculated in ABINIT are presented in tabular form in the ESI, † and visually in Fig. 4 (along with the previously reported 16 phonon spectrum of Y 2 Mo 3 O 12 ). Fig. 4 shows the common features of the phonon spectra of the four materials studied: a band of low-energy librational, translational, and bending modes separated from stretching modes with higher energy, with asymmetric stretches at lower energies than symmetric stretches. 65,66 As discussed in ref. 5, ZrMgMo 3 O 12 has a larger spread of stretching mode energies due to the splitting of the A site into a Zr site and Mg site, but otherwise its G-point phonon spectrum is similar to that of the A 2 Mo 3 O 12 materials. The energies of the low-energy phonons can be compared more easily in Fig. 4, which shows the cumulative distribution function of the modes below 500 cm À1 . The phonons below ca. 200 cm À1 are expected to have negative mode Grüneisen parameters and contribute to NTE, while those of higher energies are expected to have positive or near-zero Grüneisen parameters. 31,32 Fig. 4 shows that Y 2 Mo 3 O 12 has significantly more modes with energies below 200 cm À1 and has generally lower phonon energies in this region, which is consistent with its observed NTE over a wide temperature range. 4,9 Conversely, Al 2 Mo 3 O 12 has considerably fewer modes below 200 cm À1 , with a much broader distribution of its phonon frequencies in the region below 500 cm À1 . Of particular interest is a mode at 26 cm À1 in Al 2 Mo 3 O 12 which is anomalously low in energy by comparison to the other translational and librational modes in this material. At some point in the Brillouin zone, this mode could be the soft mode with negative frequency that causes the orthorhombic-monoclinic phase transition.

View Article Online
The calculated G-point phonon frequencies of Al 2 Mo 3 O 12 and ZrMgMo 3 O 12 are compared to experimental Raman spectra (measured at room temperature) in the ESI. † Overall, the results of this comparison, in addition to the comparison of experimental elastic moduli to calculated values presented above, show that the calculated elastic and phononic properties of these materials, while not quantitatively accurate, show reasonable agreement with the experimental data.

Axial Grüneisen parameters
The axial coefficients of thermal expansion are determined by the anharmonicity of the interatomic potentials, as expressed through Grüneisen parameters, and by the stiffness tensor, which represents the resistance of the lattice to thermal deformation. The present computational analysis does not allow the axial Grüneisen parameters (g ii ) to be determined directly, but using the calculated anisotropic elastic constants and phonon energies and the experimental CTE tensor the following expression can be invoked: 67 where V m is the molar volume, C V is the isochoric heat capacity, and c ii and a ii are stiffness and thermal expansion tensor elements, respectively. Here, C V was calculated using the Einstein and Debye models following the procedure described in ref. 68 using the calculated elastic tensors and G-point phonon frequencies. 16 The axial Grüneisen parameters were determined at T = 300 K for ZrMgMo 3 O 12 , Sc 2 Mo 3 O 12 , and Y 2 Mo 3 O 12 , 9,16 and at T = 550 K for Al 2 Mo 3 O 12 (temperatures chosen so that the materials are in the orthorhombic phase). The results are shown in Fig. 5 and Table 5. Thermal expansion along the b-and c-axes strongly correlates with the Grüneisen parameter along that axis, suggesting that vibrational modes with negative Grüneisen parameters along those axes are responsible for NTE. Essentially, this high degree of correlation is because the c ij a jj and c ik a kk terms in eqn (7) have opposite signs and nearly cancel. These terms are due to the Poisson effect, where a strain (in this case, thermally induced) along one axis causes a corresponding strain along the axes perpendicular to it. 69 As shown in Tables 2-4, some of the axial Poisson ratios in these materials are quite large, close to the theoretical limit of 0.5 for stable isotropic materials. 70 However, NTE along the b-and c-axes increases the CTE along the a-axis due to the Poisson effect. The positive thermal expansion along the a-axis in Y 2 Mo 3 O 12 is caused entirely by this effect, as g aa is negative. Similarly, the increase in thermal expansion along the a-axis in Sc 2 Mo 3 O 12 relative to ZrMgMo 3 O 12 and Al 2 Mo 3 O 12 is due to its larger NTE along the b-and c-axes, as all three materials have similar values of g aa . Therefore, the Poisson effect is responsible for some of the thermal expansion anisotropy which causes microcracking in sintered bodies of A 2 M 3 O 12 materials. 19,20 Thermal stress Models of thermal stress due to thermal expansion anisotropy showed significant thermal stress in all four modeled materials. Thermal stress extrema reached hundreds of MPa, sufficient to cause significant microcracking, as seen experimentally in Al 2 Mo 3 O 12 and Al 2 W 3 O 12 . 19,20 The thermomiotic b-and c-axes are, on average, placed in tension while the a-axis is in compression due to thermal expansion mismatch. The net stress on the material is zero, as the boundaries are unconstrained. An example cross section of the thermal stress distribution in Y 2 Mo 3 O 12 can be seen in Fig. 6. The stress extrema are found at grain boundaries where several misaligned grains meet, and can be expected to cause intergranular microcracking.
The effect of elastic anisotropy on the thermal stress distributions can be quantified by comparison to Kreher's model for elastically isotropic polycrystals. 71,72 Anisotropic models have lower mean strain energy densities, hWi, as shown in Table 6, indicating that elastic anisotropy acts to decrease thermal stress. In the elastically isotropic case, the average strain in the material is required to be equal to the thermal strain by the (unconstrained) boundary conditions and the conditions of static equilibrium. 71,72 When the coupling of elastic anisotropy and thermal expansion anisotropy is included, the thermal expansion of the stiffer axes is slightly more expressed than that of the more compliant axes, which allows some relaxation of the thermal stress. This effect is seen in consolidated polycrystals as deviations from the intrinsic linear CTE (Da c ), shown in Table 6. While Da c is small in comparison to the bulk CTEs of Al 2 Mo 3 O 12 and Y 2 Mo 3 O 12 , for ZrMgMo 3 O 12 and Sc 2 Mo 3 O 12 the deviation is of similar magnitude to the intrinsic CTE. Therefore, elastic anisotropy could explain the relatively large differences observed in dilatometric results compared with the intrinsic CTE (e.g. from XRD) for these materials. 2,5 The stress distributions within the polycrystals were fit to statistical distributions as described in the ESI. † The distributions in the different materials are generally similar, and using two lognormal distributions for the stress along each axis provides excellent fits to the data (see Fig. S4-S11, ESI †). As the size of the A site cation increases from Al 3+ to Y 3+ , increasing thermal expansion anisotropy is accompanied by decreasing stiffness (vide supra). These two factors have opposite effects on the magnitude of the thermal stresses, resulting in similar stress distributions and extrema. This similarity indicates that it will prove difficult to use chemical control of the A site cation to reduce thermal stress due to thermal expansion anisotropy. When the thermal stresses are sufficiently large to cause microcracking, a further deviation from the intrinsic CTE can be found at elevated temperatures due to crack growth and healing    causing expansion and shrinkage, respectively, of the polycrystalline body. 19,20 In all cases, the two lognormal distributions in each fit had opposite signs for their prefactors, leading to positive and negative stress extrema each being fit by a single unbounded lognormal tail. The utility of the lognormal distribution in fitting the extreme stresses is likely connected to some underlying physical phenomenon. The lognormal distribution is known to arise as a result of multiplication of a series of random variables, whereas the normal distribution arises from summation of a series of random variables. 73 Therefore, the lognormal fit to the stress distribution implies some multiplicative character to their origin. In anisotropic polycrystalline materials, thermal expansion mismatches cause strains in the material, and consequently reaction forces and stresses. In elastically isotropic materials, the interactions of reaction forces throughout the material are isotropic. This isotropy results in the stresses at any point being related to the strains caused by thermal expansion anisotropy throughout the remainder of material additively, and the normal thermal stress distribution is a logical result. However, in a material where the elastic constants are anisotropic, the transfer of stresses through the material is not so simple. The elastic anisotropy leads to the grains reacting differently to the applied strains based on their orientations, with increased stress in stiffer directions and vice versa. Therefore, elastic anisotropy introduces a multiplicative factor into the thermal stress distributions, since the relationship between the stress at a point and the strain at another point is modified by the variations in the stiffness tensor.

Conclusions
Materials in the A 2 Mo 3 O 12 and AMgMo 3 O 12 families were shown to have significant elastic anisotropy, a finding consistent with their characteristically large CTE anisotropy. The thermomiotic axes of these materials were shown to be stiffer than the PTE axes, which could act to decrease their bulk CTEs. However, the stiffness along the thermomiotic axes decreases with increasingly negative thermal expansion. The stiffness along the PTE axes did not strongly correlate with their CTEs. However, Al 2 Mo 3 O 12 , ZrMgMo 3 O 12 , and Sc 2 Mo 3 O 12 were shown to have similar axial Grüneisen parameters along their PTE axes, with differences in thermal expansion being driven by elastic factors. Calculated G-point optic phonon frequencies showed differences between the materials in the low-energy region which correlate with their thermal expansion behaviour, and the optic phonon data compared to experimental Raman spectra validated the results.
The thermal stress distributions calculated herein show that thermal expansion anisotropy can cause large thermal stresses in A 2 M 3 O 12 materials, enough to cause significant microcracking upon heating and cooling. The distributions of the thermal stresses are affected by the properties of the CTE and elastic tensors. The results indicate that the stress distributions in the four members of the A 2 Mo 3 O 12 family studied are quite similar.
In both cases, maximum tensile and compressive thermal stresses on the order of 0.5 GPa can be expected, with positive and negative extrema being well-described by lognormal distributions. The elastic tensors were found to couple with the thermal expansion tensors to produce small deviations in the bulk CTE, which were large by comparison to the intrinsic CTE in the cases of ZrMgMo 3 O 12 and Sc 2 Mo 3 O 12 .