Ternary mixed-anion semiconductors with tunable band gaps from machine-learning and crystal structure prediction

Ternary mixed-anion semiconductors with tunable band gaps from machine-learning and crystal structure prediction

Maximilian Amsler amsler.max@gmail.com Laboratory of Atomic and Solid State Physics, Cornell University, Ithaca, New York 14853, USA    Logan Ward Present address: Computation Institute, University of Chicago, Chicago, IL 60637, USA Department of Materials Science and Engineering, Northwestern University, Evanston, Illinois 60208, USA    Vinay I. Hegde Department of Materials Science and Engineering, Northwestern University, Evanston, Illinois 60208, USA    Maarten G. Goesten Department of Chemistry and Chemical Biology, Cornell University, Ithaca, New York 14853, USA    Xia Yi Center for Nanoscale Materials, Argonne National Laboratory, Argonne, Illinois 60439, USA    Chris Wolverton c-woverton@northwestern.edu Department of Materials Science and Engineering, Northwestern University, Evanston, Illinois 60208, USA
December 29, 2018

We report the computational investigation of a series of ternary \ceX4Y2Z and \ceX5Y2Z2 compounds with X={Mg, Ca, Sr, Ba}, Y={P, As, Sb, Bi}, and Z={S, Se, Te}. The compositions for these materials were predicted through a search guided by machine learning, while the structures were resolved using the minima hopping crystal structure prediction method. Based on ab initio calculations, we predict that many of these compounds are thermodynamically stable. In particular, 21 of the \ceX4Y2Z compounds crystallize in a tetragonal structure with I-42d symmetry, and exhibit band gaps in the range of 0.3 and 1.8 eV, well suited for various energy applications. We show that several candidate compounds (in particular \ceX4Y2Te and \ceX4Sb2Se) exhibit good photo absorption in the visible range, while others (e.g., \ceBa4Sb2Se) show excellent thermoelectric performance due to a high power factor and extremely low lattice thermal conductivities.

I Introduction

Computational approaches are being employed at an increasing rate to discover and design novel materials with tailored properties to tackle global environmental challenges. Traditionally, two approaches are frequently employed for these efforts: while high-throughput density functional theory (DFT) calculations (1; 2; 3; 4; 5; 6) are popular to explore the chemical space, crystal structure prediction schemes (CSP) (7) are used to determine the ground state configuration at a given composition by globally minimizing the (free) energy.

Very recently, novel methods based on materials informatics and machine learning (ML) models have emerged to assist the search for materials with improved properties in industrially relevant applications. Trained on available materials data either from experimental observations or DFT datasets (like the OQMD (1; 2), Materials Project (3), and AFLOWlib (4)), most of these ML models aim at predicting promising chemical subspaces for materials with favorable properties (8). Examples include models for melting temperatures (9), materials thermodynamics (10; 11; 12; 13; 14; 15; 16; 17), mechanical properties of alloy systems (18; 19; 20), superconductivity (21), formation of metallic glasses (22), and electronic properties of semiconductors and insulating materials (23; 24; 25; 26; 27).

Despite these promising developments, a major drawback of such ML schemes is that they often predict promising chemical compositions without providing any information on the underlying crystal structure. However, the knowledge of the atomic arrangement in a crystal lattice is crucial for any further computational assessment from first principles calculations. Hence, it is desirable to augment any ML prediction of composition alone in some way with a determination of the corresponding ground state crystal structure to (a) validate the ML model and (b) to facilitate further computational investigations in terms of materials properties.

To tackle this issue, we combine recently developed ML models by Meredig et al. (11) and Ward et al. (25) to predict the formation enthalpies and band gap energies with a sophisticated CSP scheme, the minima hopping method (MHM) (28; 29). The MHM takes as the input the chemical composition from the ML model and optimizes the potential energy, providing promising candidates for the ground state structure. We apply this approach on a subset of complex, ternary chemistries that are predicted by the ML to exhibit finite band gaps in a particularly promising range for many energy applications, including photovoltaics, photocatalysis, and thermoelectrics. We predict a class of ternary compounds with \ceX4Y2Z composition that are both thermodynamically stable and have band gaps between 0.3 and 1.8 eV. Our calculations show that several compounds are indeed promising for photovoltaic applications, most notably \ceX4Y2Te and \ceX4Sb2Se, with a strong overlap of the absorption and solar spectrum. At the same time, several candidates exhibit large power factors and low thermal conductivities due to the underlying structural complexity, rendering them promising candidates as thermoelectric materials. In fact, \ceBa4Sb2Se shows an extremely low thermal conductivity of  WmK at 300 K, close to the best thermoelectric material known to date, SnSe (30; 31).

