Chemical bonding origin of the unexpected isotropic physical properties in thermoelectric Mg3Sb2 and related materials

The Mg3Sb2 structure is currently being intensely scrutinized due to its outstanding thermoelectric properties. Usually, it is described as a layered Zintl phase with a clear distinction between covalent [Mg2Sb2]2− layers and ionic Mg2+ layers. Based on the quantitative chemical bonding analysis, we unravel instead that Mg3Sb2 exhibits a nearly isotropic three-dimensional bonding network with the interlayer and intralayer bonds being mostly ionic and surprisingly similar, which results in the nearly isotropic structural and thermal properties. The isotropic three-dimensional bonding network is found to be broadly applicable to many Mg-containing compounds with the CaAl2Si2-type structure. Intriguingly, a parameter based on the electron density can be used as an indicator measuring the anisotropy of lattice thermal conductivity in Mg3Sb2-related structures. This work extends our understanding of structure and properties based on chemical bonding analysis, and it will guide the search for and design of materials with tailored anisotropic properties.

3 type compounds were rationalized by band structure engineering via an atomic orbital scheme 24 .
Despite there are intensive theoretical studies on how crystal orbitals affect the electronic structure and electrical transport, there is very little knowledge on the quantitative description of chemical bonding, especially the interlayer interaction and how the chemical bonding affects thermal properties in CaAl2Si2-type compounds.
Here we report a quantitative analysis of chemical bonding in an archetypical compound Mg3Sb2 based on Bader's quantum theory of atoms in molecules 25 , and compare it with the structurally related layered van der Waals solid SnS2. It is found that Mg3Sb2 possesses a nearly isotropic 3D chemical bonding network with the interlayer bond being mostly ionic with partial covalent nature, and comparable to the intralayer interactions. Such a unique bonding feature in Mg3Sb2 not only challenges the well-known Zintl formalism, but also results in nearly isotropic thermal expansion coefficients, lattice compression, atomic displacement parameters, and lattice thermal conductivity. Importantly, we show how a simplified parameter based on the electron density can be used as an indicator for the anisotropy of the lattice thermal conductivity.
Furthermore, the isotropic 3D chemical bonding network is found to be widely applicable to many other Mg-containing compounds with the CaAl2Si2-type structure.
AB2X2 with the CaAl2Si2-type structure can be described by tightly bound [B2X2] 2layers sandwiched by two-dimensional layers of A 2+ ions (Fig. 1a). Besides the interlayer A-X bond (d1), two types of bonds exist in the B2X2 slabs, i.e. the tilted and vertical B-X bond, where the vertical bond (d3) is often longer than the tilted bond (d2) 18 . Three nonequivalent atoms have completely different coordination environments: A is connected to six X atoms with six equal bonds, B is tetrahedrally coordinated by X atoms with three tilted B-X bonds and one vertical B-X bond, and X is coordinated by three A atoms and four B atoms with seven adjacent bonds 4 including three interlayer A-X bonds, three tilted B-X bonds, and one vertical B-X bond ( Fig.   1c and Supplementary Fig. 1).
Mg3Sb2 (Space group: P3 � m1, a = 4.56187 (3) and c = 7.22944(6) Å at 299 K) can be considered as a special case of the CaAl2Si2 (AB2X2) structure in which A and B are Mg1 and Mg2, respectively. To have a better understanding of the interlayer interaction, a layered metal dichalcogenide SnS2 with the trigonal CdI2-type structure (Space group: P3 � m1, a = 3.6456 (4) and c = 5.8934(11) Å at 300 K) 26 was chosen for comparison since it shares many structural similarities with Mg3Sb2. Without considering the difference in lattice parameters, Mg3Sb2 can be viewed as intercalating two monolayers of Mg ions into the van der Waals gap of SnS2 and replacing Sn and S respectively by Mg and Sb (Fig. 1b). Unlike the Sb atom surrounded by seven Mg atoms in Mg3Sb2, the S atom is coordinated by three Sn atoms and three S atoms with six adjacent bonds including three intralayer Sn-S bonds ( 1 d ′ ) and three interlayer S-S bonds ( 2 d ′ ) (see Fig. 1e).
The Non-Covalent Interaction (NCI) index, based on the electron density, ρ, and its derivatives, is a powerful tool to reveal weak interlayer interactions 27 . NCI analysis is based on the reduced density gradient (RDG) as a function of sign(λ2)ρ (see methods), where sign(λ2) is the sign of the second eigenvalue of the electron density Hessian matrix 27,28 . Negative values of sign(λ2)ρ indicate attractive interactions, whereas positive values suggest repulsive interaction. Spikes induced by the significant change in RDG approaching zero at critical points within low density regions correspond to weak interactions. The density value of the spike with low RDG relates to the strength of the corresponding interaction.
3D RDG isosurfaces with blue-green-red color (BGR) scales representing sign(λ2)ρ values are given in Fig. 1d,f. Dark green RDG isosurfaces indicate that the interlayer Mg1-Sb in Mg3Sb2 is an attractive interaction, stronger than the weak interlayer S-S interaction in SnS2 with RDG isosurfaces colored in green. In order to quantitatively understand the interlayer interactions, the dependence of RDG on sign(λ2)ρ is calculated and shown in Fig. 1g,h. As expected, distinct differences can be seen between the weak interlayer S-S and intralayer Sn-S interactions in SnS2. Compared with the intralayer Sn-S interaction, interlayer S-S interaction shows a low RDG peak with a much smaller sign(λ2)ρ value approaching zero, a clear indication of weak van der Waals interaction (Fig. 1h). In contrast, the RDG distribution of the interlayer Mg1-Sb interaction in Mg3Sb2 is very similar to those of the intralayer Mg2-Sb interactions ( Fig. 1g). The density value of the low RDG peak for the interlayer Mg1-Sb interaction is just slightly lower than those of the vertical and tilted Mg2-Sb interlayer interactions. This indicates that the interlayer and intralayer interactions in Mg3Sb2 are comparable; that is, the tilted Mg2-Sb bond is slightly stronger than the vertical Mg2-Sb bond, while the vertical Mg2-Sb bond is slightly stronger than the interlayer Mg1-Sb interaction.
Static deformation electron density maps of the (110) planes in Mg3Sb2 and SnS2 are shown in Fig. 1i Supplementary Fig. 3).
However, the nature of chemical bonds cannot simply be judged by the sign of ∇ 2 ρ(rb). A more accurate analysis is using the Laplacian profile in comparison with that of the Independent Atom Model (IAM). As can be seen in Supplementary Fig. 4, for all three types of bonds in Mg3Sb2 the Laplacian values along the bond path are less positive than those of IAM. This indicates the partial covalent nature (high polarity) in these bonds. Considering the degree of differences in Laplacian profiles, the proportion of partial covalency is slightly increasing from the interlayer Mg1-Sb to the vertical Mg2-Sb to the tilted Mg2-Sb.
According to the well-established classification scheme 29 , the interlayer Mg1-Sb and intralayer Mg2-Sb bonds can be described as polar bonds according to the BCPs properties with the small density ρ(rb), positive ∇ 2 ρ(rb), |V|/G being slightly larger than unity, negative total energy density H, and G/ρ < 1 (see Table 1). Changing from the tilted Mg2-Sb to the vertical Mg2-Sb to the Mg1-Sb bond, the density value undergoes a negligible decrease, which indicates a minor decrease in covalency and interaction strength. In contrast, a remarkable difference is 7 observed between the interlayer S-S and intralayer Sn-S bonds in SnS2 (see Table 1). The interlayer S-S interaction can be described as a weak van der Waals bond as judged by the very small electron density ρ(rb), positive ∇ 2 ρ(rb), |V|/G < 1, and positive H at the BCP, whereas the intralayer Sn-S interaction can be treated as a polar covalent bond based on the positive ∇ 2 ρ(rb), 1< |V|/G < 2, negative H, and G/ρ < 1. The electron density at the BCP of the interlayer S-S bond is approximately 8 times smaller than that of the Sn-S bond, indicating the much weaker strength of the interlayer interaction compared with that of the intralayer interaction. Upon comparison of topological properties at the BCPs (see Table 1), the three similar bonds in Mg3Sb2 are clearly more polar and weaker than the Sn-S bond in SnS2, but they are stronger than the weak interlayer S-S interaction. The above results are consistent with the NCI analysis.
The topological properties at BCPs of another archetypical CaAl2Si2-type compound, CaZn2Sb2, are analyzed and compared with those of Mg3Sb2 (see Table 1 The comprehensive chemical bonding analysis paves the way to understand the structure and properties. The largely ionic feature with partial covalency (high polarity) of chemical bonds in Mg3Sb2 explains the intrinsically poor carrier density and mobility as well as the reasonably low lattice thermal conductivity 30 . Furthermore, the comparable interlayer and intralayer interactions unveil the three-dimensional chemical bonding network in Mg3Sb2, ruling out a description as a typical 2D layered structure. Importantly, the 3D bonding network is nearly isotropic in Mg3Sb2 upon quantitative comparisons of topological properties of different bonds. Such a feature is decisive for many unique properties including structural parameters and thermal properties.  Table 1). As illustrated in the figure, the lattice parameters all display linear increasing trends as the temperature increases and the thermal expansion along the c axis is slightly larger than that along the a axis. The linear thermal expansion coefficients at room temperature along the a and c directions in Mg3Sb2 are respectively 1.88×10 -5 and 2.42×10 -5 K -1 , which leads to a nearly isotropic αc/αa of 1.29, much smaller than those of typical layered materials TiS2 31 , MoS2 32 , MoSe2 32 , and Bi2Te3 33 (see Fig. 2b and Supplementary Table 2).
Furthermore, similar results are observed in the pressure-induced lattice compression. The relative lattice parameters and interlayer distance as a function of a series of hydrostatic pressures are simulated using DFT calculations and given in Fig. 2c,d. Under the same pressure, the decrease of the relative lattice parameter c/c0 in Mg3Sb2 is just slightly larger than that of a/a0, whereas the lattice parameter c is much more compressible than a in all layered van der Waals solids SnS2 26 , TiS2, MoS2, and MoSe2.
It is clear that both the lattice expansion with temperature and the lattice compression under pressure exhibit nearly isotropic features in Mg3Sb2, which can be essentially understood by the 3D chemical bonding network in this material. In both cases, the lattice parameter c in Mg3Sb2 shows slightly larger thermal expansion or pressure compression than that of a, which can be attributed to the slightly weaker interlayer bond compared with the intralayer bonds.
Moreover, the interlayer distance of Mg3Sb2 is less compressible than those of van der Waals solids, confirming the interlayer interaction being stronger than the van der Waals force.
To gain insight on how chemical bonding affects the thermal motion of the atoms, the isotropic atomic displacement parameters Uiso of Mg3Sb2 were obtained from Rietveld refinement of multi-temperature synchrotron powder X-ray diffraction (PXRD) data (Fig. 3a).
Refinement details are provided in Supplementary Information (see Supplementary Table 1  When considering the atomic displacement parameters and potential wells along the axial directions, both the Sn and S atoms in layered SnS2 display highly anisotropic characteristics with the vibration along the c direction being significantly larger than along the a direction (see It is well known that weak chemical bonding usually accompanied by strong lattice anharmonicity leads to low lattice thermal conductivity 34,35 . To probe the effect of chemical bonding on thermal transport, the lattice thermal conductivity was simulated by DFT calculations. Indeed, weak van der Waals interaction between the layers in SnS2 leads to considerably lower lattice thermal conductivity along the c axis in comparison with that along the a axis (see Fig. 4a). High anisotropic ratio κa/κc of lattice thermal conductivity is commonly observed at 300 K in typical layered materials, such as 15.8 in SnS2, 16.2 in TiS2, 16.2 in MoS2 3 , and 11.2 in MoSe2 3 (Fig. 4b). In contrast, unlike the noticeable anisotropy in van der Waals solids, the lattice thermal conductivity in Mg3Sb2 is nearly isotropic with κ along the c axis being negligibly lower than that along the a axis (κa/κc≈1.1 at 300 K, see Fig. 4).
In order to elucidate the origin of the isotropic lattice thermal conductivity in Mg3Sb2,  Table   3).
The origin of anisotropy in thermal properties can be traced to the chemical bonding. We can define a simplified parameter, the intralayer-to-interlayer bond-strength ratio ρ intra /ρ inter , which measures the degree of anisotropy of the chemical bonding network in a layered structure.
ρ intra and ρ inter denote the electron density values at BCPs of the intralayer and interlayer bonds, respectively. For the layered AB2X2 compounds with two nonequivalent intralayer bonds, ρ intra is calculated by averaging the electron density values at BCPs of the two intralayer bonds. As can be seen in Fig. 4b, a nearly linear correlation between the anisotropic ratio κa/κc of lattice thermal conductivity and ρ intra /ρ inter is revealed. This suggests that ρ intra /ρ inter can be adopted as an indicator for the anisotropy of lattice thermal conductivity.
ρ intra /ρ inter ≈1 in Mg3Sb2 indicates a nearly isotropic 3D chemical bonding network, which results in the nearly isotropic feature in phonon dispersion, group velocity, Grüneisen parameter, and eventually in lattice thermal conductivity.
It should be noted that the nearly isotropic 3D bonding network is not limited to Mg3Sb2.  Supplementary Fig. 9). Synchrotron PXRD patterns were analyzed using Rietveld refinement in the FullProf program 37 . The where ( , ) i γ q is the mode Grüneisen parameter for the phonon branch i at wave vector q, given as where ( , ) i ω q is phonon frequency, V is volume, and    plane.

