Understanding the complex electronic structure of Mg₃Sb₂ and the effect of alloying through first-principles tight-binding models

Abstract

Understanding the electronic origin of good thermoelectric (TE) transport properties is a crucial step to design and discover new TE materials. Mg₃Sb₂ is an intensely researched thermoelectric material, the high TE performance of which comes from six-fold degenerate conduction band minimum at low symmetry k points, leading simultaneously to a high Seebeck coefficient and a good carrier mobility. Furthermore, upon alloying, its low-lying conduction band structure changes, which improves or deteriorates the TE transport properties. Although extensively studied experimentally, the electronic origin of these changes is not well understood. Based on the first-principles derived Wannier function tight-binding analysis of both pristine Mg₃Sb₂ and its Bi/Ca alloyed derivatives, we reveal that the simple Mg s–Sb p orbital interactions are not adequate and sufficient to describe the electronic structure of Mg₃Sb₂. Inclusion of both higher energy Mg p orbitals and lower energy Sb s orbitals is demonstrated to be important to understand the shape of the lowest conduction band and its change upon alloying. The systematic chemical understanding developed in this work could enable us to better control the conduction band dispersion by manipulating chemical interactions in Mg₃Sb₂.

---

1. Introduction

Thermoelectric (TE) materials, with their ability to convert waste or environment heat into electricity, can help to bring a carbon neutral society.^1–4 The thermoelectric conversion efficiency can be determined by materials' figure of merit zT = T × S²σ/κ where κ is the total thermal conductivity that includes both the lattice and electronic contributions, S is the Seebeck coefficient and σ the electrical conductivity.^5,6 The term S²σ is governed by the material's electronic band dispersion. Compounds with electronic structures showing light carrier pockets as well as symmetry induced band degeneracy typically exhibit higher TE performance.⁷ Therefore, understanding the electronic structures of thermoelectric compounds, as well as their chemical bonding origin, is of crucial importance for finding new TE compounds or optimizing existing TE ones.^8,9

Among the best thermoelectric materials, Mg₃Sb₂ and its alloyed compounds exhibit relatively high thermoelectric performance in the middle temperature region.^10–12 For instance, the state-of-the-art n-type Mg₃Sb₂ alloyed with Bi can reach zT of above 1.0 over a wide temperature range. As a result, experimental modules utilizing Mg₃Sb₂ as one of the thermoelectric elements shows a promising conversion efficiency of 7.3% with a temperature difference of 315 K.¹³ Indeed, its high TE performance can be attributed to the fact that its light conduction band minimum (CBM) is situated at a low symmetry k point in the Brillouin zone (BZ), leading to six-fold-degenerate carrier pockets that lead to high Seebeck coefficients while maintaining good carrier mobility.^14–16 Furthermore, alloying with Bi can further reduce the conduction band effective mass at the CBM point and increase band degeneracy by reducing the conduction band edge energy at the Γ point to convergence. At the same time, the lattice thermal conductivity can be effectively reduced with heavier Bi atoms, further improving the thermoelectric performance of Mg₃(Sb,Bi)₂.^12,17–19

Despite the progress in improving and tuning Mg₃Sb₂'s thermoelectric properties, its electronic structure, in relation to the orbital interactions between Mg and Sb atoms, is not so well understood due to the material's complex crystal structure. Zhang et al. suggested that the interactions between Mg and Sb atoms are mostly ionic in nature based on the evidence of electron density topology analysis.²⁰ Their calculated atomic charges, based on Bader partitioning,^21,22 are relatively close to the formal charge of Mg²⁺ and Sb³⁻ in Mg₃Sb₂ (+1.5 on Mg and −2.2 on Sb) and support strong ionicity in the compound. However, Wang et al. showed that covalent interactions between Mg and Sb are also important, with electron localization function (ELF) basins observed between Mg and Sb atoms (a fingerprint for covalent bonding). It turns out, from their calculations, that an integrated crystal orbital Hamilton population (ICOHP) value over Mg and all its nearest neighbor Sb atoms is around −6 eV, further supporting their covalent interpretation.²³

Some previous research studied the effect of Bi alloying on the electronic structure of the lowest conduction band.^12,18,24,25 In a review by Li et al., Bi alloying was depicted to slightly decrease the term E^Mg_s − E^anion_p and subsequently lead to the rapid lowering of conduction band minimum at the Γ point.²⁴ However, the atomic orbital energy difference between Bi 6p and Sb 5p orbitals is only around 0.2 eV, which is unlikely to lead to a dramatic change of the conduction band. A later study by Pan et al. hypothesized that the Mg 3s bands become “wider” upon Bi alloying, therefore decreasing or even closing the band gap at Γ.²⁶ However, it is doubtful that such a view, without firm evidence, corresponds to the actual reason for the closing of the direct gap at Γ.

Other studies tried to explain the location of the conduction band minimum at the BZ edge in Mg₃Sb₂. Sun et al. observed some bonding character between the two non-equivalent Mg atoms and suggested that the conduction band minimum at a low symmetry point is caused by Mg–Mg interactions.²⁷ Later, Wang et al., when comparing Mg₃Sb₂ and CaMg₂Sb₂, suggested that the absence of Ca (3d shell) unoccupied orbitals in the conduction band drives CBM from the high symmetry M point to the low symmetry CBM point in Mg₃Sb₂ (we use CBM to indicate the k point of the conduction band minimum in the first BZ).²⁸ However, they somehow ignored the existence of unoccupied Mg 3p character, which can play a similar role. A final work to mention is that from Brod et al., which connects the conduction band minimum to the Mg–Sb nearest neighbor antibonding interactions based on the evidence from band resolved COHP between Mg 3s and Sb 5p orbitals.²⁹ The weakness of their analysis is that the model they proposed is only qualitative and is restricted to consider only s–p interactions. Since the lower conduction bands, mainly made of Mg 3s orbitals, cannot be separated from the higher energy conduction bands, understanding of the conduction band structure based only on Mg s and Sb p orbitals is not sufficient.

