Elasticity of Mg3Bi2-xSbx

Mg3Bi2-xSbx is a promising thermoelectric material working around room temperatures. Compared to electronic and thermoelectric properties, its mechanical properties are of great importance in practical applications but much less understood. Herein, we have systematically studied the elasticity of Mg3Bi2-xSbx by means of first-principles calculations with a large supercell of 40 atoms. We demonstrated that the 10-atom-unitcell is undersized with improper electronic structures. With the elastic constants, we have explored the comprehensive elastic features and the three-dimensional distribution of fundamental characteristics of Young’s modulus and Poisson’s ratio and their variation with respect to the Sb content x. We interpolate the variation in terms of the valence electron concentration. We have further examined the hardness, ductility, anisotropicity, and Debye temperatures. The elasticity exhibits strong anisotropy where the maxima are approximately three times larger than the minima for modules. A nearly linear dependence is also observed on the Sb content except x in the vicinity of 0.5. Our atomistic insights on elasticity might be helpful in the material design of thermoelectrics with desirable mechanical properties. Our work could serve as a map for tuning the mechanical properties of Mg3Bi2-xSbx and guide the possible synthesizing of novel thermoelectric material.


Introduction
Materials with the ability to convert electricity into thermal energy are required in the wide application of thermoelectrics [1]. Compared to the state-of-the-art Bi 2 Te 3 [2] or Ag 2 Se [3] alloys at around room temperature, Mg 3 Bi 2-x Sb x is much cheaper, free from the expensive elements, but has better thermoelectric performance [4][5][6] and thereby has the promising commercial prospects [7,8]. A figure of merit of zT ≈ 1.5 at 716 K for Mg 3 Bi 1.5 Sb 0.5 is reported [9]. It is comparable to but cheaper than BiTe alloys, whose, for instance, peak zT of 1.4 at 100 • C is captured [10]. So far, Mg-Bi-Sb materials have garnered increasing attention as their distinguished thermoelectric properties have been reported both by experiments and theories [11][12][13][14]. The thermal, electrical, and thermoelectric properties are extensively studied [15], such as the band topology [16][17][18][19], phonon dynamics [20][21][22][23][24], and topological thermoelectricity of nodal-line semimetal [25][26][27][28]. The intrinsic defects, such as grain boundaries, make the synthesizing of Mg 3 Bi 2-x Sb x challenging in the past few years. However, Tamaki et al. have observed the prominent of high zT, which is expected to appear with negligible grain boundaries for the high Bi/Sb ratio at three [9]. Mg 3 Bi 2-x Sb x has a higher carrier mobility and low thermal conductivity arising from the disordered structure [29]. Aside from the microstructures, the intrinsic properties of these site-disordered materials are mainly responsible for their outstanding thermoelectric behaviors.
As one of the fundamental properties, elasticity is required to completely describe the mechanical behavior of this kind of cubic material [30], detailing the elastic stability, elastic anisotropy, ductile manner, and so on, which can surely serve as a reference for the further experimental and theoretical studies and applications on these thermoelectric compounds. First-principles calculations have been widely used to explore the underling physics from the atom-scale perspective. However, the mechanical properties of Mg 3 Bi 2-x Sb x are much less studied compared to the electronic and thermoelectric properties. The variation of the mechanical properties with respect to the Sb content is still elusive.
In this work, we have systematically investigated the elasticity of Mg 3 Bi 2-x Sb x via first-principles calculations within the framework of density functional theory. We have compared the results from two models: one within 10-atom-unitcell and the other 40atom-unitcell. The mechanical response of the Bi/Sb ratio of 3 is scrutinized. Then we studied the influence of Sb potency on the elastic properties. They have a consistent linear trend except for several potencies around 0.5. The three-dimensional distribution of anisotropic Young's modulus and Poisson's ratio are illustrated for the Sb content x ranging from 0 to 2. The general trends of hardness, ductility, anisotropicity, and Debye temperatures are further examined.