Ii Methods

ii.1 Machine learning model

Machine learning models are composed of a set of training data, a representation that converts the training data into a form suitable for machine learning (i.e., tensors), and a machine learning algorithm. For this work, we use the training set, representation, and algorithms described in Ref. 25. Specifically, we use the formation enthalpy and band gap energy data from the Open Quantum Materials Database (OQMD) (1; 2) to train our models. From the entire OQMD, we only use the lowest energy structure at each composition - a total of 228,649 entries. Our representation is composed of 145 different attributes of the composition of each entry, such as the mean melting temperature of the constituent elements and whether the compound is charge-balanced. We use a hierarchical machine learning model based on decision trees that was found in Ref.  (25) to predict the band gap energies, and the Random Forest algorithm (32) to predict formation enthalpies. All models are created using the Material Agnostic Platform for Informatics and Exploration (Magpie) (25), and are available as Supplementary Information to this paper and on GitHub. (33)

ii.2 Structural searches

We employ the Minima Hopping Method (MHM) (28; 29) to explore candidate structures for the most promising compositions proposed by the ML model. The MHM implements a reliable algorithm to explore the low-lying portions of a potential energy surface solely given the chemical composition. Consecutive, short molecular dynamics (MD) escape trials are employed to overcome energy barriers, followed by local geometry optimizations. The Bell-Evans-Polanyi principle is exploited by aligning the initial MD velocities along soft-mode directions in order to accelerate the search (34; 35). In the past, the MHM has been successfully employed to predict or resolve the structure of a wide class of materials (36; 37; 38; 39; 40; 41; 42; 43).

ii.3 Formation enthalpies and electronic struture

The training data for the ML model is generated from DFT in a high-throughput fashion and is freely available in the OQMD. All DFT calculations are performed with the Vienna Ab initio Simulation Package (VASP) (44; 45; 46) within the projector augmented wave (PAW) formalism (47; 48) in conjunction with the PBE parameterization of the generalized gradient approximation to the exchange correlation functional (49). We use -centered -point meshes with  -points per reciprocal atom and a plane-wave cutoff energy of 520 eV. All structural relaxations are carried out by taking into account the atomic and cell degrees of freedom until the force components on the atoms were within 0.01 eV/Å, and stresses were within a few kbar. Similarly, the phases resulting from the structural searches are refined with the same DFT settings to obtain comparable formation energies on the same technical footing.

The band gaps used to train the ML model are based on the PBE functional. Although semi-local functionals are known to underestimate the band gaps, the overall correlation with respect to different materials chemistries is well reproduced, i.e., the error is systematic and corresponds to a constant shift of the band gap values.

To analyze the chemical bonding, we study the Crystal Orbital Overlap Population (COOP), which decomposes the electronic density of states into bonding and anti-bonding contributions. For this purpose, we use the the ADF Band package (50; 51; 52; 53), which employs basis sets of numerical atomic orbitals together with Slater-type orbitals. We use double zeta polarized ZORA basis functions (ZORA-DZP). This basis set is found to be sufficiently converged by carefully inspecting the band structures. Ab initio COOPs obtained from linear combination of atomic orbitals (LCAO) are especially sensitive to numerical issues arising from basis set over-completeness due to linear dependencies in construction of the overlaps.

ii.4 Optical properties

The optical absorption properties for selected candidate materials are evaluated by computing the frequency dependent dielectric function , using the Heyd-Scuseria-Ernzerhof (HSE06) hybrid functional (54; 55; 56; 57). We employ a cutoff energy of 270 eV and a -point mesh, and include at least 124 virtual bands (corresponding to more than twice the number of occupied states).

ii.5 Electronic transport

The electronic transport properties are computed by solving the Boltzmann transport equation within the constant relaxation time approximation using the Boltztrap code (58). This approximation is commonly applied for doped semiconductors and assumes that the electronic relaxation time varies little with energy on the scale of  (59). The PBE band structures are resolved on a -point mesh. Typical values for the relaxation time is around  s (60). Here, we use a constant relaxation time of  s in accordance with the work of Bilc et al. (61).

ii.6 Lattice dynamics

Phonons are computed by using the finite difference approach to obtain the second order interatomic force constants. Supercells of dimension are used, and all calculations are performed with the Phonopy package (62).