In this work, we use a first-principles based tight-binding model to provide a quantitative analysis of the orbital interactions in Mg₃Sb₂. The effect of alloying is also captured and compared using the same tight-binding framework within the virtual crystal approximation. In contrast to previous understandings, we show that the lower energy Sb 5s orbitals and higher energy Mg 3p orbitals play an important role in establishing the electronic structure of Mg₃Sb₂, especially the lowest conduction bands mainly of Mg 3s character and where the Mg 3s and Sb 5p interactions only play a minor role. In addition, our analysis also consistently explains the effect of alloying in Mg₃Sb₂ from the changes of orbital interactions.

2. Methods

2.1. Electronic structure calculations

The electronic structure of Mg₃Sb₂ was calculated self-consistently using the density functional theory (DFT) method using plane waves as implemented in the Quantum Espresso (QE) package.³⁰ Norm-conserving pseudopotentials in the form of generalized gradient approximation (PBE-GGA) used in the calculations were obtained from the Pseudo-dojo website.^31,32 The energy cutoff for the plane wave expansion was set to 80 Ry and the cutoff of charge density was four times higher. Fixed occupation with a k-grid of 14 × 14 × 7 was used for the Brillouin zone summation of density.³³ The lattice parameters and atomic coordinates were fully relaxed (forces less than 1 × 10⁻³ Ry per au). The relaxed cell parameters a and c of 4.595 and 7.272 Å, respectively, are slightly overestimated with respect to the experimentally measured values (a = 4.573 and c = 70.229 Å) for Mg₃Sb₂.³⁴ The calculated indirect electronic band gap is 0.168 eV, smaller than the experimentally measured value (>0.4 eV) and that calculated with TB-mBJ functionals.^10,35 The valence band (VB) maximum is located at the Γ point and the conduction band minimum is located at k = (0.425, 0.0, 0.333). The coordinate of the CBM point agrees very well with previous calculations.³⁵ The calculated secondary conduction band minimum at the K point is 0.2 eV above the conduction band minimum. Although the band gap is underestimated, important features of the electronic structure are well produced by the PBE-GGA calculations, including the trend of alloying. As a result, we do not expect that the underestimation of the band gap undermines the discussion here. For simplification, the spin orbital coupling (SOC) effect was not considered in the calculations. The effect of SOC in Mg₃Sb₂ was studied by Zhang and Iversen.¹⁹ It was found that the SOC splits the electronic bands near the VB edge at the Γ point, which mainly derive from Sb p orbitals, but its effect on the conduction band is negligible. Since our study does not focus on the quantitative aspect of the VB edge, it should be safe to omit SOC in the DFT calculations and subsequent Wannier tight-binding model.

Virtual crystal approximation (VCA), in which an averaged pseudopotential approximates the effect of “alloying” of a specific site, was used for the calculation of the electronic structure of alloyed compounds.³⁶ In this way, we were able to compare the Wannier tight-binding model of pristine Mg₃Sb₂ and its alloyed compounds. We notice that the limitation of VCA is that the effects of local atomic configuration/relaxation are neglected. (For example, both Mg–Sb and Mg–Bi bonds with different distances exist in actual Bi alloyed Mg₃Sb₂, but not in VCA calculations.) The averaged pseudopotential was obtained using the virtual.x code that comes with the QE package. All structures using VCA pseudopotentials were fully relaxed within the same crystal symmetry as that of Mg₃Sb₂. The relaxed lattice parameters from VCA calculations are shown in Fig. S1.† Bader analysis was performed from the computed electron density using the Bader code.³⁷ High symmetry points in the Brillouin zone of space group P3̄m₁ are shown in Fig. 2.³⁸

2.2. Wannier tight-binding model

Orthogonal Wannier functions (that resemble atomic orbitals) were used as a basis to obtain a first-principles tight-binding model for Mg₃Sb₂ and its alloys, as implemented in the Wannier90 code.³⁹ It must be noted that Wannier functions used here are not maximally localized, but close to the set of initial atomic orbitals apart from orthogonalization.⁴⁰ Thus, we denote these Wannier functions using the initial atomic orbitals. The real space Hamiltonian operator was obtained by projecting eight valence bands and twelve disentangled conduction bands to the localized orbitals of Sb and Mg s/p orbital functions.⁴¹ All states lying 4 eV above the CBM were preserved exactly in the Wannier projection. For alloyed Mg₃Sb₂ derivatives calculated using VCA, the exact same projection scheme and initial orbital functions were used to obtain their corresponding tight-binding Hamiltonian operator. The band structures obtained from the Wannier tight-binding model, compared to the DFT calculated ones, are shown in Fig. S2.† More details about the Wannier projections are given in Appendix A.

3. Results

3.1. Electronic structure of Mg₃Sb₂

The structure of Mg₃Sb₂, which crystalizes in the layered CaAl₂Si₂ structure type with space group P3̄m₁ (no. 164), is shown in Fig. 1a. The Mg atoms occupy two different Wyckoff positions: interlayer Mg1 atoms are located at the 1a site and are octahedrally coordinated with six neighboring Sb atoms (at 3.12 Å apart), and intralayer Mg2 atoms are located at the 2d site (1/3, 2/3, z) in a tetrahedral environment somewhat elongated along the z direction. Surprisingly, although the crystal structure features a plane of Mg separating edge-sharing Mg2–Sb₄ tetrahedra, the two Mg atoms do not differ much chemically, as evidenced from the computed atom-in-molecular (AIM) Bader charges and the ICOHP values reported previously.^20,23

