I. INTRODUCTION
Materials with low thermal conductivity have long been investigated for their broad applications, especially as thermal insulation materials. Researchers usually follow a selection guideline proposed by D.R. Clarke Reference Clarke1 to search for candidate materials. It says that materials with low thermal conductivity should satisfy four principal conditions: a large molecular weight, a complex crystal structure, nondirectional bonding, and a large number of different atoms per molecule. Accordingly, many studies have been focused on candidates, such as rare-earth silicates, Reference Sun, Liu, Wang, Wang, Zhou and Hu2 zirconates, Reference Liu, Wang, Zhou, Liao and Li3,Reference Feng, Xiao, Zhou and Pan4 and phosphates, Reference Zhou and Liu5,Reference Du, Wan, Qu and Pan6 since introducing rare-earth elements is an effective way to increase the atomic mass as well as bonding heterogeneity. However, in view of the sustainable development and cost minimization of thermal insulators, it is imminent to explore a new strategy in which rare-earth elements are no longer crucial to reduce thermal conductivity.
As an aluminosilicate free of a rare-earth element, magnesium cordierite (Mg2Al4Si5O18) was taken as the cutting-in point to enrich the scope of material selection. It has been reported to own very low thermal conductivity at room temperature (RT) [∼3 W/(m·K)] Reference Kingery, Bowen and Unhlmann7,Reference Cao, Vassen and Stoever8 compared with other candidates (details summarized in Table I Reference Kingery, Bowen and Unhlmann7–Reference Sun, Zhou, Wang and Li10 ). At the same time, the density of Mg2Al4Si5O18 is much lower (∼2.51 g/cm3) than most thermal insulation materials. Recently, some lightweight materials with low thermal conductivity like MPO4, Reference Zhou and Liu5 MP2O7, Reference Xiang, Feng and Zhou11,Reference Xiang, Feng and Zhou12 and MSiO4 Reference Xiang, Feng, Li and Zhou13 have attracted more attention, because of their broad applications as thermal insulators, flame resistance, and fireproofing materials in fields of transportation and construction industry. Their lattice thermal conductivities were comprehensively predicted using first-principles calculations.
Density (g/cm3) and thermal conductivity [W/(m·K)] at room and high temperature of candidates for low thermal conductivity applications.

Besides the low density and low thermal conductivity, Mg2Al4Si5O18 engages our attention for its specific complex crystal structure. Many studies have shown that one promising approach to reduce thermal conductivity is to introduce guest atoms into the crystalline framework, e.g., clathrate and skutterudite, since the loosely coupled guests are associated with high thermal resistivity. Reference Fulmer, Lebedev, Roddatis, Kaseman, Sen, Dolyniuk, Lee, Olenev and Kovnir14 As is shown in Fig. 1, either of the two polymorphs Reference Winkler, Dove and Leslie15 of Mg2Al4Si5O18 provides a specific structure to host Mg atoms at the interstitial sites of the framework of [AlO4] and [SiO4] tetrahedra. It is expected that this characteristic structure may play a pivotal role in reducing the thermal conductivity of Mg2Al4Si5O18.
Crystal structures of (a) α-Mg2Al4Si5O18 and (b) β-Mg2Al4Si5O18, as well as (c) and (d) their projections along the c axis (Mg in gray, Al in red, Si in blue and O in green).