The thermal transport calculations are carried out by taking into account three-phonon interactions. We use the recently developed compressive sensing lattice dynamics (63) (CSLD) technique to obtain the third order force constants (FC). For every compound, the FC are fitted to atomic forces from 40 randomly displaced supercells containing 756 atoms, using a cutoff energy of 400 eV and a -only -points setting. Cutoff radii of 14 Å and 5 Å are used to truncate the 2- and 3-body interactions, respectively. The FC are subsequently fed into the ShengBTE package (64) to iteratively solve the linearized Boltzmann phonon transport equation.

Iii Results and Discussion

iii.1 Discovery of New Ternary Semiconductors

Our initial search for new semiconducting materials is guided by the compounds predicted to be stable by the ML model of Meredig et al. (11). One of the compounds predicted to be stable is \ceBa2As2S5. As shown in Fig. 1(a), we are able to confirm this earlier prediction using the newer ML models from Ward et al. (25). We then evaluate all Ba-As-S compositions with less than 12 atoms per formula unit to detect other potentially-stable materials in this system. We find that there exists a second region of stable compounds in the proximity of the Ba–As binary system. We then evaluate the Ba–As–S system with the band gap energy ML model from Ward et al. (25) and find that a large fraction of this ternary system is predicted to have band gaps of approximately 1 eV (see Fig. 1(b)). Both the stability and band gap energy predictions indicate that the Ba–As–S system contains promising, yet-undiscovered semiconducting materials.

Figure 1: (a) Stability and (b) band gap energy of various compositions in the Ba–As–S ternary system predicted using ML models. The stability is determined by first predicting the formation enthalpy of a composition using ML, and then measuring the hull distance obtained from the OQMD. We find two regions with negative formation enthalpies, which indicates the presence of a potential compound currently missing from the convex hull of the OQMD. The band gap energy model predicts that a significant fraction of the Ba–As–S compounds have band gap energies of approximately 1 eV. Panel (c) shows the convex hull construction based on DFT energies, including all experimentally observed phases and the putative ground states from MHM simulations for \ceBa4As2S and \ceBa5As2S2. Blue and orange circles denote thermodynamically stable and unstable phases, respectively.

Given the two composition regions of interest, we then explore the literature to further refine a list of candidate materials. We find that \ceBa2As2S5 is an already known material that was however not present in the OQMD at the time our training set was assembled. This result provides strong confidence that our ML models indeed yield reliable predictions. On the other hand, to the best of our knowledge, there exist no experimental evidence of a ternary compound near the Ba– and As–rich portion of the phase diagram. To generate a first list of candidates, we enumerate compositions that correspond to the common oxidation states of \ceBa^2+, \ceAs^3- and \ceS^2-. Starting with the smallest possible formula units (f.u.) to render CSP computationally tractable, we are left with the ternary compositions \ceBa4As2S and \ceBa5As2S2, with number of atoms per f.u. of 7 and 9, respectively.

Figure 2: The structural details of the \ceX4Y2Z compounds with I-42d symmetry, predicted through the MHM simulations. Here, \ceBa4As2S serves as a representative compound. Panels (a), (b) and (c) show the crystal structure from three different perspectives with the distorted, corner sharing tetrahedra of \ceXY4 (\ceBaAs4, blue) and \ceXZ4 (\ceBaS4, yellow) in a polyhedral representation. Panel (d) shows the two enantiomeric networks formed by the X (Ba) atoms.

Based on this initial assessment, we perform MHM calculations with two f.u. for the two systems \ceBa4As2S and \ceBa5As2S2. In addition, we also explore the combination of all elemental substitutions in \ceX4Y2Z and \ceX5Y2Z2 with X={Mg, Ca, Sr, Ba}, Y={P, As, Sb, Bi}, and Z={O, S, Se, Te}, i.e., 64 distinct compositions. Each of the MHM runs is terminated after finding some 100 distinct structures (corresponding to roughly 200 MHM iterations). For many systems of similar size, this number of iterations might not be sufficient to conclusively determine the ground state structure. Therefore, in a second step, we adopt a high-throughput substitutional search to validate our predictions. For this purpose, we select the lowest energy candidates within each of the MHM runs in the 64 distinct compositions as prototype structures. These prototypes are then decorated with the elements of the remaining 63 compositions before performing a single, local geometry relaxation with refined parameters. Among those, the lowest energy structure \ceX4Y2Z and \ceX5Y2Z2 are selected as the final candidate phases for the putative ground states.