Fig. 1 Crystal structure and electronic structures of Mg₃Sb₂. (a) Crystal structure of Mg₃Sb₂ with the two non-equivalent Mg sites shown in their octahedral and tetrahedral coordination, and (b) calculated electronic band structure of Mg₃Sb₂ and the projected density of states (pDOS). The black lines in the band structure show the DFT-calculated band structure while the red line (visible in the upper conduction band) show the band structure obtained with the Wannier tight binding model.

Fig. 2 Upper panel: Brillouin zone of Mg₃Sb₂ with labeled special symmetry points. Lower panel: energy dispersion of the Bloch functions constructed from the different Wannier orbitals as a function of the wave vector k. The dotted lines are the three lowest conduction bands of Mg₃Sb₂. The zero energy is set to be the calculated energy of the conduction band minimum in Mg₃Sb₂.

An analysis of the electronic structure (Fig. 1b) indicates that the occupied states, from the two bands at −10 eV up to the Fermi level, are mainly of Sb s and p character. However, it should be noted, from the tight-binding projected density of states (pDOS), that Mg orbitals, both 3s and higher energy 3p orbitals, somewhat mix into the valence band region. Above the Fermi level, the conduction bands derive mainly from the unfilled Mg s and p orbitals with some mixing with Sb p orbitals. Due to the substantial orbital mixing both in valence and conduction bands, it can be reasoned that covalent interactions between Mg and Sb atoms are important for binding in Mg₃Sb₂. Note that looking at electron density alone, as done previously by Zhang et al., may wrongly lead to the neglection of such a covalent character.^20,42

Through tight-binding parameterization, we can gain important insights of the orbital interactions in Mg₃Sb₂ that allow us to evaluate previous assumptions about the bonding nature of Mg₃Sb₂ and build a solid ground for discussing the effect of alloying, i.e., substituting a large part of Mg or Sb. First, through the Wannier function tight-binding model, we can recognize the importance of long-range interactions that control the electronic structure in Mg₃Sb₂. In contrast to high symmetry cubic PbTe or perovskites the electronic structure of which can be reproduced well by only considering the first or secondary neighbor interactions,^43,44 interactions up to 10 Å are required to obtain a fair description of the conduction band dispersion in Mg₃Sb₂ (shown in Fig. S3 and S4†). As a result, understanding of the role of orbital interactions in Mg₃Sb₂ and its alloyed derivatives is required to separate the effects of long and short-ranged chemical interactions.

According to perturbation theory, the “ionic” interaction between two orbitals typically occurs when their energies strongly differ. For example, in NaCl, the atomic orbital energy of Na 3s is more than 8 eV higher in energy than the Cl 3p atomic orbitals, in their respective elemental form, causing electron transfer from the former to the latter and consequently inducing strong ionic character when the solid compound is formed. In contrast, the energy of the Mg s atomic orbital is only around 0.8 eV higher than that of Sb p orbitals in their elemental form, suggesting that the Mg s atomic orbital is more likely to interact with Sb p orbitals in other ways.⁴⁵ In solids, due to the periodicity and the modulation of the wave vector k, it is more meaningful to discuss Bloch functions, constructed from the localized, atomic-centered Wannier functions by:

The energies of these Bloch functions thus include both the “on-site” energy of the atomic orbitals and orbital interactions, modulated by k, between their orbital sublattices.^40,46 Surprisingly, in Mg₃Sb₂, the tight-binding parameterization reveals that the Bloch functions constructed from the Sb p orbitals lie above those from Mg s orbitals in energy in a large portion of the BZ, as shown in Fig. 2. Therefore, only considering the Mg s and Sb p orbitals and their interactions cannot explain the electronic structure in Mg₃Sb₂ where the valence bands are dominantly of Sb p orbital character and the conduction band of Mg s orbital character.

3.2. Orbital interactions

Fig. 3 illustrates the interaction diagrams between Bloch orbitals at different high-symmetry k points in the Brillouin zone obtained from the tight-binding calculations. These diagrams allow us to construct the electronic band structure in terms of interactions between different Bloch functions. At different k points, these interactions are determined by:
and are different due to the complex phase factor e^ikR that determines the accumulation and cancellation of interactions.^45,47,48 The orbital interactions in this form also include long-range terms due to the summation over R. The band energies ε_ik can then be obtained after diagonalization of the Hamiltonian matrix on the basis of periodic Bloch functions.

Fig. 3 Interaction diagrams between Bloch functions in Mg₃Sb₂ at different k points in the Brillouin zone. The transparency of the interaction lines is proportional to the coefficient of the orbital component in the final energy eigenstates. The resulting valence and conduction bands are colored in black and red, respectively.

To understand why the valence band is mainly from the Sb p character, it is important to notice, from Fig. 3, that the Sb p bands interact more significantly with the higher energy Mg p bands than with the Mg s bands. It is thus clear that the polar covalent interaction between Mg p and Sb p drives down the energy of the Sb p bands to dominantly form the valence band, leading to the band gap that separates the valence bands from the lowest conduction band (mainly from the Mg s contribution). Strong Mg p and Sb p interactions are reflected by the mixing of orbitals in the projected DOS in Fig. 1b, as well as from the previous COHP calculations.²³ Furthermore, in Fig. S5,† we show that the main valence and conduction band structure is approximately correct even if only Mg p, Sb p orbitals and their interactions are considered in the tight-binding model.