In addition, Mg2Al4Si5O18 is well-known for its ultralow coefficient of thermal expansion, good thermal shock resistance, excellent stability up to high temperatures, good chemical stability, and dielectric properties. Reference Mirwald16–Reference Dong, Feng, Dong, Wang, Yang, Gao, Liu and Meng23 It is usually used as refractory materials, electronic packaging materials, bio-ceramics, infrared radiation materials, filtering material for molten metals, and honeycomb catalyst carrier in automobile exhaust system. All these great properties and performance of Mg2Al4Si5O18 support it as a strong candidate for precisely evaluating the intrinsic lattice thermal conductivity.
In this work, first-principle calculations are performed to study the lattice thermal conductivity of two cordierite polymorphs, α-Mg2Al4Si5O18 and β-Mg2Al4Si5O18. The equilibrium crystal structure, Mulliken population and bonding characteristics are calculated as well. The full set of elastic constants or stiffness matrix, polycrystalline shear modulus, bulk modulus, Young's modulus, and Poisson's ratio are computed, from which sound velocities, Debye temperature and Grüneisen parameter are also derived. Anisotropic elastic stiffness, such as directional Young's modulus is discussed. Finally, temperature dependence of thermal conductivity as well as the minimum value are predicted by the Slack model Reference Slack24 and Clarke model, Reference Clarke1 respectively. Our results demonstrate that the low thermal conductivity of Mg2Al4Si5O18 is closely related to its complex crystal structure, and that crystal structure engineering provides us an alternative opportunity to discover lightweight thermal insulators.
II. CALCULATION METHODS
A. Computational details
The theoretical calculations were performed using first principles based on density functional theory Reference Ceperley and Alder25 as implemented in the CASTEP code. Reference Segall, Philip, Probert, Pickard, Hasnip, Clark and Payne26 Vanderbilt-type ultrasoft pseudo potentials (US-PP) Reference Vanderbilt27 and a plane-wave expansion of the wave functions were used. The local density approximation Reference Ceperley and Alder25 was selected to describe the exchange correlation functional. The energy cutoff was fixed to be 450 eV. The k-point sampling of the Brillouin zone was constructed by using the Monkhorst–Pack method Reference Monkhorst and Pack28 with 2 × 2 × 2 grids in primitive cells. All atomic positions and lattice parameters were fully relaxed using Broyden–Fletcher–Goldfarb–Shanno Reference Pfrommer, Côté, Louie and Cohen29 minimization scheme until the difference on total energy is within 5.0 × 10−6 eV/atom, maximum ionic Hellmann–Feynman force within 0.01 eV/Å, maximum stress within 0.02 GPa and maximum ionic displacement within 5 × 10−4 Å.
As color coded in Fig. 1, the T1 tetrahedral site is disorderly occupied by Al/Si in the ratio of 2 to 1 while the T 2 site have the ratio of 1 to 2. Thus, virtual crystal approximation (VCA) Reference Bellaiche and Vanderbilt30 within the US-PP method was used to model the fractional occupancy of Al/Si atoms of α-Mg2Al4Si5O18. VCA assumes that on each potentially disordered site there is a virtual atom which interpolates between the characters of the actual components. This method is reasonably accurate to model Al/Si disorder in silicates, such as hollandite and gehlenite, Reference Winkler, Pickard and Milman31 and it could well reproduce the crystal structure, structural stability, and bulk modulus.
After applying two strain patterns and calculating the resultant stress tensors, the second-order elastic coefficients of α- and β-Mg2Al4Si5O18 were derived by linear fitting of the stress–strain relationship, as implemented by Milman and Warren. Reference Milman and Warren32 Three positive and three negative amplitudes were applied for each strain component with a maximum strain value of 0.3% to avoid the nonlinear stress–strain relationship at higher strain values. Convergence tolerance for optimization on atomic internal freedoms was selected as the difference in the total energy within 1 × 10−6 eV/atom, the ionic Hellmann–Feynman force within 0.002 eV/Å and the maximum ionic displacement within 1 × 10−4 Å. The choice of present computational parameters ensures good agreement with experimental values of the crystal structure and elastic properties.
B. Structural distortion
For an equilibrium crystal structure, definition of the degree of distortion Δd of the coordination polyhedron Reference Makovicky and Balic-Zunic33 yields:
where d(M–O)n and <d(M–O)n> are the M–O bond length for the nth O atom and the average length of M–O bonds in the polyhedron, respectively. It is expected that the higher the value of Δd is, the larger the degree of distortion of the polyhedron is.
C. Elastic properties
Bulk modulus B and shear modulus G for crystalline aggregates can be estimated from single-crystal elastic constants with the well-established Voigt–Reuss–Hill approximation. Reference Hill34 According to the conversion formulas between elastic moduli, Reference Nye35 the Young's modulus E, and Poisson's ratio ν can be calculated according to the bulk modulus B and shear modulus G:
D. Lattice thermal conductivity
Slack suggested the following description for lattice thermal conductivity κ [W/(m·K)] in a temperature range in which heat conduction is dominated by Umklapp scattering among acoustic phonons: Reference Slack24,Reference Morelli, Jovovic and Heremans36
$$\kappa = {{{B_\gamma }\bar M\delta \theta _{\rm{D}}^3} \over {{\gamma ^2}N_{{\rm{cell}}}^{2/3}T}}\quad ,$$
where
$\bar M$
is the average mass of atoms in g/mol, δ3 is the average volume of atoms in m3, θD is the Debye temperature in K, γ is the Grüneisen parameter, N
cell is the number of atoms per primitive unit cell, and B
γ is a coefficient related to γ:
Due to the significant complexity of the crystal structure and chemical bonding in Mg2Al4Si5O18, it is difficult to access the mode Grüneisen parameter Reference Zwanziger37 from the phonon dispersion. Therefore, in this work, the acoustic γ Reference Zwanziger37–Reference Liu, Wang, Li and Zhou40 is calculated to represent the anharmonicity of phonon scattering. The method is based on an average of lattice vibrational frequency, i.e., equalizing vibration energy of acoustic phonon to the energy of elastic wave propagation in a continuum approximation of polycrystalline. Hence, the acoustic Grüneisen parameter γ can be directly derived from Poisson's ratio ν or sound velocities as Reference Sanditov, Mantatov, Darmaev and Sanditov38,Reference Leont'ev41 :
Further, we use the definition of Debye temperature Reference Slack24,Reference Anderson39 from elastic constants as follows:
where h is the Planck constant in J·s, k B is the Boltzmann constant in J/K, N A is the Avogadro constant, ρ is density in kg/m3, M cell is the mass per primitive cell in kg, and v m is the mean sound velocity determined by density and elastic constants Reference Liu, Wang, Li and Zhou40 :
$$\eqalign{ & {v_{\rm{l}}} = {\left( {{G \mathord{\left/ {\vphantom {G \rho }} \right. \kern-\nulldelimiterspace} \rho }} \right)^{1/2}} \cr & {v_{\rm{t}}} = {\left[ {{{\left( {{{B + 4G} \mathord{\left/ {\vphantom {{B + 4G} 3}} \right. \kern-\nulldelimiterspace} 3}} \right)} \mathord{\left/ {\vphantom {{\left( {{{B + 4G} \mathord{\left/ {\vphantom {{B + 4G} 3}} \right. \kern-\nulldelimiterspace} 3}} \right)} \rho }} \right. \kern-\nulldelimiterspace} \rho }} \right]^{1/2}}\quad , \cr}$$
here v l and v s are the longitudinal and transverse sound velocity, respectively.
At a high temperature limit, lattice thermal conductivity approaches its minimum value and remains independent to temperature thereafter. Clarke Reference Clarke1 assumed that the minimum phonon mean free path converges to the average interatomic distance in the material at high temperature and approximated different atoms in a molecule with an equivalent atom having a mean atomic mass of M cell/n cell. Therefore, the theoretical minimum thermal conductivity can be calculated by:
III. RESULTS AND DISCUSSION
A. Crystal structure and chemical bonding
Cordierite mainly exists in two thermodynamically stable polymorphs, including an Al/Si disordered hexagonal phase α-Mg2Al4Si5O18 (P6/mcc, Z = 2) and an Al/Si ordered orthorhombic phase β-Mg2Al4Si5O18 (Cccm, Z = 4). Reference Winkler, Dove and Leslie15,Reference Cohen, Ross and Gibbs42,Reference Meagher and Gibbs43 Figure 1(a) is the crystal structure of α-Mg2Al4Si5O18 and Fig. 1(c) is its projection along the c axis. Al and Si atoms are illustrated in red and blue so that the Al/Si disordered sites are represented as two-colored spheres with a fixed red/blue ratio. Figs. 1(b) and 1(d) are the crystal structure and c-axis projection, respectively, of β-Mg2Al4Si5O18. It is clear that both crystal structures essentially consist of the framework of [AlO4] and [SiO4] tetrahedra with [MgO6] octahedral coordination. To be specific, the [AlO4] and [SiO4] tetrahedra can be divided into two types: corner-sharing T2 tetrahedra forming a distorted hexagonal ring which is contra rotated by ∼30° in the adjacent T2 ring layer; and T1 tetrahedra connecting neighboring T2 rings as well as sharing two opposite edges with octahedra that coordinate Mg atoms. The Mg atoms are located in the distorted octahedra that share edges with three T1 tetrahedra. The main difference between α- and β-Mg2Al4Si5O18 is whether there exists Al/Si ordering. Specifically speaking, in the high-temperature phase α-Mg2Al4Si5O18 that is stable above 1450 °C, the Al and Si atoms are distributed over two sets of tetrahedral sites: three T1 and six T2 sites per formula unit. The Al and Si distribution involves: two Al and one Si disorderly occupying the T1 tetrahedra; and two Al and four Si disorderly present on the T2 tetrahedra. In the low-temperature phase β-Mg2Al4Si5O18 that is stable under 1450 °C, the tetrahedra at T1 sites are grouped into two non-equivalent types: T11 and T16, while the tetrahedra at T2 sites are denoted into three non-equivalent types: T21, T23, and T26. The third symbol of each label represents the symmetry operation (1 for identity map, 3 for the third-turn, and 6 for the sixth-turn) that correlates the atoms in equivalent sets. Within this new structure, Al and Si atoms are able to form a completely ordered distribution Reference Meagher and Gibbs43 (see details in Table II).
Tetrahedral sites in the Mg2Al4Si5O18 polymorphs.