To evaluate the thermodynamic stability of the candidates, we rely on the DFT phase stability by taking into account all phases that have been computed in a high-throughput fashion in the OQMD (1; 2). The database includes both experimentally observed phases reported in the ICSD as well as a wealth of hypothetical compounds, created by new decorations of frequently observed structural prototypes from the ICSD. The energies of all phases are subsequently used to construct the convex hulls for each of the X–Y–Z systems. The Gibbs triangle convex hull for the Ba–As–S system is shown in Fig. 1(c), and analogous figures for all other ternary systems can be found in the Supplementary Materials.

In several systems we find that either \ceX4Y2Z or \ceX5Y2Z2, and in some cases even both compositions, are predicted to be thermodynamically stable. The only exceptions are Ba–Bi–S, Ba–Bi–Se, Ba–Sb–S, Ca–Bi–S, Ca–Bi–Se, Ca–Sb–S, Sr–Bi–S, Sr–Bi–Se, Sr–Sb–S, and all magnesium containing compounds, Mg–Y–Z. Hence, phases where the Y and Z elements are two or more periods apart have a tendency to be unstable. Further, we notice a general trend that \ceX4Y2Z compounds are more stable than their \ceX5Y2Z2 counterparts, i.e., a \ceX5Y2Z2 compound only lies on the convex hull of stability if also the corresponding \ceX4Y2Z compound is thermodynamically stable. The only exception is \ceCa5Sb2Se2 which is sufficiently low in energy to push \ceCa4Sb2Se away from the hull. Therefore, we henceforth focus on the characterization of the \ceX4Y2Z phases only.