For n-type thermoelectric transport properties, the band structure of the lowest conduction band is the most important. The same band also exhibits complex changes upon alloying, indicating peculiar chemical interactions. Interestingly, orbital interaction diagrams in Fig. 3 reveal that the lowest conduction band interact weakly with the orbital components that form the valence band. Therefore, the energy of the conduction band edge is somewhat close to the energy of the Mg s Bloch orbitals (Fig. 2) with dominant Mg s contribution. Additional information regarding orbital interactions can be provided by the summed, off-site COHP curves which reveal the overall bonding or antibonding nature of the interactions (defined in Appendix B). In Fig. 4a, the dispersion of the lowest conduction band in Mg₃Sb₂ together with the off-site COHP and orbital components is shown. Furthermore, in Fig. 5, S6 and S7,† the entire COHP matrix at three k points for the lowest conduction band is shown in the form of a heatmap which provides the sign and magnitude of orbital-wise COHPs. Based on these results, the interpretation of the bonding nature of the lowest conduction in Mg₃Sb₂ can be made.

Fig. 4 Band dispersion (top), off-site COHP (middle) and orbital component (bottom) of the lowest conduction band in (a) Mg₃Sb₂, (b) Bi-alloyed Mg₃Sb₂ on Sb site, and (c) Ca-alloyed Mg₃Sb₂ on Mg1 site. The alloyed systems correspond to 50% substitution of the corresponding sites. In the bottom figure, Mg s and p indicate the combined components from all three Mg atoms present in the unit cell. Note that Mg s/p includes both Mg1 s/p and Mg2 s/p orbitals in the orbital component figure.

Fig. 5 Off-diagonal (orbital pair-wise) COHP heatmap for the lowest conduction band at the CBM point. The diagonal elements (on-site COHP) are not shown in the heatmap. Negative values of COHP indicate bonding character while positive values indicate antibonding character.

At the Γ point (and similarly at A and A′ points), the lower Sb s orbital plays an important role. At these points, the orbital component of the lowest conduction band from the Sb s contribution is large and the off-site COHP curve reveals an overall antibonding character (also see Fig. S6†). The strong s–s interaction can be understood since at Γ, all the Sb s orbitals are in phase with the Mg s orbitals while the s–p interactions are cancelled. As a result, the dominant orbital interaction at Γ can be considered to be that of a two-orbital problem involving the Mg and Sb s orbitals.

In Mg₃Sb₂, the conduction band minimum at the K point (1/3, 1/3, 0) is very close to the true conduction band minimum and can be considered to be nearly degenerate at medium temperature in terms of electron transport.³⁵ From Fig. 3 and 4a, it is found that the conduction band minimum state at K only contains Mg1 s and Mg2 p orbital contributions. Although the large component of higher energy Mg p orbitals at K would lead to higher band energy, this conduction band minimum is stabilized by the interaction between the different Mg p orbitals as shown in Fig. S7.† As a result, the conduction band minimum energy at K is located only slightly higher than the energy of Mg1 s Bloch orbitals (Fig. 2).

The true CBM is located somewhere near the L point (1/2, 0, 1/2) and M point (1/2, 0, 0). At L, M and the CBM points, the conduction band minimum states are dominantly of Mg s character (see Fig. S8† for the real space electronic density corresponding to the carrier pockets of CBM). Although other orbitals participate in the interaction at the conduction band minimum state, the net interaction COHP is very small, according to the results shown Fig. 4a. We argue that the small net COHP of CBM is a result of the cancellation of weak bonding/antibonding interaction between Mg s and other higher/lower energy orbitals (as in a typical three-orbital system in a molecule according to perturbation theory which leads to a non-bonding intermediate molecular orbital). This can be seen from the pair-wise COHP heatmap shown in Fig. 5 where both positive and negative off-site COHPs are presented. It is important to note that the small off-site COHP does not only occur at the CBM point, but also in the vicinity of the CBM point in the BZ. Since the change in COHP in the vicinity of the CBM point is much smaller compared to the change in energy of the conduction band, the covalent interactions (as revealed by off-site COHP) do not directly lead to band dispersion. Therefore, our interpretation differs from the mechanism suggested by Brod et al. that attribute band dispersion near CBM mainly to the antibonding interaction between the Mg s and Sb p orbitals.²⁹ In contrast to the overall small off-site COHP values, the ratio of orbital components changes significantly in the vicinity of the CBM point, suggesting that the band dispersion near CBM is due to the change in on-site COHP (see Appendix B). At the CBM point, the lower energy Mg s orbital reaches its largest orbital components, as shown in Fig. 4a. Moving away from the band minimum, the component of higher energy Mg p orbitals increases significantly. This shift in orbital components from lower energy orbitals to higher energy orbitals leads to band dispersion (see Fig. S9† for the on-site COHP curve, which follows the band energy closely).

Even though the antibonding interactions between the nearest neighbors Mg and Sb do not lead directly to the band dispersion, we found that they are likely to play an important role in determining the position of the CBM point. Wang et al. suggested that the position of the CBM point is a result of a competition of interactions between Sb and Mg s orbitals at different sites.²⁸ If this is indeed the case, then it can be expected that alloying would shift this conduction band minimum. However, during our analysis, we observed that the position of the CBM point, approximately at (0.42, 0, 0.33) in the BZ, is relatively stable upon alloying. Although alloying with Ca on Mg1 site shifts the conduction band minimum to the M point, a local minimum persists even at 50% substitution (see Fig. 6 and Table S1† for the coordinates of local minima for all the alloyed compounds). Furthermore, truncating the tight-binding Hamiltonian operator reveals that the local minimum can be observed when only nearest neighbor interactions are considered, even though the energy of the lowest conduction band is not systematically correct (Fig. S4†). The fact that the position of the CBM minimum is robust with changing orbital interaction values implies that the position of this minimum is related to the unique local crystal structure configuration of Mg₃Sb₂ and nearest neighbor Mg–Sb interactions.

Fig. 6 Band structure of Mg₃Sb₂ alloyed with Bi or Ca at different non-equivalent sites. Bands of pristine Mg₃Sb₂ are shown in brown while the bands of the alloyed compounds are shown in green. Top figures show band structures calculated using the VCA approximation. The bottom figures show the band structures calculated using the Wannier tight-binding model that incorporates the effect of alloying by changes in the on-site energies and the nearest neighbor interactions. A close-up view is shown for the Bi doped band structures.