The optimized lattice parameters of α- and β-Mg2Al4Si5O18 are compared with the experimental data Reference Winkler, Dove and Leslie15 in Table III. Based on the results, the deviations of theoretical lattice parameters from the experimental values for α-Mg2Al4Si5O18 are 1.08% for a and b, and 1.54% for c; and the deviations for β-Mg2Al4Si5O18 are 0.90% for a and b, and 1.10% for c.
Theoretical lattice parameters (a, b, c in Å and α, β, γ in °) of α- and β-Mg2Al4Si5O18, together with experimental data for comparison.

Table IV lists the degrees of distortion Δd of coordination polyhedra for α- and β-Mg2Al4Si5O18. It is observed that the bond lengths of T1 (T11 and T16) sites are slightly larger than those of T2 (T21, T23 and T26) sites. Besides, Si–O bonds are shorter than Al–O bonds, which are normally identified in aluminosilicates. Interestingly, the bond lengths at both T1 and T2 sites in α-Mg2Al4Si5O18 are slightly smaller than their counterparts (T11, T16 and T21, T23, and T26) in β-Mg2Al4Si5O18. Accordingly, the degree of distortions for all three kinds of polyhedra in α-Mg2Al4Si5O18 is smaller than those in β-Mg2Al4Si5O18. When modeling the Al/Si disorder in α-Mg2Al4Si5O18, the VCA method eliminates the difference between Si–O and Al–O bonds by the virtually averaged Si1/3Al2/3–O and Si2/3Al1/3–O bonds. Therefore, it neglects such effects as local distortions around atoms in α-Mg2Al4Si5O18. Furthermore, distortions of T1 (T11 and T16) and M polyhedra are almost zero while the degrees of distortion of tetrahedra at T2 (T21, T23 and T26) sites are higher than those at T1 (T11 and T16) and M sites. Combined with the research studies on structure dependence on temperature or pressure, Reference Milberg and Blair20,Reference Miletich, Gatta, Willi, Mirwald, Lotti, Merlini, Rotiroti and Loerting44,Reference Mirwald, Malinowski and Schulz45 the conclusion is easily drawn that T1 (T11 and T16) tetrahedra and M octahedra play more active roles than T2 (T21, T23 and T26) tetrahedra in responding to deformation. It is supposed that since T1 (T11 and T16) tetrahedra are less constrained by the symmetry compared to T2 (T21, T23 and T26) tetrahedra, they can accommodate distortion by rotating or bending about adjoining tetrahedrons. This mechanism results in near-zero distortion in T1 (T11 and T16) tetrahedra. On the contrary, T2 (T21, T23 and T26) tetrahedra restricted by the six-member rings are less free for response to strain. Therefore, the degrees of distortions of T2 (T21, T23 and T26) and M sites are much larger than those for T1 (T11 and T16) sites. And it is also noted that T1 (T11 and T16) and M sites are closely related to the interstitial Mg atoms. Therefore, the weakly bound Mg atoms, which can be regarded as the guest atoms within the framework of Al and Si tetrahedra, may reduce thermal conductivity by effectively scattering phonons. In addition, the distortion of Si tetrahedra is smaller than that of Al tetrahedra since Si–O bonds are stiffer than Al–O bonds. Mulliken population Reference Mulliken46 of Si–O and Al–O bonds in β-Mg2Al4Si5O18 is analyzed in Fig. 2 to support this explanation. The population of Si–O bond is larger than that of Al–O bond, indicating Si–O bond is stronger than the latter one. The Mulliken population of ionic M–O bond is negative and is not shown here.
Bond lengths (Å) and degree of distortion of polyhedra of α- and β-Mg2Al4Si5O18.