Materials and Methods
Both Mg 3 Bi 2 and Mg 3 Sb 2 have a crystal structure in the trigonal P-3m1 space group with a group number of 164. The primitive unitcell of Mg 3 Bi 2 contains 3 Mg atoms and 2 Bi atoms. The crystalline structures of Mg 3 Bi 2-x Sb x (0 ≤ x ≤ 2) adopt the same space group of Mg 3 Bi 2 . The 2 Bi sites are replaced by Bi or Sb atoms according to their stoichiometry, as shown in Figure 1a. The symmetric repetitive within all enumerated 2 × 2 × 2 supercell structures of 40 atoms is eliminated by the open-source program Disorder [31]. These structures are then optimized by rough convergence criteria of first-principles calculations. The structure of minimum energy within all optimized structures at each Sb content x is the representative configuration. For the elastic property calculations, more precise structural optimizations are then used for configurations at each Sb content x. Twelve forced deformations, corresponding to positive and negative micro-strains for six independent components, are applied respectively to the configuration. Elastic constants can be obtained by the corresponding elastic response, manifested as six stresses components from the first-principles results.
Our first-principles calculations are within the framework of Kohn-Sham density functional theory [32] as implemented in the Vienna Ab initio Simulation Package (VASP) [33]. The interactions between atoms are described by Projector Augmented-Wave (PAW) [34] method. Generalized gradient approximation (GGA) developed by Perdew-Burke-Ernzerhof (PBE) [35] was adopted for the exchange-correlation functional. K-meshes are gamma-centered Monkhorst-Pack scheme [36] of 4 × 4 × 2 for structures of 40-atom unitcells and 9 × 9 × 2 for those of 10-atom unitcells. The cutoff energy for plane wave expansion is 260 eV. The order of Methfessel-Paxton method [37] is 1, which is responsible for the partial occupancies for each orbital, and its smearing width is 0.2 eV. The convergence accuracy is set to be 10 −5 and 10 −8 eV in the structure searching and configuration optimizing for the electronic self-consistent iterations. Energy convergence criteria of 10 −4 eV are used for the fast optimization in structure searching, while the force convergence criteria of 10 −3 eV/Å are set for that of accurate structural optimization. The three-dimensional distribution data is generated by VASPKIT package [38] such as Young's modulus, Poisson's ratio, and shear modulus, which is simply expressed in terms of the compliance matrix [39]. The compliance matrix is the inverse matrix of the rotated elastic constant whose coordinates is parallel to the appointed orientations. Materials 2022, 15, x FOR PEER REVIEW 3 of 12 Our first-principles calculations are within the framework of Kohn-Sham density functional theory [32] as implemented in the Vienna Ab initio Simulation Package (VASP) [33]. The interactions between atoms are described by Projector Augmented-Wave (PAW) [34] method. Generalized gradient approximation (GGA) developed by Perdew-Burke-Ernzerhof (PBE) [35] was adopted for the exchange-correlation functional. K-meshes are gamma-centered Monkhorst-Pack scheme [36] of 4 × 4 × 2 for structures of 40-atom unitcells and 9 × 9 × 2 for those of 10-atom unitcells. The cutoff energy for plane wave expansion is 260 eV. The order of Methfessel-Paxton method [37] is 1, which is responsible for the partial occupancies for each orbital, and its smearing width is 0.2 eV. The convergence accuracy is set to be 10 −5 and 10 −8 eV in the structure searching and configuration optimizing for the electronic self-consistent iterations. Energy convergence criteria of 10 −4 eV are used for the fast optimization in structure searching, while the force convergence