In summary, the valence and conduction band structure of Mg₃Sb₂ is more complex than that expected simply from Mg s–Sb p interactions. Indeed, the Sb p orbitals mostly form the valence band with bonding interactions from both Mg s and Mg p orbitals. On the other hand, the lowest conduction band structure is determined as much by Sb s and Mg p as by Sb p orbitals. Near the CBM point, the lowest conduction band is almost non-bonding with a rather flat COHP curve, and the band energy reaches a minimum when the orbital component of Mg s reaches a maximum value.

3.3. Understanding the effect of alloying

We proceed to show how the effect of alloying at different k points can be understood based on the corresponding orbital interaction diagrams given by the tight-binding analysis. The effect of alloying (up to 50% substitution) is approximated using VCA in our DFT calculations and the resulting electronic structures are shown at the top of Fig. 6. The band structure changes by alloying from VCA calculations agree well with the calculations of the end point compositions performed in previous studies.^20,24,28

Previously, the effect of alloying was thought to be the result of disruption of the nearest neighbor interactions.²⁹ However, since it is observed (Fig. S3 and S4†) that long-range interactions, at least beyond nearest neighbors, are important to produce a quantitative tight-binding electronic structure in Mg₃Sb₂, it should be investigated whether the alloying effect can be understood by changes of the nearest neighbor interactions alone. To this end, we combined the tight-binding parameters from both pristine Mg₃Sb₂ and its alloyed compounds so that in the tight-binding Hamiltonian, the nearest neighbor interaction parameters (including the on-site energy) are those of the alloyed systems while any interactions beyond the nearest neighbors are kept as that of pristine Mg₃Sb₂. For simplicity, we refer to these combined tight-binding models “short-range substituted tight-binding” (SR-TB) and the band structure calculated based on these combined tight-binding models are shown in the bottom of Fig. 6. The changes observed in DFT calculations that are reproduced by the SR-TB band structures can be attributed to changes of the local, nearest neighbor interactions while those that cannot be reproduced should be attributed to changes of the relatively long-range interactions. Thus, through the SR-TB model, we can attribute the band structure changes of alloying to changes in the nearest-neighbor interactions or long-range interactions beyond the nearest neighbors.

We observed that the following band structure changes can be understood by changes of nearest neighbor interactions upon alloying: (i) the changes in valence band energies, (ii) the rapid closing of the direct band gap at the Γ point and the relatively fixed energies near the CBM and the K points upon Bi alloying, (iii) the upshift of the lowest conduction band except at the Γ point for Ca substituting Mg1 site, and (iv) the shift of the conduction band minimum to the three-fold degenerate M point also for Ca on Mg1 site. We show in the following that these changes can be consistently understood from peculiar orbital interactions in Mg₃Sb₂.

However, some other changes in electronic structure upon alloying cannot be attributed only to nearest neighbor interactions. From the VCA calculations, we see that the band energy at M and M′ moves up slightly upon Bi alloying while the energy at the CBM point goes down, leading to a more dispersive band at CBM. The change in conduction band effective mass upon Bi alloying in Mg₃Sb₂ is known from previous calculations and experimental observations,^19,26 but it is not reproduced by changes in short range interactions, as shown by the enlarged band dispersion in the middle panel of Fig. 6. Previously, Brod et al. argued that the disruption of nearest neighbor Mg–Sb interactions was responsible for the more dispersive conduction band minimum.²⁹ However, our quantitative result reveals that the increased conduction band dispersion is related to change of interactions beyond the nearest neighbors, perhaps due to the more delocalized nature of the Bi p orbitals compared to that of Sb. We also note that the band structure changes introduced by Ca substitution on Mg2 site also cannot be understood by just considering changes in short range interactions (see Fig. S10†). The reason could similarly be a more delocalized p orbital behavior within the covalent Mg–Sb network.

Model orbital interaction diagrams shown in Fig. 7 can, in turn, serve as a guide to understand the effect of alloying for the lowest conduction bands at important k points. These simplified orbital interaction diagrams are derived partly from the orbital interaction diagram of Mg₃Sb₂ shown in Fig. 3, partly from the COHP results based on the tight-binding model. We first discuss the case of Bi alloying. At the Γ point, the Mg s and Sb s antibonding interactions largely determine the energy of the lowest conduction band. Upon Bi alloying, the anion s energy is reduced most notably (by ∼1 eV upon 50% of Bi alloying at the anion site, Table S2†). The lowering of the anion s energy substantially reduces the energy of the antibonding conduction band minimum at Γ according to the two-orbital interaction scheme (for orthogonal orbitals, ε⁰_1/2 are the orbital energies before interaction, V₁₂ is the interaction parameter and ε_1/2 are the eigen-energies from the two-orbital interaction):

Fig. 7 Simplified pictures of orbital interactions in Mg₃Sb₂ at important k points for the lowest conduction band, in terms of the interactions between Mg s and other orbitals. (Note that Mg s and Mg p interactions in general do not occur at the same site.)

Since the anion orbitals do not participate in the lowest conduction band state at the K point, as shown in Fig. 4 and 7, the orbital components at the K point almost remain the same upon Bi alloying (Fig. 4b). A minor upshift in energy is observed and could be attributed to a slightly decreased bonding overlap between the Mg p orbitals as a result of an increase in the lattice parameter (Fig. S1†). Finally, at the CBM point, there is a downshift in band energy. At this point, multiple weak interactions of Mg s with higher/lower orbitals cancel each other, according to our previous analysis. The lowering of the anion s energy would reduce the antibonding interactions for the CBM states, leading to a slightly lower energy of CBM. These changes described here are consistent with the previous calculations, particularly those depicted by Zhang et al.³⁵