Mulliken Populations of Si–O and Al–O bonds in β-Mg2Al4Si5O18.

B. Elastic properties
The calculated elastic constants for the polymorphs are listed in Table V and compared with the experimental values and other theoretical data. Reference Bubeck47,Reference Toohill, Siegesmund and Bass48,Reference Vinograd, Perchuk, Gerya, Putnis, Winkler and Gale49 The number of independent elastic constants is 5 for a hexagonal system or 8 for an orthorhombic system. All calculated results of α- and β-Mg2Al4Si5O18 are used to validate the Born–Huang criteria Reference Born and Huang50 for mechanical stability: a) c 44 > 0, c 11 > |c 12|, (c 11 + c 12) c 33 > 2c 13 2 for the hexagonal system; and b) c 11 > 0, c 22 > 0, c 33 > 0, c 44 > 0, c 55 > 0, c 66 > 0, c 11 + c22 + c 33 + 2(c 12 + c13 + c 23) > 0, c 11 + c 22 − 2c 12 > 0, c 11 + c 33 − 2c 13 > 0, c 22 + c 33 − 2c 23 > 0 for the orthorhombic system. The results ensure the structural stabilities of the polymorphs under elastic strain perturbations.
Theoretical elastic constants (GPa) of α- and β-Mg2Al4Si5O18 together with experimental data and another calculation for comparison.