Results
Resembling the bi-components solid solution, the initial-generated structures of the larger supercell are duplicated from the crystalline units, shown in Figure 1a. Mg atoms can be divided into two kinds, denoted as Mg1 (yellow) and Mg2 (brown), according to distinctive polyhedrons of neighbor Bi sites. For Mg1, it is surrounded by six Bi sites, forming an octahedron environment while that of Mg2 is four, possessing a diamond-like tetrahedron environment. Among all enumerated structures of 40 atoms, the symmetric repetitive ones are eliminated. The rest structures are then optimized to search for an energyminimized structure. The final chosen structure of the lowest energy is used to represent real-world states because the most stable signifies the highest possibility of existence.
For instance, we take the structures of Mg 3 Bi 1.5 Sb 0.5 as an example since this compound possesses the best thermal-electronic effect. In general, electronic charge density is the most direct measurement of atomic interaction. The difference can be visualized intuitively from the electron location function (ELF) [40,41] of a horizontal plane crossing the most Sb atoms ( Figure 1b) and Sb atoms distribution for the configuration of 10-atom unitcell and 40-atom unitcell, respectively ( Figure 1c). The maximum of ELF only appears around the Sb ions owing to the smaller radius compared to Bi ions. Thereby, the stronger interactions around the Sb ions are typically deduced from the mechanistic perspective.
Due to the periodic images along the x 1 and x 2 axes, which is the intrinsic shortcoming of a small supercell, the Sb atoms take a stratiform distribution for a 10-atom unitcell. The structural stability of different grid arrangements with different ordered ions but the same Sb content, as well as the possible different lattice distortion patterns even with the same ordered ions, cannot be simulated. While in the 2 × 2 × 2 duplicated structure of a 40-atom unitcell, Sb ions accumulate with other Sb ions, giving a trend of precipitation that may trigger differences in the prediction of vast properties, especially the significant thermoelectric properties, to some extent.
We have investigated the elastic properties for Mg 3 Bi 1.5 Sb 0.5 using a 40-atom unitcell since the structure of the larger supercell is screened from the thermal stability. The 40-atom unitcells are more reasonable than the 10-atom unitcell owing to the larger configurational entropy. The elastic constants are firstly obtained from the DFT calculations. With these elastic constants, the distribution of anisotropic Young's modulus, shear modulus, and Poisson's ratio in the three-dimensional Cartesian coordinates are then computed and displayed in Figure 1d. The axes x 1 and x 3 are parallel to the lattice vector a and c defined in Figure 1a, respectively.
These elastic properties are extremely orientation-dependent. The maximum Young's modulus is 60 GPa along the direction of the x 3 axis, while the minima is 24 GPa within the cubic-shaped faces. For shear modulus and Poisson's ratio, another investigated direction is needed, which is perpendicular to the normal vector of the plane provided. Thus, a range of values for shear modulus and Poisson's ratio can be obtained in a specific direction. Merely the maximum and minimum values are what we are interested in, and they are plotted together in Figure 1d (outer transparent surface to be the maxima and inner surface to be the minima) to be easily acquired and clearly compared. For Poisson's ratio, not only the maximum and minimum values, from −0.01 to 0.70, but also the value ranges for different orientations vary greatly.
For shear modulus, it exhibits six local maximum values, the directions of which are vertical to each other, where the shear resistance is strong in all directions. The maximum value is 21.9 GPa, and the minimum value is 18.0 GPa in these directions. Three of the upper halves make the same angle with x 3 -axis. While on the x 1 -x 2 plane, it satisfies the trigonal symmetry that midpoints of six edges are in the x 1 -x 2 plane, corresponding to six apexes of a honeycomb structure. Conversely, the direction, which is in the middle of three adjacent maximum directions that are at 90 • to each other, has the weakest stiffness in all directions, with the maxima and minima being 10.7 GPa and 8.4 GPa. For directions not in the vicinity of the above 14 specific directions, they have large shear anisotropy, which represents a large maximum and small minimum shear resistance.
Based on the ability to search for the most stable structures at any Sb content x, mechanical properties strongly depend on the relative concentration of Sb and Bi. We scrutinized the elastic constant combinations at all available Sb contents in our model. The crystalline structure belongs to the trigonal system, possessing eight independent components: C 11 , C 12 , C 13 , C 33 , C 14 , C 15 , C 45 , and C 44 . Owing to the distortion of asymmetry, the lattice is disturbed to a triclinic one considering differences that are not negligible but minor. Thereby, merely the critical elastic constants are plotted in Figure 2a nized the elastic constant combinations at all available Sb contents in our model. The crystalline structure belongs to the trigonal system, possessing eight independent components: C11, C12, C13, C33, C14, C15, C45, and C44. Owing to the distortion of asymmetry, the lattice is disturbed to a triclinic one considering differences that are not negligible but minor. Thereby, merely the critical elastic constants are plotted in Figure 2a. Across the board, the elastic constants possess a positive correlation function of Sb content x. Overall trends are linear except for the proportion of 0.25 and 0.375, which exhibit abnormal bias. The comprehensive elastic properties, such as isotropic estimated bulk modulus and shear modulus, can be derived from the elastic constant combination based on the empirical relationship, which is described as: The comprehensive elastic properties, such as isotropic estimated bulk modulus and shear modulus, can be derived from the elastic constant combination based on the empirical relationship, which is described as: where Voigt-Reuss-Hill approximation [42][43][44] of the bulk modulus is given as the function of elastic constant C ij and its compliance S ij : and that of shear modulus is:  [45]. The bulk modulus and shear modulus, estimated from the elastic constants combinations as well, are 42.0 GPa and 25.9 GPa, respectively. It has double resistance to shearing but a similar bulk modulus compared to Mg 3 Bi 2-x Sb x . The causa essential for these variances of electronic, thermal, and thermoelectric properties is the arrangement of the electron cloud caused by the different occupation of Bi/Sb. Whatever the x is, the total number of valence electrons does not change for the ternary compound Mg 3 Bi 2-x Sb x . Therefore, from the statistical perspective, average valence electron concentration (VEC) has an inverse relationship with the volume, i.e., VEC is of a higher level for the smaller volume.
The volume of supercell at each Sb/Bi ratio is also given in Figure 2c, which possesses an almost consummate linear declining relation of x. The bulk modulus, considered as the representative of the comprehensive elastic property, is thereby intuitively induced to possess a linear degressive function of x. The real tested bulk modulus curve with average VEC, represented by the reciprocal of the structural volume, is plotted in Figure 2d. Same with the trend of elastic constants, anomalous low bulk modulus appears at the VEC in the vicinity of x = 0.5, the special one with prominent thermoelectricity, likely associated with this phenomenon. In this manner, we can consciously enlarge the VEC by doping some other ions to tune up the properties, such as improving the bulk stiffness. This also explained the better thermoelectric performance observed around x = 0.5 after Te and Co doping [29,46], and some similar method can be further applied to the design of functional materials.
In the like manner of Mg 3 Bi 1.5 Sb 0.5 , we investigated the stiffness property ( Figure 3) and Poisson's effect (Figure 4) in all directions and all available Sb content x of 40-atom unitcells.     The evolution of direction-dependent Young's modulus in the three-dimensional space with the growth of Sb's proportion is shown in Figure 3a-o. For each proportion, Young's modulus exhibits a distinctive anisotropy where the maximum value is about three times larger than the minimum. Distribution in space is like a cubic shape, whose corner occurs along the x 3 direction, which denotes the stiffest property. Similarly, the cubic structure is stretched along x 3 direction. Six planes of the cubic are concave faces because of the low rigidity, whereas, vice versa, eight corners are protuberant. The lattice vector in the x 1 -x 2 plane is across the midpoints of the edges in this plane where the trigonal symmetry still exists. These distribution characteristics arise from similar lattice structures where the Mg2 ions take the 3 m and −3 m symmetry sites of the original generated structure.
A smooth transition is displayed from purely Bi ions occupied to purely Sb ions occupied with the increasing value of the Sb ions occupying ratio from the Young's modulus along x 3 in Figure 3p. In general, Mg 3 Bi 2 and Mg 3 Sb 2 exhibit similar Young's modulus distribution features except for the specific value at eight corners. The maxima value of Young's modulus grows from 55 GPa to 68 GPa, while the minima value has little change with a fluctuation within the range from 23 to 25 GPa. It is reasonable that Bi and Sb have a similar outer layer valence electron distribution but different atom sizes, which is responsible for the interionic distance. For instance, the lattice parameter of Mg 3 Bi 2 along the x 1 direction is 9.45 Å while that of Mg 3 Sb 2 is 9.20 Å and that along the x 3 are 14.8 Å and 14.5 Å, respectively. Naturally, larger stiffness is owing to the shorter atomic distance, where higher electronic charge density exists, leading to the stronger interaction between ions.
Poisson's ratio describes the reverse transverse deformation caused by mechanical loading. For Mg 3 Bi 2-x Sb x materials, there is no obvious change relationship from the 3D figures with the increase of x, which is different from Young's modulus and shear modulus. This is because Poisson's ratio is related to the relative value of the elastic constant, while each component of the elastic constant has an approximately linear increasing trend. It is worth noting that Poisson's ratio for Mg 3 Bi 2-x Sb x varies widely for different angles. For the plane vertical to the x 3 direction, the Poisson's ratios of all directions on it are small, minima of which are close to 0 or even negative, where the transverse shrinking deformation does not change with the tensile strain, which is rare in conventional alloy materials. In contrast, for all directions within the x 1 -x 2 plane, the Poisson's effect spans a lot, although the minimum values are also close to zero.
For six trigonal symmetry directions, two of which are the lattice vectors, exhibit the largest Poisson's ratio of 0.8. To conclude, as shown in Figure 4p, with the increase of Sb content, the Poisson's effect does not obtain an obvious improvement except for the proportion of 0.25 and 0.375, as well as other mechanical properties. The innate physical images between the prominent thermoelectricity and well-characterized elastic properties at around 0.5 are still elusive and pending issues.
The Vicker's hardness can be estimated by empirical models from Jiang's model [47], Teter's model [48], or Chen's model [49,50], expressed by: and the universal Young's modulus E is given by E = 9GB/(3B + G), whose results are illustrated in Figure 5a. The negative values of the two materials (x = 0 and x = 0.125) predicted from the H v and H G are unphysical and, therefore, not shown here. It implies that the empirical relationship of H v and H G is improper for these compounds. Predictions of Jiang's models show that the largest VEC has the largest hardness and vice versa. Bi 2 Te 3 has a much higher hardness of 399 HV as a reference.
Materials 2022, 15, x FOR PEER REVIEW 9 of 12 that the empirical relationship of and is improper for these compounds. Predictions of Jiang's models show that the largest VEC has the largest hardness and vice versa. Bi2Te3 has a much higher hardness of 399 HV as a reference. Similar to hardness, ductility could be presumably predicted by the ratio between bulk and shear modulus / , which is named as Pugh ratio [51]. It is empirically induced that the crystalline is ductile if the Pugh ratio is larger than 1.75 (dashed line in Figure 5b), while the Puff ratio of Bi2Te3 is 1.62, which is predicted to be brittle in contrast. As a general trend, excluding that of x around 0.5 in Figure 5b, the ductility of Mg3Bi2-xSbx exhibits a Similar to hardness, ductility could be presumably predicted by the ratio between bulk and shear modulus B/G, which is named as Pugh ratio [51]. It is empirically induced that the crystalline is ductile if the Pugh ratio is larger than 1.75 (dashed line in Figure 5b), while the Puff ratio of Bi 2 Te 3 is 1.62, which is predicted to be brittle in contrast. As a general trend, excluding that of x around 0.5 in Figure 5b, the ductility of Mg 3 Bi 2-x Sb x exhibits a declining dependence of Sb content x since the Pugh ratio is positively corrected to ductility. Mg 3 Bi 2-x Sb x behaves in a more brittle manner for x = 0.25 than that of Mg 3 Sb 2 .
The anisotropicity can be described by universal anisotropies A U , shear anisotropies A G and bulk anisotropies A B [52,53], whose expression is: which quantitatively provides insights into the orientation-dependence anisotropic properties. As shown in Figure 5c, the calculated A B is close to 0. Shear anisotropies A G are much larger for Mg 3 Bi 2-x Sb x, which reveals that the bulk properties are isotropic and shear resistances are on the contrary, which agrees with the aforementioned orientation-dependent shear modulus in Figure 1d; however, Bi 2 Te 3 does not, whose A B is 1.58, larger than 1 as well, declaring its stronger anisotropy. The universal anisotropies A U suggest the strong directional dependence of Mg 3 Bi 2-x Sb x compared to other alloys such as Zr 5 Sn 3 and Zr 5 Sn 3 X (X = B, Nb, and Sn), whose holistic A U are much less than 0.5 [53]. The anisotropy is strengthened with more Sb content whose A U increases slightly with respect to an increase in VEC. Debye temperature, an important evaluation of thermodynamic properties, is calculated by the average sound velocity, which can be estimated by elastic constants [54]. Proportional to average sound velocity v m , the Debye temperature Θ is estimated by: where the h, N, and k B are the Planck constant, Avogadro constant, and Boltzmann constant, respectively, ρ is the mass density, q is the number of atoms in the supercell, and M is the molecular weight of all atoms in the supercell. The average sound velocity is proposed by shear sound velocity v s and longitudinal sound velocity v l : where v s = G/ρ and v l = B + 3 4 G /ρ. It is obtained from Figure 5d that Debye temperature grows from the minimum 16.8K of Mg 3 Bi 2 to the maximum 22.7 K of Mg 3 Sb 2 , similar to that of Bi 2 Te 3 , which is predicted to be 23.2 K. Debye temperature possess a positive function of VEC indicating that the melting point is higher with tighter bonding.