With Ca alloying on Mg1 site, the band gap is increased with the upshift in energy for the lowest conduction band. This is not surprising since Mg s orbitals play an important role in this band. However, the reason for the very minor changes at the Γ point and the shifts of the conduction band minimum from the six-fold degenerate CBM point to the three-fold degenerate M points is worth discussing. The higher energy cation p orbitals are important here. At the Γ point, although the conduction band energy is expected to increase with the increase of the on-site energy of cation s, there are, at the same time, stronger bonding interactions with the higher energy orbitals of cation p due to the reduced energy difference E_p − E_s on the cation. The stronger bonding interaction stabilizes the lowest conduction band states (see Fig. S11†). This result is supported by the fact that the overall antibonding interactions at Γ in pristine Mg₃Sb₂ are reduced significantly for Ca alloyed on Mg1 site, with a much larger Mg p component, as shown in Fig. 4c.

The same bonding interactions with higher energy cation p orbitals also explain the change of conduction band minimum from low symmetry point to M point. Near the CBM point in Mg₃Sb₂, the cation s components reach their maximum, so the conduction band energy increases in correspondence with the change of cation s energy. However, at M point, there is a stronger bonding participation of higher energy cation p orbitals (from different atomic sites), which can be observed when we compare the off-site COHP at M′ point with that at the CBM point, in Fig. 4c. The bonding interactions from higher orbitals drive down the conduction band energy at M′ and M compared to the bonding minimum at the CBM point, shifting the conduction band minimum from the six-fold degenerate CBM point to three-fold degenerate M point. This interpretation of the changes of conduction band also explains why the conduction band minimum does not move to the L point with Ca alloying, since the lowest conduction band states at L (L′) point contain less cation p contribution than at M (M′), as can be seen in Fig. 4c. The above analysis is in line with previous studies, except that previous work considered the Ca 3d shell but does not consider the higher energy orbitals in the pristine case.^28,29

In summary, understanding orbital interactions in Mg₃Sb₂ allowed us to interpret the effect of alloying which affects the lowest conduction band structure whose carrier pockets are important for thermoelectric transport properties. Although it was previously thought that the effect of alloying could be mainly attributed to the Mg s and Sb p nearest neighbor interaction, our tight-binding based analysis demonstrates instead that changes of interactions between Mg s orbitals with higher energy Mg p and lower energy Sb s orbitals are more important. Furthermore, some subtle band structure changes upon alloying may require effects beyond nearest neighbor interactions.

3.4. Strategy for band structure engineering

Let us discuss here the implication of the chemical insights developed by suggesting a band structure engineering strategy to converge the conduction band minimum at the K and CBM points. In Mg₃Sb₂, the lowest conduction band energy at the K point is slightly above the conduction band minimum. According to our analysis, the energy of local conduction band minimum at the K point is related to the bonding interaction between the p orbitals of different Mg atoms. As a result, it can be expected that if the orbital overlap between these Mg p orbitals increases, their bonding interaction will be stronger and will lead to an energy lowering at the K point. This can be achieved by a compressive strain. At the same time, due to the relative weak interaction for the lowest conduction band at the CBM point, its energy should not be affected so much. In Fig. S12,† we show, with first-principles band structure calculations, that such a convergence is indeed realized under a −2% strain. We note that this phenomenon was observed in a previous report by Xie et al.⁴⁹ However, its origin was not explained, and this approach has not been experimentally attempted to the best of our knowledge. A recent work illustrated that strain engineering can tune the valence band and thus the p-type thermoelectric performance in Mg₃Sb₂ films. According to our analysis, the effect of strain on the n-type thermoelectric performance is also worth investigating.

4. Conclusion

In conclusion, our first-principles based tight-binding study provides a detailed analysis of orbital interactions occurring in Mg₃Sb₂. Somewhat in contrast to previous studies which considered the formation of valence and conduction bands with a Mg s–Sb p (ionic) interaction framework, we demonstrated here that the covalent interactions are important with higher energy Mg p orbitals participating importantly in the crystal bonding (valence band) and that both the higher energy Mg p and lower energy Sb s orbitals shape the dispersion of the lowest conduction band. Although the chemical interactions are different, the conduction band minima at the CBM and K points lie very close in energy to the Mg s orbitals, leading to near-degenerate conduction valleys. Important band structure changes upon alloying on Sb or Mg1 site can be understood as changes in the nearest neighbor interactions. However, the change in effective mass of CBM upon Bi alloying should be attributed to changes in the long-range interactions. These chemical interactions determine the electronic structure and correspondingly the electrical transport properties. We hope that the study presented here can provide helpful insights for further improvement of the thermoelectric transport properties in n-type Mg₃Sb₂ by offering a guideline for experimentally manipulating chemical interactions in this system.

Appendix

A. Wannier tight-binding electronic structure calculations

In this section we describe the projection scheme to obtain localized Wannier functions and the details of the tight-binding calculations employed in this work. Given a set of trial orbitals , a set of non-orthogonal Bloch wavefunctions can be obtained by projecting a set of J bands that are the eigenstates of the crystal Hamiltonian operator at k (for clarity, we use k to indicate wave vector in the following):⁴⁰

An orthogonal set of Bloch functions, which are used to obtain the localized orthogonal Wannier functions, can be obtained by the transformation:

where S_k^−1/2 is the overlap matrix between the non-orthogonal functions . The superscript (W) indicates “Wannier” gauge.⁵⁰ Finally, a set of localized Wannier functions can be obtained by the Fourier transformation:
with the reverse transformation:

It should be noted that a disentangle method should be applied to obtain a set of optimally connected subspaces since the higher conduction bands are also considered in the projection. This disentangle method selects a subspace by unitary transformation that minimizes the spread of the wavefunction.^40,41