Among all elastic constants, c 11, c 22, and c 33 represent the stiffness against principal strains, whereas c 44, c 55, and c 66 correspond to the resistance to shear deformations. A common feature for both phases is that c 11, c 22, and c 33 are much higher than c 44, c 55, and c 66. In addition, the ratios c 23/c 44, c 13/c 55, and c 12/c 66 deviating greatly from unity mean the existence of noncentral interatomic forces in the compounds. Moreover, the six main diagonal elements c ii of α-Mg2Al4Si5O18 are slightly larger than those of β-Mg2Al4Si5O18.
The elastic properties derived from elastic constants are summarized in Table VI, together with experimental values Reference Miletich, Gatta, Willi, Mirwald, Lotti, Merlini, Rotiroti and Loerting44,Reference Bubeck47,Reference Toohill, Siegesmund and Bass48,Reference Ogiwara, Noda, Shoji and Kimura51,Reference Koepke and Schulz52 for comparison. The polycrystalline bulk modulus B, shear modulus G and Young's modulus E of α-Mg2Al4Si5O18 are higher than those of β-Mg2Al4Si5O18. The low shear moduli G of cordierite probably originates from the character that tetrahedra at T1 sites are relatively softer than those at T2 sites (restrained by the six-member ring in the xy plane along z axis). To evaluate the trend of ductility or brittleness for cordierite, the ratio of B/G is calculated by Reference Pugh53 :
Theoretical polycrystalline bulk modulus B (GPa), shear modulus G (Gpa), Pugh's ratio G/B, Young's modulus E (GPa) and Poisson's ratio ν.