Conclusions
We have systematically investigated the elasticity of Mg 3 Bi 2-x Sb x (0 ≤ x ≤ 2) by means of DFT calculations. For the same concentration, thus the same chemical formula and stoichiometry, the atomistic configurations are different at unicells that differ in size. The distinct configurations can be viewed directly from the structures and electron location function distributions for Mg 3 Bi 1.5 Sb 0.5 with 10 atoms and 40 atoms. The structure has impacts on the properties, and we have examined the effect of the unitcell size on the elasticity of the Mg 3 Bi 2-x Sb x , taking Mg 3 Bi 1.5 Sb 0.5 as an example. Anisotropic elastic properties are discovered from the three-dimensional visualization of modules, whose maxima values are about three times larger than the minima, and Poisson's ratios possess a wide variation range with different orientations.
The elasticity variety of the Sb/Bi ratio is also systematically explored for the ternary compound Mg 3 Bi 2-x Sb x . From the single elastic constant component to comprehensive bulk modulus, the elastic properties exhibit a linear dependency of the valence electric concentration VEC, except for that of x around 0.5, which is reported to have the merit zT 1.5 at room temperature. The distribution of Young's modulus and Poisson's ratio in the three-dimensional Cartesian coordinates of all available x are analyzed, as well as other estimated material properties, including hardness, ductility, anisotropicity, and Debye temperature. The abnormal distribution of Poisson's ratio at x = 0.25 and x = 0.375 can be visualized directly. These phenomena may be of great importance in finding the underling physical rules for Mg 3 Bi 1.5 Sb 0.5, which might give important references for the further discovery of novel functional materials with higher zT.
Author Contributions: Q.P. and X.-J.C. conceived the project; S.Z., X.Y. and Q.P. performed the first-principle-based calculation and analysis. All authors analyzed data, discussed the results, wrote the draft, and reviewed and edited the manuscript. All authors have read and agreed to the published version of the manuscript.