The electronic structure can be determined from the obtained Wannier function tight-binding Hamiltonian matrix elements (where m and n denote different Wannier functions in the unit cell and R is lattice translation) by Fourier transformation:

Finally, the diagonalization of the Hamiltonian matrix gives the eigenstates at k: , where the coefficients u^k are determined by:⁴⁶

The orbital coefficients that determine the nth orbital components of the eigenstates for ith band are given by c_in(k) = |u_in^k|².

In this work, the trial Wannier functions used for the above-mentioned projection were constructed from the real radial and angular atomic functions. The radial part of the wavefunction is given by:

R(r) = 2α^3/2 exp(−αr)

while the angular part is given by the real spherical harmonic function for s and p orbitals:

B. COHP calculations

To understand the bonding and antibonding character of the electronic bands, the off-site crystal orbital Hamilton population (COHP) was used to measure quantitatively the orbital interaction between atoms. Dronskowski and Blöch defined the COHP between different local atomic centered orbitals.⁵¹ Instead of following their definition, a simpler quantity that measures the orbital interaction between the periodic Bloch functions at k was used in this work, denoted as COHP^ik_nm. The band structure energy ε_ik of the ith band at k can be written as:

Furthermore, “on-site” and “off-site” (non-local) terms can be defined in such a way that ε_ik = COHP^ik_on-site + COHP^ik_off-site.⁵¹ They are given by:

where {n′} is the set of orbital components n′ that is on the same atoms with component n. In the example of Fig. 5, the off-site COHP at the CBM point is the sum of all the orbital pair-wise (off-diagonal) COHP values.

Author contributions

Wenhao Zhang: conceptualization, methodology, analysis, investigation and writing. Jean-François Halet: review and editing, supervision. Takao Mori: review and editing, supervision, funding acquisition.

Conflicts of interest

All authors declare no financial or non-financial competing interests.

Acknowledgements

We acknowledge support from JST Mirai Program JPMJMI19A1. The authors thank Longquan Wang (NIMS & University of Tsukuba) for helpful discussions.

References