In materials with the B/G ratio higher than 1.75, local strain in the region of the fracture extends easily and thus it increases the ductility. On the contrary, decrement of the B/G ratio indicates a higher stress concentration at the root of a crack since plastic flow is relatively difficult. In the present case, the B/G ratios for α- and β-Mg2Al4Si5O18 are 2.44 and 2.49, respectively. It means that both phases are possibly quasi-ductile and that Al/Si-disordered α phase is a little more brittle than the ordered β phase.
The α- and β-Mg2Al4Si5O18 with hexagonal and orthorhombic crystal systems are expected anisotropic elasticity. We herein use the classical Zener Reference Zener54 index to measure the degree of elastic anisotropy. Three shear anisotropic parameters A 1, A 2, A 3 or namely Zener ratios are calculated for the (100), (010), and (001) shear planes:
$$\eqalign{ & {A_1} = {{4{c_{44}}} \mathord{\left/ {\vphantom {{4{c_{44}}} {\left( {{c_{22}} + {c_{33}} - 2{c_{23}}} \right)}}} \right. \kern-\nulldelimiterspace} {\left( {{c_{22}} + {c_{33}} - 2{c_{23}}} \right)}} \cr & {A_2} = {{4{c_{55}}} \mathord{\left/ {\vphantom {{4{c_{55}}} {\left( {{c_{11}} + {c_{33}} - 2{c_{13}}} \right)}}} \right. \kern-\nulldelimiterspace} {\left( {{c_{11}} + {c_{33}} - 2{c_{13}}} \right)}} \cr & {A_3} = {{4{c_{66}}} \mathord{\left/ {\vphantom {{4{c_{66}}} {\left( {{c_{11}} + {c_{22}} - 2{c_{12}}} \right)}}} \right. \kern-\nulldelimiterspace} {\left( {{c_{11}} + {c_{22}} - 2{c_{12}}} \right)}}\quad . \cr}$$
For a material with elastic anisotropy, all Zener ratios would obviously deviate from unity. As shown in Table VII, both phases have nonunity Zener ratios and the results demonstrate anisotropic elastic stiffness.
The shear anisotropy factors or Zener ratios of α- and β-Mg2Al4Si5O18 obtained from theoretical calculation.

Consequently, directional Young's moduli Reference Nye35 of α- and β-Mg2Al4Si5O18 are calculated by the equation expressed in the Einstein summation convention:
Here, n i , n j , n k , and n l are the direction cosines of angles with the three principal directions and S ijkl is the elastic compliance (inverse of stiffness tensor).
Figure 3 and Fig. 4 show the directional Young's moduli of both phases by the color map as well as three projections along the x, y, z axes. Directional dependences of Young's Modulus for α- and β-Mg2Al4Si5O18 display weak anisotropy. Hexagonal α-Mg2Al4Si5O18 are elastically isotropic in the xy plane, which means that the elastic constants are invariant to arbitrary rotations around the z axis [see Fig. 3(d)]. This condition is called transverse isotropy. Reference Ledbetter55 The root of this consequence is the symmetry restriction (c 66 = (c 11 − c 12)/2) in hexagonal crystals. In the meantime, orthorhombic β-Mg2Al4Si5O18 without the restriction is anisotropic in the xy plane [see Fig. 4(d)]. Young's modulus of β-Mg2Al4Si5O18 is 154.5 GPa along the x axis and 145.7 along the y axis. The maximum value of the Young's modulus is 165.6 GPa along all directions in the xy plane for α-Mg2Al4Si5O18 and 154.5 GPa along the x axis for β-Mg2Al4Si5O18. The minimum values of Young's modulus for α- and β-Mg2Al4Si5O18 are 137.2 and 128.6 GPa, respectively, along the z axis.
(a) Direction-dependent Young's modulus of α-Mg2Al4Si5O18 and its projection along (b) x, (c) y, and (d) z axes.
(a) Direction-dependent Young's modulus of β-Mg2Al4Si5O18 and its projection along (b) x, (c) y and (d) z axes.