Table 1 | Topological properties of the bond critical points (r b ). d is the bond length. ρ(rb)
and ∇ 2 ρ(rb) are the charge density and its Laplacian at the BCP, respectively. G and V denote the kinetic and potential energy at the BCP, respectively. H is the total energy density (H = G + V). G, V, H and G/ρ are in a.u.     Supplementary Figure 11. Convergence of the lattice thermal conductivity of Mg3Sb2 with respect to the cutoff distance (rcutoff) corresponding to the interaction range from the third to the seventh nearest neighbors for the anharmonic calculations.

Supplementary Tables
Supplementary Table 1. Rietveld refinement details of the synchrotron PXRD data of Mg3Sb2 at 770-299 K upon cooling. The R factors and χ 2 shown here are the data from the Mg3Sb2 phase.
Tactual represents the actual temperature calibrated by the thermocouple.

Supplementary Note 1. Topological analysis of electron density
The quantitative analysis of chemical bonding in this work is based on Bader's quantum theory of atoms in molecules (QTAIM) 3 where n(r) is the a unit vector perpendicular to the surface S at r. An atom is described as the union of an attractor and its basin. The atomic property of an atom is then calculated by the integration within its atomic basin.

Supplementary Note 2. Synchrotron PXRD patterns and Rietveld refinement
From multi-temperature synchrotron PXRD patterns of Mg3Sb2 shown in Supplementary Fig.   9, we can see that the Sb phase appears as the temperature increases. To keep consistence, the secondary phase Sb is included in the refinements of all temperature points. Although there is small difference between the thermal expansion coefficients calculated by the data upon heating and cooling, the anisotropic ratio of thermal expansion coefficient upon heating and cooling is nearly the same (see Supplementary Tables 8 and 2). As shown in Supplementary Table 7, the amount of Sb phase is gradually increasing with increasing temperature from 299 K to 810 K, whereas the amount of Sb phase becomes very stable upon cooling from 770 K to 299 K (see Supplementary Table 1). To avoid the possible influence from the increasing Sb phase upon heating, we thereby use the cooling data at 770-299 K with the stable amount of Sb in the main text. The underlying mechanism of the appearing Sb phase upon heating in Mg3Sb2 powder will be discussed in our future work.
The Debye temperature was extracted by fitting the isotropic atomic displacement parameters Uiso based on a Debye model: 5-7 where ħ is the reduced Planck constant, T the absolute temperature, m is the mass of the atom, kB is the Boltzmann constant, ΘD is the Debye temperature, and d is a disorder parameter 8

Supplementary Note 3. Theoretical calculation
The Grüneisen parameters for the acoustic modes might be negative, which will lead to the cancellation between the acoustic and optical modes. In addition to the average method used in the main text, the average Grüneisen parameter is calculated using the sum over the squared  Supplementary Table 10.
The lattice thermal conductivity of TiS2 was calculated by ShengBTE code 10 based on a full iterative solution to the Boltzmann transport equation for phonons. The second-order and thirdorder interatomic force constants were computed in the 4×4×2 supercells. The displacement 51 amplitude of 0.09 Å was adopted for harmonic force constants calculations to ensure the wellconverged properties 11 . Van der Waals functional optB86b-vdW 12 in VASP 13 code was used for all calculations with an energy convergence criterion of 10 −6 eV and a plane wave energy cutoff of 600 eV.