Many of the oxide compounds \ceX4Y2O have been experimentally reported, and are known since the early 1970’s. \ceCa4Y2O (65; 66; 67), \ceSr4Y2O (68; 69; 70), and \ceBa4Y2O (68; 69; 71) with Y={P, As, Sb, Bi} crystallize all in a \ceK2MgF4 structure type with symmetry (some with minor distortions). Our MHM simulations recover this structure type for most of the oxides, or slight distortions thereof, giving us confidence that the structural searches are well converged. For the sulfides, selenides and tellurides, however, we find a diverse set of ground state structures for the distinct chemical compositions (the putative ground state structures at every composition can be found in the Supplemental Materials. This difference in behavior is expected due to the very strong metal-oxide bonding (giving rise to well defined, deep thermodynamic wells around the ground states) compared to the weaker metal-sulfide, metal-selenide, and metal-telluride interactions which lead to a higher density of configurational states in sulfides, selenides and tellurides.

iii.2 Structure and Bonding

One particular structure with I-42d symmetry appears especially frequently in our search. In fact, the majority of the phases crystallize in this structure, which is shown in Fig. 2. Here, we use \ceBa4As2S as a representative prototype to discuss the detailed structural features on behalf of all ground state phases in this particular structure type. The tetragonal structure has cell dimensions of  Å  and  Å. The Ba cation occupies the Wyckoff position at , while the As and S atoms occupy the and sites at and , respectively. Each S atom has four nearest neighbors with a Ba–S bond length of 3.231 Å, which is slightly larger than in \ceBaS with the rocksalt structure (3.1935 Å). On the other hand, each As atom has four Ba atoms in close proximity, two at a distance of 3.242 Å, and two at 3.265 Å. These values compare well with the bond distance of 3.292 Å  in the experimentally reported \ceBa4As_2.67 compound with an anti-\ceTh3P4 structure type. The four Ba atoms surrounding each cation form strongly distorted, almost flat tetrahedra, which are shown as yellow and blue polyhedra in the panels (a), (b) and (c) of Fig. 2. These tetrahedra are linked by sharing corners to form a network structure.

The I-42d structure of \ceBa4As2S is closely related to the experimentally reported, cubic \ceBa4As_2.67 phase with I-43d symmetry (72). Note that Li et al. classify \ceBa4As_2.67 as a Zintl compound (72), with isolated As anions (73), but another way to interpret this compound is as a simple salt: A nominal composition assuming \ceBa^2+ \ceAs^3- would be \ceBa3As2, but the conventional cell contains 16 Ba atoms, resulting more conveniently in (\ceBa4As_2.67). Li et al. interpret the structure as Ba atoms forming two enantiomeric networks, isostructural to the ones shown in panel (d) of Fig. 2. The As atoms on the other hand occupy the total 12 combined As and S sites of the corresponding \ceBa4As2S structure in the voids of the Ba network. These sites however only have \sfrac89 occupancy (0.866 in experiment) to achieve the correct charge balance. In other words, the \ceX4Y2Z compounds can be interpreted as \ceTh3X4 in the inverse \ceTh3P4-type structure, where the average charge of -2.67 of the anions is achieved by substituting the Th-site with two and one elements of mixed oxidation states -2 and -3 (here, Y and Z), respectively. Hence, the \ceX4Y2Z compounds are mixed-anion analogs of \ceBa4As_2.67. These anion substitutions result in the tetragonal distortion away from the cubic symmetry.

Figure 3: The electronic structure of \ceBa4As2S. The band structure on the left is colored based on the band character from atomic projections, while the Crystal Orbital Overlap Population (COOP) is shown on the right. The grey line denoted with “Total” represents the total density of states (DOS).

Following the classification of Li et al., we can also interpret the \ceX4Y2Z compounds as Zintl (valence) phases with isolated anions Y and Z (according to the stoichiometry and assumed oxidation states), where the alkaline earth metals donate their valence electrons to the semimetals (here the pnictogens and chalcogens) to satisfy the octet rule. Hence, we can expect semiconducting/insulating properties throughout. Fig. 3 shows the electronic band structure of \ceBa4As2S as a representative prototype, with colors indicating the band character based on atomic projections. The band gap is direct, with the valence band maximum (VBM) and conduction band minimum (CBM) located at . The valence bands stem dominantly from As -type orbitals, with slight contributions from Ba with and -type character. The S states are buried further below the Fermi level, at around -2.5 eV. The conduction bands stem primarily from the Ba states.

Figure 4: The band gaps of the \ceX4Y2Z compounds computed with the PBE functional. The three panels represent values for X=Ca, Sr and Ba. In every panel, the x-axis corresponds to the Z elements, here {S, Se, Te}. The y-values of the center of every circle corresponds to the band gap maginitude, while the diameter represents the chemical element on the Y-site, here {P, As, Sb, Bi}. The color of the circles indicates the stability of the ground state structures at every composition (i.e., the convex hull distance), taking into account all relevant phases available in the OQMD. Finally, the presence of a black, filled dot at the center of the circles indicates that the ground state crystallizes in the tetragonal I-42d structure.

We do not expect strong (covalent) bonding in these I-42d structures, but a predominantly ionic behavior. Nevertheless, we inspect the orbital interactions via the COOP, since they will still play a role to some extent. In the right panel of Fig. 3 we show the COOP for the \ceBa4As2S compound computed with the ADF Band package (50; 51; 52; 53). These COOP are obtained by directly analyzing the contributions from the LCAO basis set, which are employed by the ADF Band package. In contrast to plane wave codes, no atomic projections are required in ADF Band, and therefore one can expect more accurate results. The bonding behavior is very similar to the one for \ceSr4Bi3 in the cubic anti-\ceTh3P4 structure type reported by Li et al. (72), however with a shifted Fermi energy: here, the bonding states are completely filled below the Fermi level, while the anti-bonding states lie above the band gap. Again, note how the valence bands stem predominantly from the Ba–As interactions, while the Ba–S leads only to significant bonding contributions at energies of  eV.

Figure 5: The absorption spectra of the \ceX4Y2Z compounds in the I-42d structure. The imaginary part of the dielectric function is shown with respect to the incident photon energy on top of the solar spectral irradiance (74)

Although the overall features of the band structure are very similar for all \ceX4Y2Z compounds (e.g., they all exhibit direct gaps at ), the magnitude of the band gaps changes significantly with the chemistry. Fig. 4 shows a graphical representation of the band gaps together with the thermodynamic stability. Each panel corresponds to an element for X, namely X={Ca, Sr, Ba}. Within each panel, the y-axis defines the band gap, while the x-axis and the size of the circles indicates the element on the Z and Y sites, respectively. Finally, the color of the circles shows the hull distance, a measure of thermodynamic stability.

In the literature we find few reported synthesis of some of these ternary Zintl compounds. In an attempt to react \ceCa5Sb3 with S, Hurng et al. (65) synthesized \ceCa4Sb_2.4S_0.4 in a defective cubic anti-\ceTh3P4 structure type, very close to the nominal stoichiometry of \ceCa4Sb_2S. However, based on our calculations, the ground state structure at the composition \ceCa4Sb_2S is not the corresponding, ternary tetragonal I-42d structure, but a monoclinic structure with P21/m symmetry, which is energetically favorable by 27 meV/atom. Further, we see from the color scheme in Fig. 4 that \ceCa4Sb_2S, even in its monoclinic structure, is thermodynamically unstable with respect to decomposition into competing phases, as indicated by the positive hull distance of 12 meV/atom. Hulliger (75) reported the synthesis of ternary europium compounds in the inverse \ceTh3P4 structure type (\ceEu4Y2Z with Y={P, As, Sb, Bi} and X={S, Se, Te}) as well as \ceCa4Bi2Te, \ceSm4Bi2Te and \ceYb4Bi2Te.

Based on these findings we can draw two conclusions: First, although the DFT formation energies at 0 K allows us to assess if a certain phase is thermodynamically stable or not, it does not uniquely provide us with an answer if it is synthesizeable, especially at finite temperatures. In particular, \ceCa4Sb2S is thermodynamically unstable at 0 K based on DFT, but nevertheless experimentally synthesizeable. Second, some phases that we determine to be unstable in the tetragonal I-42d structure might in fact be accessible experimentally, stabilized due to entropic effects and structural disorder. E.g., the energetic penalty of  meV/atom between the I-42d phase of \ceCa4Sb2S and the convex hull corresponds to a temperature of  K, an energy difference that can be readily overcome at synthesis conditions.

We also identify three general chemical trends. First, the band gap decreases as we move down in the periodic table for the X site. Similarly, the band gap decreases as we move down in the periodic table for the Y site. This behavior can be directly attributed to the decreasing electronegativity of the Y elements: the valence bands, which stem from their corresponding states, are pushed up in energy with respect to the conduction bands. Finally, there is a less prominent tendency of converging band gap values as we move down in the periodic table for the Z site. Most importantly, we see that understanding these trends from chemical substitutions on the three distinct sites allows the tuning of the gap in a range between 0.3 and 1.8 eV.

Note that we analyze here the band gaps using the semi-local PBE functional, which is known to systematically underestimate the true gap energies. However, we also neglect the relativistic effects of spin-orbit coupling (SOC), which overall lowers the gap energies, especially for the compounds containing very heavy elements like bismuth (see Supplemental Materials). Hence, the two effects somewhat cancel each other out. Further, we show PBE results for consistency, since the training data for the ML models are trained on data generated within the same DFT footing.

iii.3 Properties for Energy Applications

As semiconductors, this range of band gaps is especially of interest for energy applications, and we will specifically discuss two of them in detail, namely photovoltaics and thermoelectrics. For the former, the Shockley-Queisser limit determines the ideal gap for a single p-n junction photovoltaic cell at a value of 1.34 eV, resulting in a maximal (theoretical) efficiency of 33.7%. Further, the direct nature of the band gap (both the VBM and the CBM are located at , see Fig. 3) would allow for an efficient photo absorption, one of the major drawbacks of silicon solar cells that currently dominate the photovoltaic market (diamond silicon has an indirect band gap and hence requires phonon-assisted photo-absorption).

Figure 6: The power factors of the \ceX4Y2Z compounds in the I-42d structure as a function of temperature and electronic chemical potential . Solid, dark lines represent values computed at a temperature of 300 K, while the color gradient shows the results at higher temperatures in intervals of 100 K. The maximal temperature (lightest color) corresponds to 800 K.

To assess how well the \ceX4Y2Z compounds are suited for photovoltaic applications, we compute their absorption spectra and compare them with the solar spectrum. Semi-local functionals frequently underestimate the band gap, thereby also underestimating the absorption edge. We therefore compute the absorption using the hybrid HSE06 functional, which has been shown to improve the band gaps for many materials classes (76). In Fig. 5 we show the imaginary part of the frequency dependent dielectric function, which is computed on a -point mesh. Although this mesh might seem rather coarse, we performed careful convergence tests and determined that the overall features of the spectrum are well captured. In particular, we compare our results to results and by solving the Bethe-Salpeter equation (see Supplemental Materials). Note, however, that we did not include the effect of spin-orbit coupling, which lowers the band gaps, especially for the bismuth-containing compounds (see Supplemental Materials).

As expected, the compounds that exhibit the best absorption behavior over a wide photon energy range in the solar spectrum are the ones with overall low band gaps, mainly containing Te. The best candidates are \ceX4Y2Te and \ceX4Sb2Se, with an absorption edge in the range of 1.2-2.0 eV and a high absorption peak around 1.8-3.0 eV, well below the UV regime. Hence, these materials might be well suited for photovoltaic applications, especially as thin films due their direct band gaps and the associated high absorption efficiency.

Figure 7: The lattice thermal conductivities of the \ceX4Y2Z compounds in the I-42d structure as a function of temperature. Note that at very high temperatures, the validity of the Boltzmann transport equation breaks down as the phonon mean free paths become shorter than the smallest interatomic distances.

Zintl phases have been recently intensely studied as promising candidates for thermoelectric applications (77; 78). In general, the Carnot conversion efficiency of thermoelectric materials is governed by the figure of merit , where is the absolute temperature, is the thermopower, is the electrical conductivity, while e and are the electronic and lattice thermal conductivities, respectively (79). Optimal values for the power factor (PF) can be achieved in heavily doped semiconductors with a carrier concentration in the range of  cm (80; 81), with a narrow band gap which gives rise to large values of , while a sharp increase in the DOS around the Fermi level leads to an enhanced according to Mott’s theory (82). Due to their chemical and structural complexity, Zintl compounds exhibit inherently low lattice thermal conductivities together with favorable power factors (83; 84; 85). For example, the so-called 9-4-9 system represents Zintl compounds with \ceA9M_4+xPn9 stoichiometry (A={Yb, Eu, Ca, Sr}; M={Mn, Zn, Cd}, Pn={Sb, Bi}) and was originally discovered in the 1970s, but only recently has it drawn attention for thermoelectric applications (86). Through careful tuning of the carrier concentration in the inexpensive \ceCa9Zn_4+xSb9 compound, a high value of has been recently reported up to at 875 K (87).

In terms of the electronic properties, the band structure of \ceBa4As2S in Fig. 3 shows that many local extrema exist in a small energy range below the Fermi level at low-symmetry -points. These degenerate energy levels give rise to the sharp increase in the DOS around the Fermi energy. -doping could be readily achieved either through partially substituting the Y-site with group 13 or 14 elements, or by carefully tuning the Y/Z ratio on the anion sites towards values larger than 2. Similar arguments apply for the other \ceX4Y2Z compounds, where band degeneracies are found not only for the valence but also the conduction bands (e.g., \ceCa4Sb2Te, see Supplemental Materials). Analogously, -doping can be achieved by tuning the Y/Z ratio towards values below 2. The experimentally reported room-temperature Seebeck coefficients of 140 and 230 V/K (both p-type) for \ceEu4As2Te and \ceSm4Bi2Te in the anti-\ceTh3P4 structure type (75) indicate that similar values can be achieved for the compounds studied in the present work.

Within the relaxation-time approximation of the electronic Boltzmann equation, the electrical conductivity and the electronic thermal conductivity depend in the constant relaxation time . An accurate computation of from first principles is challenging, and depends on the scattering mechanisms with phonons, carriers, defects, and grain boundaries. Nevertheless, we can gain some (qualitative) insight without explicitly computing , but by approximating it with values commonly observed in doped semiconductors, in the range of  s. Here, we use  s throughout, a value which was initially obtained by Bilc et al. by fitting the electrical resistivity to experimental data for \ceFe2VAl_1-xM_x, M ={Si, Ge} (61). The electronic bands are resolved on a dense -point mesh, with the PBE exchange-correlation functional. Fig. 6 shows the power factors computed with the Boltztrap package as a function of doping and temperature. The values of range between 2 and 5 mWmK at 300 K and -doping, comparable to other state-of-the art thermoelectric materials, while values of up to 10 mWmK are predicted for -doping. These values further increase with temperatures that are close to operating conditions of thermoelectric generators.

Due to the above-mentioned issues in accurately estimating , it is difficult to directly determine the figure of merit . However, we can use the Wiedemann-Franz law to rewrite in terms of the Lorentz factor by factorizing into a purely electronic part, , and the electronic fraction of the total thermal conductivity,

In this way, the relaxation time cancels between the numerators and denominators, assuming that is wave vector independent. If we additionally neglect the effect of , the value of provides a convenient measure for the upper limit of . Fig. S1 in the Supplemental Materials shows as a function of doping and temperature for all compounds with the I-42d symmetry ground state structure. Overall, the results show that values of can be achieved with a careful doping in the appropriate regime.

In terms of lattice dynamics, Fig. S2 in the Supplemental Materials shows the phonon dispersion and the partial phonon DOS of \ceBa4As2S as a representative prototype, where every band is colored according to the amplitudes of the atomic eigendisplacements. There are three clearly distinct frequency regimes that are separated by pseudo gaps in the DOS. (a) The very low frequencies up to about 2 THz are dominated by Ba vibrations (on the X-sites), with some mixing of As (Y-site) contributions at 2-3 THz. (b) Between 4 and 5 THz we observe phonon bands that stem mainly from As (Y-site) vibrations. (c) The highest frequency modes at above 5 THz stem almost entirely from the S (Z-site) vibrations.

In addition to the phonon spectrum, we require higher order anharmonic contributions to the potential energy which give rise to phonon-phonon interaction to assess the lattice thermal conductivity . Using the recently developed CSLD technique, we fit the third order force constants, and use them to solve the linearized Boltzmann phonon transport equation using the ShengBTE package. The resulting lattice thermal conductivities are shown in Fig. 7 as a function of temperature.

Typical for Zintl compounds, the values of are overall very low, and less than  WmK for all compounds at 300 K. We notice two general chemical trends:

  1. decreases with the increasing mass of the X-element. This behavior can be readily attributed to the overall decrease in the vibrational frequencies of the acoustic modes, which carry the majority of the heat. The phonon dispersion and the partial phonon DOS in Fig. S2 (Supplemental Materials) shows that the low-frequency modes are dominated by the vibration of the X-elements, hence it is not surprising that it strongly affects the thermal conductivity.

  2. Given fixed X and Z elements, the value of decreases with increasing mass of the Y-element. We can attribute this correlation to the strong influence of the Y-elements on the low-frequency, acoustic phonons: the phonon band structures shows a mixing of Y and X-element vibrations. The higher the mass of the Y element, the stronger the mixing.

The lowest thermal conductivity is observed for \ceBa4Sb2Se with  WmK at 300 K, which is close to the value of the best currently known thermoelectric material, \ceSnSe ( WmK)(30; 31). Note that we only report the intrinsic, bulk thermal conductivities here. Doping and structural engineering through nano structuring can further reduce the values of . In fact, we expect that, when synthesized, the anionic sites might exhibit some degree of disorder, which could further reduce the thermal conductivity while preserving the favorable electronic properties.

Iv Conclusions

In summary, we present the properties of a class of ternary semiconducting compounds, which we discovered by effectively leveraging a machine learning model to predict the composition/band gap and a structure prediction scheme to assess their underlying ground states. Using ab initio calculations, we study in detail the properties of this class of materials with \ceX4Y2Z composition. The most prominent ground state structure that we identify has I-42d symmetry and is closely related to the binary anti-\ceTh3P4 structure type with cubic symmetry.

Our calculations show that the electronic structure of these compounds can be tuned through chemical substitution, leading to band gaps that are suited for a wide range of energy applications. The direct band gaps and the favorable absorption spectra of \ceX4Y2Te and \ceX4Sb2Se are especially promising for photovoltaic applications. Further, the high band degeneracy in the valence/conduction bands lead to high power factors, which, together with the low lattice thermal conductivities, renders these materials potentially attractive in thermoelectric generators.

V Acknowledgments

We thank Prof. R. Hoffmann for valuable expert discussions. M.A. (DFT calculations) acknowledges support from the Novartis Universität Basel Excellence Scholarship for Life Sciences and the Swiss National Science Foundation (project No. P300P2-158407, P300P2-174475). M.G.G. (COOP calculations) acknowledges support from the Rubicon Research Programme (project 019.161BT.031), which is (partly) financed by the Netherlands Organization for Scientific Research (NWO). L.W. (ML model) and C.W. acknowledge the financial assistance Award 70NANB14H012 from the US Department of Commerce, National Institute of Standards and Technology as part of the Center for Hierarchical Materials Design (CHiMaD). V.I.H (OQMD calculations) acknowledges support from the National Science Foundation Materials Research Science and Engineering Center (MRSEC) of Northwestern University (DMR-1720139). The computational resources from the Swiss National Supercomputing Center in Lugano (projects s621,s700, and s861), the Extreme Science and Engineering Discovery Environment (XSEDE) (which is supported by National Science Foundation grant number OCI-1053575), the Bridges system at the Pittsburgh Supercomputing Center (PSC) (which is supported by NSF award number ACI-1445606), the Quest high performance computing facility at Northwestern University, and the National Energy Research Scientific Computing Center (DOE: DE-AC02-05CH11231), are gratefully acknowledged.


Comments 0
Request Comment
You are adding the first comment!
How to quickly get a good reply:
  • Give credit where it’s due by listing out the positive aspects of a paper before getting into which changes should be made.
  • Be specific in your critique, and provide supporting evidence with appropriate references to substantiate general statements.
  • Your comment should inspire ideas to flow and help the author improves the paper.

The better we are at sharing our knowledge with each other, the faster we move forward.
The feedback must be of minimum 40 characters and the title a minimum of 5 characters
Add comment
Loading ...
This is a comment super asjknd jkasnjk adsnkj
The feedback must be of minumum 40 characters
The feedback must be of minumum 40 characters

You are asking your first question!
How to quickly get a good answer:
  • Keep your question short and to the point
  • Check for grammar or spelling errors.
  • Phrase it like a question
Test description