C. Thermal conductivity
Lattice thermal conductivity is determined by various phonon scattering mechanisms, from boundary scattering to phonon–phonon Umklapp scattering with increasing temperature. Reference Tritt56,Reference Toberer, Baranowski and Dames57 Slack Reference Slack24 proposed the widely accepted model, i.e., Eq. (4), to depict the lattice thermal conductivity at a moderately high temperature. In such a case, the predominant scattering mechanism is the intrinsic phonon–phonon Umklapp scattering and the heat is transported primarily by acoustic phonons.
Table VIII presents the parameters used for the calculation of thermal conductivity from Slack and Clarke models [see Eqs. (4) and (10)]. As mentioned above, the sound velocities v l, v s, v m and Debye temperature θD are estimated from elastic constants [see Eqs. (7)–(9)]. Also, the Grüneisen parameter γ is obtained from Poisson's ratio ν [see Eq. (6)]. Figure 5 shows the temperature-dependent thermal conductivities predicted by the Slack model for both phases, together with some complex ceramics with low thermal conductivity. Reference Feng, Xiao, Zhou and Pan4,Reference Luo, Wang, Wang, Li and Hu58,Reference Tian, Sun, Wang and Wang59 The thermal conductivity at RT is very low, 2.95 and 2.64 W/(m·K) for α- and β-Mg2Al4Si5O18. As the temperature increases, thermal conductivity decreases with reciprocal T as described by Eq. (4). Thermal conductivity approaches its minimum value estimated by the Clarke model [see Eq. (10)], about 1.33 W/(m·K) for α-Mg2Al4Si5O18 at 605 K and 1.29 W/(m·K) for β-Mg2Al4Si5O18 at 558 K. As discussed above, the VCA method increases the symmetry and decreases the freedom of α-Mg2Al4Si5O18 compared to its real local structure. So the predicted thermal conductivity of α-Mg2Al4Si5O18 would be slightly overestimated. It is also expect that the actual value will to be very close to that of β-Mg2Al4Si5O18. Generally, the polymorphs have quite similar physical and chemical properties so they are not normally distinguished in the ceramic society. As the temperature increases further, thermal conductivity will no longer decrease because phonon mean free path reaches its lower bound limited by the average atomic distance δ. In other words, theoretical lattice thermal conductivity below the calculated minimum value is of no practical concern. From room to high temperature, the thermal conductivity of cordierite shows low magnitude and just changes within a narrow range, 2.9–1.3 W/(m·K), which is preferred in applications as thermal insulators.
Theoretical values of average atomic mass
$\bar M$
(g/mol), average distance between atoms δ (Å), number of atoms per primitive unit cell N
cell, density ρ (g/cm3), sound velocities (v
l, v
s, v
m in km/s), Debye temperature θD (K), Grüneisen parameter γ, and minimum thermal conductivity κmin [W/(m·K)] for polycrystalline cordierite polymorphs, together with Lu2SiO5 for comparison.

Temperature dependences of thermal conductivities of α- and β-Mg2Al4Si5O18, together with other complex ceramics for comparison. The theoretical minimum thermal conductivity is also shown as dashed line.