1. L. E. Bell, Science, 2008, 321, 1457–1461.
2. T. Mori and S. Priya, MRS Bull., 2018, 43, 176–180.
3. D. Beretta, N. Neophytou, J. M. Hodges, M. G. Kanatzidis, D. Narducci, M. Martin- Gonzalez, M. Beekman, B. Balke, G. Cerretti, W. Tremel, A. Zevalkink, A. I. Hofmann, C. Müller, B. Dörling, M. Campoy-Quiles and M. Caironi, Mater. Sci. Eng., R, 2019, 138, 100501.
4. R. Freer and A. V. Powell, J. Mater. Chem. C, 2020, 8, 441–463.
5. K. Koumoto and T. Mori, Thermoelectric Nanomaterials Materials Design and Applications, Springer-Verlag GmbH, 2015.
6. H. J. Goldsmid, Introduction to Thermoelectricity, Springer Berlin Heidelberg, Berlin, Heidelberg, 2016, vol. 121.
7. Z. M. Gibbs, F. Ricci, G. Li, H. Zhu, K. Persson, G. Ceder, G. Hautier, A. Jain and G. J. Snyder, npj Comput. Mater., 2017, 3, 8.
8. J. Mao, Z. Liu, J. Zhou, H. Zhu, Q. Zhang, G. Chen and Z. Ren, Adv. Phys., 2018, 67, 69–147.
9. T. Hendricks, T. Caillat and T. Mori, Energies, 2022, 15, 7307.
10. J. Zhang, L. Song, A. Mamakhel, M. R. V. Jørgensen and B. B. Iversen, Chem. Mater., 2017, 29, 5371–5383.
11. X. Shi, X. Wang, W. Li and Y. Pei, Small Methods, 2018, 2, 1800022.
12. K. Imasato, S. D. Kang and G. J. Snyder, Energy Environ. Sci., 2019, 12, 965–971.
13. Z. Liu, N. Sato, W. Gao, K. Yubuta, N. Kawamoto, M. Mitome, K. Kurashima, Y. Owada, K. Nagase, C.-H. Lee, J. Yi, K. Tsuchiya and T. Mori, Joule, 2021, 5, 1196–1208.
14. Y. Pei, A. D. LaLonde, H. Wang and G. J. Snyder, Energy Environ. Sci., 2012, 5, 7963.
15. Y. Pei, H. Wang and G. J. Snyder, Adv. Mater., 2012, 24, 6125–6135.
16. Y. Tang, Z. M. Gibbs, L. A. Agapito, G. Li, H.-S. Kim, M. B. Nardelli, S. Curtarolo and G. J. Snyder, Nat. Mater., 2015, 14, 1223–1228.
17. Z. Liu, W. Gao, H. Oshima, K. Nagase, C.-H. Lee and T. Mori, Nat. Commun., 2022, 13, 1120.
18. R. Shu, Y. Zhou, Q. Wang, Z. Han, Y. Zhu, Y. Liu, Y. Chen, M. Gu, W. Xu, Y. Wang, W. Zhang, L. Huang and W. Liu, Adv. Funct. Mater., 2019, 29, 1807235.
19. J. Zhang and B. B. Iversen, J. Appl. Phys., 2019, 126, 085104.
20. J. Zhang, L. Song, M. Sist, K. Tolborg and B. B. Iversen, Nat. Commun., 2018, 9, 4716.
21. G. Frenking, The Chemical Bond, Wiley-VCH, Weinheim, 2014.
22. F. R. Wagner and Y. Grin, in Comprehensive Inorganic Chemistry III, Elsevier, 2023, pp. 222–237.
23. J. Wang, J. Mark, K. E. Woo, J. Voyles and K. Kovnir, Chem. Mater., 2019, 31, 8286–8300.
24. A. Li, C. Fu, X. Zhao and T. Zhu, Research, 2020, 2020, 2020/1934848.
25. H. Shang, Z. Liang, C. Xu, J. Mao, H. Gu, F. Ding and Z. Ren, Research, 2020, 2020, 2020/1219461.
26. Y. Pan, M. Yao, X. Hong, Y. Zhu, F. Fan, K. Imasato, Y. He, C. Hess, J. Fink, J. Yang, B. Büchner, C. Fu, G. J. Snyder and C. Felser, Energy Environ. Sci., 2020, 13, 1717–1724.
27. X. Sun, X. Li, J. Yang, J. Xi, R. Nelson, C. Ertural, R. Dronskowski, W. Liu, G. J. Snyder, D. J. Singh and W. Zhang, J. Comput. Chem., 2019, 40, 1693–1700.
28. R. Wang, Z. Guo, Q. Zhang, J. Cai, G. Liu, X. Tan and J. Jiang, J. Mater. Chem. A, 2022, 10, 11131–11136.
29. M. K. Brod, S. Anand and G. Jeffrey Snyder, Mater. Today Phys., 2023, 31, 100959.
30. P. Giannozzi, O. Andreussi, T. Brumme, O. Bunau, M. Buongiorno Nardelli, M. Calandra, R. Car, C. Cavazzoni, D. Ceresoli, M. Cococcioni, N. Colonna, I. Carnimeo, A. Dal Corso, S. de Gironcoli, P. Delugas, R. A. DiStasio, A. Ferretti, A. Floris, G. Fratesi, G. Fugallo, R. Gebauer, U. Gerstmann, F. Giustino, T. Gorni, J. Jia, M. Kawamura, H.-Y. Ko, A. Kokalj, E. Küçükbenli, M. Lazzeri, M. Marsili, N. Marzari, F. Mauri, N. L. Nguyen, H.-V. Nguyen, A. Otero-de-la-Roza, L. Paulatto, S. Poncé, D. Rocca, R. Sabatini, B. Santra, M. Schlipf, A. P. Seitsonen, A. Smogunov, I. Timrov, T. Thonhauser, P. Umari, N. Vast, X. Wu and S. Baroni, J. Phys.: Condens. Matter, 2017, 29, 465901.
31. J. P. Perdew, K. Burke and M. Ernzerhof, Phys. Rev. Lett., 1996, 77, 3865–3868.
32. D. R. Hamann, Phys. Rev. B: Condens. Matter Mater. Phys., 2013, 88, 085117.
33. H. J. Monkhorst and J. D. Pack, Phys. Rev. B: Condens. Matter Mater. Phys., 1976, 13, 5188–5192.
34. E. Zintl and E. Husemann, Z. Phys. Chem. B, 1933, 21, 138–155.
35. J. Zhang, L. Song and B. B. Iversen, npj Comput. Mater., 2019, 5, 76.
36. L. Bellaiche and D. Vanderbilt, Phys. Rev. B: Condens. Matter Mater. Phys., 2000, 61, 7877–7882.
37. W. Tang, E. Sanville and G. Henkelman, J. Phys.: Condens. Matter, 2009, 21, 084204.
38. M. I. Aroyo, D. Orobengoa, G. de la Flor, E. S. Tasci, J. M. Perez-Mato and H. Wondratschek, Acta Crystallogr., Sect. A: Found. Adv., 2014, 70, 126–137.
39. G. Pizzi, V. Vitale, R. Arita, S. Blügel, F. Freimuth, G. Géranton, M. Gibertini, D. Gresch, C. Johnson, T. Koretsune, J. Ibañez-Azpiroz, H. Lee, J.-M. Lihm, D. Marchand, A. Marrazzo, Y. Mokrousov, J. I. Mustafa, Y. Nohara, Y. Nomura, L. Paulatto, S. Poncé, T. Ponweiser, J. Qiao, F. Thöle, S. S. Tsirkin, M. Wierzbowska, N. Marzari, D. Vanderbilt, I. Souza, A. A. Mostofi and J. R. Yates, J. Phys.: Condens. Matter, 2020, 32, 165902.
40. N. Marzari, A. A. Mostofi, J. R. Yates, I. Souza and D. Vanderbilt, Rev. Mod. Phys., 2012, 84, 1419–1475.
41. N. Marzari and D. Vanderbilt, Phys. Rev. B: Condens. Matter Mater. Phys., 1997, 56, 12847–12865.
42. R. O. Jones, S. R. Elliott and R. Dronskowski, Adv. Mater., 2023, 2300836.
43. M. K. Brod, M. Y. Toriyama and G. J. Snyder, Chem. Mater., 2020, 32, 9771–9779.
44. S. Boyer-Richard, C. Katan, B. Traoré, R. Scholz, J.-M. Jancu and J. Even, J. Phys. Chem. Lett., 2016, 7, 3833–3840.
45. T. A. Albright, J. K. Burdett and M.-H. Whangbo, Orbital Interactions in Chemistry, Wiley, Hoboken, New Jersey, 2nd edn, 2013.
46. J. M. Ziman, Principles of the Theory of Solids, Cambridge Univ. Press, Cambridge, 2nd edn, 1999.
47. R. Hoffmann, Angew Chem. Int. Ed. Engl., 1987, 26, 846–878.
48. J. C. Slater and G. F. Koster, Phys. Rev., 1954, 94, 1498–1524.
49. C. Xia, J. Cui and Y. Chen, ACS Appl. Electron. Mater., 2020, 2, 2745–2749.
50. G. Pizzi, D. Volja, B. Kozinsky, M. Fornari and N. Marzari, Comput. Phys. Commun., 2014, 185, 422–429.
51. R. Dronskowski and P. E. Bloechl, J. Phys. Chem., 1993, 97, 8617–8624.

---

Footnote

† Electronic supplementary information (ESI) available. See DOI: https://doi.org/10.1039/d3ta04192a

---