To better understand the physical origin of the very low thermal conductivity of Mg2Al4Si5O18, Lu2SiO5 Reference Tian, Sun, Wang and Wang59 reference material with the lowest thermal conductivity in Fig. 5 is selected as a comparison because: (i) Its crystal structure consists of similar polyhedral framework to Mg2Al4Si5O18 except that Mg2Al4Si5O18 contains loosely bound Mg atoms; (ii) Both of Lu2SiO5 and Mg2Al4Si5O18 are silicates, but the former contains heavy Lutetium atoms while the latter is free of the rare-earth element. These are representative materials for weighting the possible strategies to develop a novel candidate with very low thermal conductivity.
As a well-established method for evaluating thermal conductivity among numerous materials, the Slack model by itself provides important information for material selection. Equation (4) can be mathematically factorized into three parts, the number of atoms per cell N cell, Grüneisen parameter γ, and Debye temperature θD, by which thermal conductivity of a material is determined. The power values of N cell, γ and θD are −2/3, −2, and 3, respectively, in Eq. (4), so that a tiny increment of γ or N cell or decrement of θD would obviously decrease thermal conductivity. Mg2Al4Si5O18 has a higher Debye temperature θD (652 K) than Lu2SiO5 (423 K) but it still holds lower thermal conductivity (2.64 W/(m·K) for β-Mg2Al4Si5O18) than Lu2SiO5 [4.98 W/(m·K)] at RT. Therefore, θD is not the determinant factor of the low thermal conductivity of Mg2Al4Si5O18. In addition, they both have nearly the same atomic distance, δ = 2.35 Å, in the sparse crystal structures. Parameter δ is also screened out the key factor for the divergence of thermal conductivity.
The other two parameters, the Grüneisen parameter γ and number of atoms per cell N cell of Mg2Al4Si5O18, are the predominant reasons for its low lattice thermal conductivity. The Grüneisen parameter of Mg2Al4Si5O18 is about 30% higher than Lu2SiO5 (1.92 for Mg2Al4Si5O18 and 1.5 for Lu2SiO5), which may lead to more than 40% decrease in the magnitude of thermal conductivity. The larger γ indicates a higher level of anharmonic phonon scattering in Mg2Al4Si5O18, which is the origin of higher thermal resistance. The huge anharmonicity relates to the significant distortions of polyhedra and strong noncentral interatomic forces in Mg2Al4Si5O18, as well as the loosely bound Mg atoms within the AlO4 or SiO4 tetrahedra framework.
Furthermore, it is worthy of discussing the physical scheme on how the number of atoms in the primitive cell, N
cell, affects the lattice thermal conductivity of Mg2Al4Si5O18. It qualitatively characterizes the effects of structural complexity on lattice thermal conductivity. The N
cell is 58 for Mg2Al4Si5O18 (in contrast, only 32 atoms for Lu2SiO5). As the N
cell increases, the lattice parameters increase in real space and the position of the Brillouin zone boundary in reciprocal space moves inward and cuts the acoustic phonon branch with lower frequency. Simultaneously, the phonon dispersion “folds in” on itself as the N
cell increases.
Reference Luo, Wang, Li, Hu and Wang60,Reference Toberer, Zevalkink and Snyder61
These two physical mechanisms will decrease the group velocity v of acoustic phonons as well as the acoustic Debye temperature
Reference Toberer, Zevalkink and Snyder61,Reference Roufosse and Klemens62
(defined as
${\theta _{\rm{a}}} = {{{\theta _{\rm{D}}}} /\mathord{\left{\vphantom {{{\theta _{\rm{D}}}} {N_{{\rm{cell}}}^{{1 \mathord{\left{\vphantom {1 3}} \right\kern-\nulldelimiterspace} 3}}}}} \right\kern-\nulldelimiterspace} {N_{{\rm{cell}}}^{{1 /\mathord{\left{\vphantom {1 3}} \right\kern-\nulldelimiterspace} 3}}}}$
). Hence the lattice thermal conductivity is significantly decreased by increasing the number of atoms in the primitive cell.
The last discussion is addressed to the character of low density, or lightweight of Mg2Al4Si5O18. Mg2Al4Si5O18 has smaller average mass of atoms or lower density than Lu2SiO5, namely the density is only 1/3 of the value of Lu2SiO5. However, the lattice thermal conductivity of Mg2Al4Si5O18 yields a much lower value compared with its competitor Lu2SiO5 with a heavy rare-earth element. This result provides exciting information that it is possible to achieve very low thermal conductivity in lightweight materials. This is in contrast to our general knowledge on increasing the average atomic mass or incorporating rare-earth element to effectively block heat transportation. Reference Feng, Xiao, Zhou and Pan63–Reference Xu, Pan, Wang, Wan, Qi, Miao, Mori and Torigoe65 We herein highlight a perspective on developing a promising candidate with low thermal conductivity through tailoring the complexity of the crystal structure. To summarize, the lattice thermal conductivity is predominantly determined by Debye temperature, Grüneisen parameter, the number of atoms per primitive unit cell, as well as the mean atomic mass. To select an alternative strategy, giant anharmonicity of phonon scattering, giant structural complexity, and low density are believed to bring out the extremely low thermal conductivity in some lightweight materials like Mg2Al4Si5O18.
IV. CONCLUSION
Equilibrium structures, bonding characteristics, mechanical properties, and thermal conductivities of two cordierite polymorphs, α-Mg2Al4Si5O18 and β-Mg2Al4Si5O18, are predicted by first-principle total energy calculations. VCA method is verified by comparison between theoretical and experimental results. Crystal structure and chemical bonding analysis indicate that isolated Mg atoms are loosely bound with surrounding framework of Al/Si tetrahedra, generally resulting in the decrease of thermal conductivity. Elastic constants, polycrystalline moduli, Poisson's ratio, elastic anisotropy, and directional Young's modulus are discussed with respect to the crystal structure. Generally speaking, α-Mg2Al4Si5O18 has slightly higher mechanical properties than β-Mg2Al4Si5O18. Sound velocities, Debye temperature, and Grüneisen parameter are derived from the above parameters and are used to predict the temperature-dependent lattice thermal conductivity of cordierite and its minimum value by Slack's and Clarke's models. The intrinsic lattice thermal conductivities of α- and β-Mg2Al4Si5O18 decrease with temperature as 804.6/T and 719.7/T, respectively, and then approach the theoretical minimum values at high temperatures. The present results and related discussions suggest that the unique complex crystal structure of cordierite plays a key role in decreasing the thermal conductivity. Finally it is addressed that Mg2Al4Si5O18 is a promising lightweight ceramic with extremely low thermal conductivity.
ACKNOWLEDGMENTS
This work was supported by the Natural Sciences Foundation of China under Grant Nos. 51032006 and 51372252.











