Computational Design of Thermoelectric Alloys Through Optimization of Transport and Dopability

Alloying is a common technique to optimize the functional properties of materials for thermoelectrics, photovoltaics, energy storage etc. Designing thermoelectric (TE) alloys is especially challenging because it is a multi-property optimization problem, where the properties that contribute to high TE performance are interdependent. In this work, we develop a computational framework that combines first-principles calculations with alloy and point defect modeling to identify alloy compositions that optimize the electronic, thermal, and defect properties. We apply this framework to design n-type Ba 2(1 x)Sr2xCdP2 Zintl thermoelectric alloys. Our predictions of the crystallographic properties such as lattice parameters and site disorder are validated with experiments. To optimize the conduction band electronic structure, we perform band unfolding to sketch the effective band structures of alloys and find a range of compositions that facilitate band convergence and minimize alloy scattering of electrons. We assess the n-type dopability of the alloys by extending the standard approach for computing point defect energetics in ordered structures. Through the application of this framework, we identify an optimal alloy composition range with the desired electronic and thermal transport properties, and n-type dopability. Such a computational framework can also be used to design alloys for other functional applications beyond TE.


Introduction
Alloys are disordered solid solutions of two or more materials (metals or insulators), which could be iso-or hetero-structural. Such solid solutions may span the entire composition range (complete solid solutions) or may form in limited composition ranges with miscibility gap(s). Alloying has been utilized to design materials with superior structural and functional properties since time immemorial. From alloying copper with tin to improve the mechanical strength of tools and weaponry in the Bronze age to alloying of iron to produce steel during the Industrial Revolution, alloy design and optimization has become the cornerstone of human advancement. More recently, In x Ga 1 x N alloys have enabled tunable light emitting diodes, wherein the band gap of GaN is tuned by controlling the In content 1 . Alloying is extensively used to optimize functional properties of materials for photovoltaics (e.g., CdSe x Te 1 x , 2 CsSn x Pb 1 x I 3 hybrid perovskites 3 ), energy storage (e.g., Li 6 PS 5 x Se x Br argyrodite solid electrolytes 4 ), ferroelectrics (e.g., Sc x Al 1 x N 5 ), and thermoelectrics (e.g., PbSe x Te 1 x 6 ) etc.
Unlike ordered and stoichiometric compounds, alloys represent a continuous compositional space and structural disorder adds to the complexity. Designing alloys to optimize functional properties in this continuous chemical space while also accounting for the effects of structural disorder is a formidable task. Empirical approaches have been quite successful for alloy design, especially high-throughput combinatorial thin film synthesis tech-  7 However, applications such as thermoelectrics that require synthesis of bulk (vs. thin-film) samples, alloy optimization remains a tedious experimental process relying on trial-and-error and chemical intuition.
Designing thermoelectric (TE) alloys is particularly challenging because it is a multi-property optimization problem, where the properties that contribute to high TE performance are interdependent. Alloying has been employed to optimize a variety of properties that directly or indirectly affect the TE performance: (1) enhance phonon scattering to reduce lattice thermal conductivity, 8,9 (2) promote band convergence to improve power factor, 8,[10][11][12] (3) tune band gap to suppress bipolar conduction, 8,10,11 and (4) optimize dopability 13 . Figure 1 summarizes some of the beneficial effects of alloying on the TE performance, including enhanced power factor through band convergence, optimized doping by manipulating the native defects, and reduced lattice thermal conductivity due to increased phonon scattering and lattice softening. Notably, automated high-throughput experimental procedures are being developed to synthesize bulk samples with the goal of discovering new TE materials and optimizing TE alloys. 14,15 Computations have greatly facilitated new TE materials discovery, [16][17][18][19] and their optimization through electronic doping 13,20 . At the same time, first-principles modeling of TE alloys has provided useful mechanistic insights 10,21 ; however, a systematic computational framework for designing and optimizing TE alloys has not been demonstrated yet. In concert with highthroughput bulk synthesis, computations can provide guidance to navigate the continuous compositional space, and identify the optimal compositions in the multi-property phase space. Outside of thermoelectrics, computations have been successfully em- Reduced κ L due to phonon scattering Reduced κ L due to phonon scattering ω Fig. 1 Beneficial effects of alloying on the electronic and thermal transport and doping of thermoelectric materials. Band convergence leads to higher power factor (S 2 s ), increased phonon scattering and lattice softening lowers lattice thermal conductivity (k L ), and changes in defect chemistry may promote dopability.
ployed to design high-entropy [22][23][24] and functionally-graded alloys 25,26 , suggesting that it is possible to develop a combined theory-experiment framework for rational design and optimization of TE alloys.
In this work, we develop a computational framework for designing TE alloys by combining first-principles calculations with alloy and point defect modeling. We believe that this is the first computational work within thermoelectrics that simultaneously considers the effects of alloying on the transport properties and doping. We apply this framework to identify specific compositions of n-type Ba 2(1 x) Sr 2x CdP 2 Zintl alloys that optimize the electronic and phonon transport properties, as well as ntype dopability. The computational framework builds upon wellestablished alloy modeling techniques, including special quasirandom structures (SQS) 27 and effective band structure (EBS) 28 to quantify the effect of alloying on the transport properties. To assess the dopability of alloys, we extend a standard approach for computing point defect energetics 29,30 in ordered structures, akin to our prior work on defect modeling of complex disordered structures. 31 Through the application of this computational framework, we identify Ba 2(1 x) Sr 2x CdP 2 alloy compositions that minimize alloy scattering of electrons, promote conduction band convergence, lower lattice thermal conductivity and are n-type dopable. We validate the predicted crystallographic and disorder properties by synthesizing and characterizing a series of alloy compositions.

Results and Discussion
We demonstrate the application of the computational alloy optimization framework by identifying Ba 2(1 x) Sr 2x CdP 2 Zintl phase alloys that simultaneously achieve: (1) conduction band convergence, (2) reduced alloy scattering of electrons, (3) reduced lattice thermal conductivity, and (4) suitable n-type dopability. As schematically illustrated in Figure 1, band convergence enhances the power factor while reduced alloy scattering of electrons minimizes reduction of the intrinsic carrier mobility. A reduction in the lattice thermal conductivity is expected due to increased phonon scattering and lattice softening. Finally, it has to be examined how alloying affects n-type dopability. The details of the firstprinciples calculations, alloy models, and defect calculations can be found in Section 4.

Crystal Structure and Site Selectivity
Zintl phases A 2 CdP 2 (A = Sr, Ba) crystallize in a noncentrosymmetric orthorhombic space group (Cmc2 1 , no. 36). There are 5 unique atomic positions, all of which are 4a Wyckoff sites (site symmetry m). The structures belong to a family of compounds, 9,30,32,33 which crystallize in the well-established Yb 2 CdSb 2 structure type. 34 This collection of compounds have been the subject of investigation in the past, including a detailed study 30 conducted by our team on the parent compounds (Sr 2 CdP 2 and Ba 2 CdP 2 ) of the Ba-Sr-Cd-P quaternaries that we presently discuss. These parent compounds feature a layered structure comprised of CdP 4 tetrahedra, forming [CdP 2 ] 4 layers ( Figure 2), with Ba/Sr atoms occupying sites between the [CdP 2 ] 4 layers (interlayer site 1) and within the [CdP 2 ] 4 layers (intralayer site 2). There are an equal number of sites 1 and 2 available in the structure. Since Ba 2(1 x) Sr 2x CdP 2 are isostructural alloys with no observed changes in the layering and coordination, the reader is referred to our previous work 30 for the overall structural description.
We model Ba 2(1 x) Sr 2x CdP 2 as random alloys using special quasirandom structures (SQS). The details of SQS generation and convergence tests are provided in the Methods (Section 4.2). The free energy of mixing (DG m ) of Ba 2(1 x) Sr 2x CdP 2 alloys is calculated from their DFT total energy and configurational entropy, as detailed in Section 4.3. We predict based on the computed DG m that Ba 2(1 x) Sr 2x CdP 2 forms a complete solid solution i.e., Sr 2 CdP 2 and Ba 2 CdP 2 are soluble for all x. In Figure 3, DG m is negative for all values of x suggesting the formation of a complete solid solution. More interestingly, our calculations reveal site selectivity of the disorder depending on whether the alloy composition is Ba-rich (x < 0.5) or Sr-rich (x > 0.5). Specifically, we find that for Ba-rich compositions (x < 0.5), it is energetically more favorable to fully occupy site 1 with Ba while confining Ba/Sr disorder to site 2 as compared to fully occupying site 2 and confining disorder to site 1. Disorder on both sites 1 and 2 is energetically even more unfavorable for all x < 0.5 and x > 0.5. Similarly, for Srrich compositions (x > 0.5), the full occupancy of site 2 by Sr and confining Ba/Sr disorder to site 1 is energetically favorable. This site selectivity is a direct consequence of the relative ionic size of Ba and Sr and the space available in the interlayer and intralayer  2) and Sr-rich (x = 0.8) compositions are shown in (B) and (D), respectively. Illustrative drawing of the structure of the 50-50 alloy is represented in (C). Ba and Sr occupy two different Wyckoff sites, with site 1 located between the layers (interlayer) and site 2 within the layers (intralayer). In Ba-rich compositions, Ba atoms fully occupy interlayer site 1 and disorder is observed at site 2. The intralayer site 2 is fully occupied by Sr in Sr-rich compositions with Sr/Ba disorder on site 1.
sites of the layered structure of Ba 2(1 x) Sr 2x CdP 2 . Figure 2 shows the site selectivity for representative Ba-rich (x = 0.2) and Sr-rich (x = 0.8) compositions.
The 50-50 alloy (x = 0.5) is unique in that the Ba/Sr disorder must occur on both sites 1 and 2 because they are present in equal numbers in the structure. The site selectivity of Ba and Sr cannot be simultaneously accommodated for x = 0.5. Our calculated DG m (Figure 3) suggests that it is possible to form the 50-50 alloy (DG m < 0) but energetically much less stable than compositions with x = 0.5 d and x = 0.5 + d , where d is a small positive number. As such, we predict that it will be challenging to experimentally synthesize the exact 50-50 alloy. our predictions of site selectivity, we synthesized a series of alloy compositions (Section 4.6) ranging from x = 0 (Ba 2 CdP 2 ) to x = 1 (Sr 2 CdP 2 ) and performed single-crystal XRD measurements (Section 4.7) to determine the crystallographic parameters, which are summarized in Table 1. Table 1 summarizes the Ba and Sr occupancy of sites 1 and 2, the cell parameters (a, b, c), and unit cell volume (V ) as a function of the alloy composition x. Our predicted site selectivity is confirmed by experiments. For Ba-rich alloys (x < 0.5), site 1 is fully occupied by Ba i.e., no Ba/Sr disorder while for Sr-rich alloys (x > 0.5), site 2 is fully occupied by Sr. Experimentally, it has proven difficult to obtain structures with composition representing perfect Ba 1.0 Sr 1.0 CdP 2 . The computed lattice parameters and unit cell volume are in excellent agreement with experimental measurements, as shown in Figure 4. Both computed and experimental values in Figure 4 are normalized by the values for Ba 2 CdP 2 (x = 0). The good agreement between theory and experiments provides confidence in the Ba 2(1 x) Sr 2x CdP 2 alloy models created using special quasi-random structures (Section 4).
With increasing Sr content (increasing x), we observe a nearlinear decrease in the unit cell volume consistent with Vegard's law (Figure 4d). Similarly, the cell parameter a exhibits an almost linear decrease with x, but cell parameters b and c change nonmonotonically (Figures 4b, 4c). Theoretically, we predict that the 50-50 alloy (x = 0.5) represents an inflection point in b and c as a function of x, which is consistent with experiments. The nonmonotonic changes and inflection in b and c are due to the site selectivity of Ba and Sr and the related structural distortions, as discussed above.
The addition of small amounts of Ba, starting from the parent compound Sr 2 CdP 2 (x = 1), causes a sharp increase in b, which is also the layer stacking direction ( Figure 2). This is explained by the fact that under Ba-poor conditions, Ba occupies the interlayer positions (site 1), a site preference which has also been observed in the Sb analogue of these quaternary structures 32 . This introduction of the larger Ba cation into the interlayer site (site 1) forces the [CdP 2 ] 4 polyanionic layers further apart from one another to accommodate the larger cation. The direct evidence for this expansion can be found in the increase in the Ba1/Sr1-P1/P2 bond lengths (Table S3), which represent the base of the square pyramidal coordination.
The sharp changes in c can also be rationalized via the larger ionic radius of Ba when compared to Sr. Under Ba-rich conditions (x < 0.5) Ba enters site 2, causing local distortions of the CdP 4 tetrahedra, including and increase in the the P1-Cd-P1 the Cd-P1-Cd bond angles (Table S4). Such distortions of the layers are indicative of the structure accommodating a larger cation in site 2. These local distortions then lead to an overall increase in the global c unit cell parameter.

Optimization of Electronic Transport
In the context of alloy optimization, there are two aspects of the electronic structure that affects carrier transport, namely, band convergence and alloy scattering. For thermoelectrics, convergence of light bands (low effective mass) is desirable to enhance the Seebeck coefficient and the power factor ( Figure 1). In contrast, alloy scattering of electronic carriers (electrons, holes) is undesirable because it reduces carrier mobility and therefore, electronic conductivity. Our previous work 30 has shown that Sr substitution on the interlayer site 1 ( Figure 2) in Ba 2 CdP 2 induces convergence of two conduction bands located at G and Z ( Figure  5). Our goal in this work is to computationally identify the optimal Ba 2(1 x) Sr 2x CdP 2 alloy for n-type transport. As such, we focus on conduction band (CB) convergence and alloy scattering of electrons.
However, random substitutional alloys lack formal translational symmetry and therefore, cannot be described by the standard bandstructure E(k) dispersion. We use the effective band structure (EBS), proposed by Popescu and Zunger, 28 to interpret the conduction band electronic structure of the alloys. The computational details of EBS calculations are provided in the Methods (Section 4.4). The EBS is obtained by unfolding the SQS supercell bands to the primitive cell Brillouin zone using spectral decomposition. The resulting EBS reveals the extent to which the band characteristics are preserved or lost. Band convergence can be determined by examining the EBS. In addition, the linewidth or band broadening at each E(k) is reflective of the electronic carrier lifetimes due to alloy scattering; the larger the band broadening (more "smeared" bands), the shorter the carrier lifetimes.
The computed EBS of Ba 2(1 x) Sr 2x CdP 2 alloys are shown for 12 different compositions in Figure 5, ranging from Ba-rich (x = 0.1) to Sr-rich (x = 0.9). We find that the relative position (on the energy axis) of the two low-lying CBs at G and Z changes with x. In the parent compounds, Ba 2 CdP 2 and Sr 2 CdP 2 , the energy separation (DE G Z ) between the two CBs are 0.26 eV and 0.35 eV, respectively ( Figure 6a). The small DE G Z in the parent compounds provides an opportunity to achieve band convergence through al- loying. In the composition range x = 0.0 0.30, the CB dispersion near G and Z closely resemble that of the parent compound, Ba 2 CdP 2 . DE G Z decreases with increasing Sr content, as seen in the regime R1 in Figure 6(a). In the regime R2, which spans x = 0.35 0.80, we observe convergence of the CBs at G and Z within 30 meV. Interestingly, for x = 0.35 and x = 0.40, the conduction band minima is located at Z, resulting in indirect gap unlike the parents compounds and the other alloy compositions. The 50-50 alloy (x = 0.5) is unique in that the two converged CBs have similar dispersion ( Figure 5) making it a candidate for optimal composition. However, we know that the 50-50 alloy is expected to be thermodynamically less stable based on the free energy of mixing (Section 2.1). Finally, for compositions with x = 0.8 1.0 in regime R3 (Figure 6a), DE G Z again increases until it reaches 0.35 eV for the parent compound, Sr 2 CdP 2 . In summary, CB convergence is achieved for compositions x = 0.35 0.70 (regime B in Figure 6). We demonstrate with this example how band convergence can be determined through the calculation of EBS.
Next, we examine band broadening in the calculated EBS to assess the extent of alloy scattering of electronic carriers. We observe significant broadening (smearing) in the valence band beyond ⇠0.5 eV from the band edge (shown only for x = 0.1 in Figure 5). Since we are interested in n-type transport, we focus on the broadening in the CB, particularly around the band edges at G and Z. For compositions between x = 0.35 0.65 (EBS for alloy composition x = 0.65 is shown in Figure S2), we observe weak band broadening, which suggests that the electronic carrier mobility in this composition range will remain largely unaffected due to alloy (disorder) scattering. As such, alloy compositions x = 0.35 0.70, which preserve the Bloch character of the CBs around G and Z, are optimal for minimizing alloy scattering of carriers. In particular, for x = 0.35 and x = 0.45, the CBs near G and Z are almost unperturbed by the alloy disorder, for reasons that are not immediately clear to us. Fortuitously, the compositions where both band convergence is achieved and alloy scattering is minimized overlap.

Optimization of Thermal Properties
Alloying generally reduces lattice thermal conductivity (k L ) due to increased phonon scattering. Further reduction in k L may result from lattice softening, as observed in SnTe-AgSbTe 2 alloys 35 and Na-doped PbTe. 36 k L can be directly determined from first principles by calculating the second and third-order interatomic force (IFCs) constants 37 . In these calculations, the effects of both phonon scattering and lattice softening are captured. However, in the finite-displacement method, the evaluation of the third-order IFCs requires pairwise atomic displacements to be performed in a chosen supercell, which rapidly scales with the supercell size. For alloys, where the "unit cell" is already a supercell, such calculations are prohibitively expensive. In contrast, the effect of alloying on lattice softening can be assessed from first principles by computing the speed of sound. In this work, we focus on lattice softening. Alloy scattering of phonons will only further reduce k L .
We computed the speed of sound (v s ) of each alloy composition from the bulk modulus (B) using the approximate relation v s = p B/r, where r is the density. 38 Figure 6b shows the variation of v s and B, both of which linearly increases with the alloy composition x following Vegard's law. Replacing Ba in the parent compound, Ba 2 CdP 2 , with Sr, which is less electropositive and has smaller atomic mass, makes the alloy stiffer (higher B). However, due to the small change in B (< 2 GPa) and v s (< 0.15 km/s) between x = 0.35 0.70, we do not expect a significant variation in k L in this range due to lattice softening/stiffening. As a result, after considering electronic transport and k L , the optimal composition range is still x = 0.35 0.70. Therefore, we choose one of the compositions within this range, Ba 0.75 Sr 1.25 CdP 2 (x = 0.625) to check if this optimal compositions is n-type dopable.

Assessment of Alloy Dopability
The thermoelectric performance of a material critically depends on whether it can be doped with the desired carrier type (electrons or holes). Furthermore, the carrier concentration of the desired carrier type needs to be tuned to optimize the TE performance. 39 The presence of high concentrations of "killer" chargecompensating native defects can prevent the desired doping. For example, high concentrations of acceptor cation vacancies in Mg 3 Sb 2 and BiCuOSe under cation-poor growth conditions, respectively, prevent n-type doping. 13,40 In the context of alloy design, it is crucial that the compositions that optimize the electronic and thermal transport are also dopable.
A material is considered n-type dopable if the most favorable native acceptor defect has high formation energy. In such a case, electrons generated by a sufficiently soluble donor dopant are not compensated by the holes created by the native acceptor. In other words, a "killer" acceptor defect is absent under the chosen condi-tions. Accordingly, an n-type dopability window (DE n ) may be defined at the CBM (Figure 7), where a large positive DE n indicates a highly n-type dopable material while a negative window would suggest difficulty in n-type doping. Similarly, a p-type dopability window (DE p ) is set by the formation energy of the favorable native donor defect at the VBM and provides a measure of the potential for p-type doping.
In this study, we are concerned with the n-type dopability of Ba 2(1 x) Sr 2x CdP 2 alloys. We chose to verify the n-type dopability of Ba 0.75 Sr 1.25 CdP 2 (x = 0.625), which lies in the optimal composition range of x = 0.35 0.70 based on the electronic and thermal transport properties (Sections 2.2, 2.3). The periodic supercell approach utilized in first-principles calculations of defect energetics is well established for ordered and stoichiometric materials. 13,29,41 A direct extension of the supercell approach to alloys is not straightforward because of the presence of different local environments in the disordered structures. As such, DFT relaxations of hundreds of supercells may be required to compute the defect formation energetics in alloys, which makes this approach prohibitively expensive. For example, in a 100-atom SQS supercell of Ba 0.75 Sr 1.25 CdP 2 , there are 25 unique Sr sites. Therefore, calculating the formation energy of Sr vacancies (V Sr ) in a single charge state (typically, multiple charge states of each defect are calculated) amounts to 25 different DFT supercell relaxations. It quickly becomes apparent that considering all types of defects, including vacancies, anti-sites, and interstitials, in different charge states requires hundreds of DFT supercell calculations. Furthermore, in a disordered material, the different local environments can lead to variations in the formation energies of the same defect. Therefore, it is fundamentally important to determine the statistical variation in the defect formation energies.
We adapt the standard supercell approach to calculate the range of formation energies of each defect type. First, we identify for each defect type, the specific sites in the chosen supercell (see Section 4.5 for details) that result in the highest and lowest formation energy in the neutral charge state (Figure 7a). We assume that, when different charge states of the defect are considered, these two specific sites will allow us to sample the range of formation energies. Next, for each defect type, we calculate the defect formation energy in multiple charge states but only at the two sites identified in the first step. Ultimately, we are able to calculate the range of formation energies as shown in Figure  7b and 7c and assess dopability of alloys. Overall, this adapted approach significantly reduces the computational overhead while providing a reliable assessment of alloy dopability. Figure 7a shows the range of formation energy for each defect type in the neutral charge state. The atomic site indices corresponding to the highest and lowest formation energy are labeled in Figure 7a. These sites are labeled in the supercell structure shown in Figure S4. The range of defect formation energies when all plausible charge states are considered are shown in Figures 7b  and 7c   discrete alloy compositions that are marked in Figure 3. The overall conclusions about the dopability will remain unchanged even if a continuous composition space is considered in the convex hull analysis. In the context of n-type dopability, cation vacancies (V Ba , V Sr , V Cd ) are found to be the "killer" electron-compensating acceptor defects that may limit n-type doping. The formation of cation vacancies will be suppressed under cation-rich growth conditions. As such, we identified the four-phase stability region that corresponds to the most cation-rich conditions such that the ntype dopability window DE n is maximized. Figure 7b shows the range of formation energy (shaded regions) of each native defect under the most favorable chemical potential conditions that maximizes DE n (P-1 in Table S9 (Figure 7c). This is atypical for Zintl phases that are commonly p-type seminconductors. 42 Even under the most favorable chemical potential conditions (P-16 in Table S9, where Ba 0.75 Sr 1.25 CdP 2 is in equilibrium with SrCd 2 P 2 , SrP 3 , Ba 0.9 Sr 1.1 CdP 2 ), the p-type dopability window DE p is large and negative, which means that holes generated by extrinsic dopants will be compensated by the electrons from the native donors. Nevertheless, our focus in this work is to optimize the n-type transport of Ba 2(1 x) Sr 2x CdP 2 alloys. We have previously studied 30  The decrease in the n-type dopability window could be attributed to the increase in the band gap with increase in Sr -Ba 2 CdP 2 has a smaller band gap (1.43 eV) 30 compared to Ba 0.75 Sr 1.25 CdP 2 (1.81 eV). If one were to assume that a fraction of the band gap increase results in a shift of the conduction band minimum to higher energy, then a simple extrapolation of the V Cd formation energy (line with negative slope -2 in Figure  7b) could explain the decrease in DE n . Based on this observation, we formulate that the alloy compositions with lower Sr content (smaller x) are optimal for maximizing n-type dopability, as quantified by DE n . Taking the electronic and thermal transport, and dopability factors into consideration, we computationally predict that Ba 2(1 x) Sr 2x CdP 2 alloys with x = 0.35 0.70 are optimal, with lower x in that range preferred for larger dopability (DE n ).

1-11 | 7
In summary, we have developed a computational framework for systematic optimization of the transport properties and dopability of thermoelectric alloys. We simultaneously consider the effects of alloying on electronic and thermal transport, and dopability, which makes this computational framework unique compared to other computational studies of thermoelectric alloys. We apply this framework to identify the optimal compositions of n-type Ba 2(1 x) Sr 2x CdP 2 Zintl thermoelectric alloys by achieving conduction band convergence, minimization of alloy scattering of electrons, lowering of k L , and ensuring n-type dopability. Our computational framework can be used to accelerate the design of thermoelectric alloys with desired transport and doping properties. Such a framework can be used to design alloys for other functional applications beyond thermoelectrics.
We understand that the computational framework presented in this work represents a step towards creating a more comprehensive framework that also takes into account the effects of alloying on phonon scattering, for example, by considering effective phonon band structures 43 and by developing computationally cheaper models to directly estimate the carrier mobility in alloys, akin to the semi-empirical models developed by one of the co-authors. 38 Looking forward, data-driven methods will play a key role in computational alloy design. Among the promising developments are the application of machine learning (ML) to design of high-entropy alloys, [22][23][24] and the development of new ML algorithms that optimize fractional compositions. 44 In concert with computational predictions, high-throughput bulk and thin-film synthesis will play a role in validation of the computational predictions. We may envision the synergy of first-principles alloy modeling (demonstrated in this work), data-driven methods, and high-throughput bulk synthesis to tackle the formidable challenge of alloy design and optimization.

Acknowledgements
JQ is funded by the NSF DIGI-MAT program, grant no. 1922758. AB acknowledges support from the University of Delaware's undergraduate research program. The experimental part of this research, led by SB and conducted at the University of Delaware, received financial support from the US Department of Energy through a grant DE-SC0008885. P.G. acknowledges support from NSF through award DMR-2102409. The research was performed using computational resources sponsored by the Department of Energy's Office of Energy Efficiency and Renewable Energy and located at the NREL. The views expressed in the article do not necessarily represent the views of the DOE or the U.S. Government.

Structure Relaxation and Bulk Modulus Calculation
First-principles calculations, including structural relaxations, were performed with density functional theory (DFT) using the Vienna Ab Initio Simulation Package (VASP) software package. 45 The generalized gradient approximation (GGA) of Perdew-Burke-Ernzerhof (PBE) was used as the exchange correlation func-tional. 46 The core and valence electrons were treated with the projector-augmented wave (PAW) method. 47 A plane-wave energy cutoff of 340 eV was used, and an automatically generated G-centered regular k-point mesh was used to sample the Brillouin zone. Bulk modulus (B) was determined by fitting the Birch-Murnaghan equation of state to a set of total energies computed at different volumes that were expanded and contracted around the equilibrium volume). 48

SQS Structures and Convergence Tests
Alloy supercell models were constructed using special quasirandom structures (SQS). 27 The concept of SQS was developed to mimic random alloys without using very large supercells or using many configurations. The design of SQS involves a stochastic search over many possible configurations of local environments within a given supercell to best reproduce the pairwise correlation of a completely random alloy. An optimal SQS is the one that minimizes the root-mean-square deviation from the random pairwise correlation.
The SQS supercells for each alloy composition were constructed using the Alloy Theoretic Automated Toolkit (ATAT). 49 In addition to minimizing the deviation from the random pairwise correlation, we also considered different supercell sizes to check for convergence of electronic structure features (band gap, effective masses). The results of the convergence tests for Ba 0.4 Sr 1.6 CdP 2 are summarized in Table S7. We found that a 100atom SQS supercell with a cutoff distance of 6 Å for pairs and triplets provides sufficiently converged electronic structure properties. We used a 100-atom SQS for all compositions (x) to calculate their free energy of mixing ( Figure 3) and effective band structures ( Figure 5).
The cations (Ba, Sr) can occupy two distinct Wyckoff positions in the Ba 2(1 x) Sr 2x CdP 2 structure (sites 1 and 2 in Figure 2). To model the site selectivity of the cation disorder, for each composition, we considered different SQS supercells with the disorder constrained to a specific Wyckoff site. Given the equal stoichiometric ratio of sites 1 and 2, for the 50-50 alloy, only one SQS supercell with disorder on both sites was considered. Each SQS supercell was fully relaxed (volume, cell shape, and atomic positions) with DFT using the functionals and parameters described in Section 4.1.

Calculation of Free Energy of Mixing
The free energy of mixing DG mix in Figure 3 was calculated as DG mix = DH mix T DS mix , where DH mix and DS mix are the enthalpy and entropy of mixing, respectively. DH mix is calculated as, where w i is the fraction of each parent compound in the alloy. H alloy is calculated from the DFT total energy of the SQS supercell and H i from the total energy of the unit cell of the parent compounds. We assume that the configurational entropy due to disorder in the alloy is the primary contribution to DS mix . The configurational entropy is calculated using the Boltzmann's for-mula as, DS config = Â j y j ln y j (2) where j is the chemical species with fractional occupation, and y j is the fraction of the disordered sites occupied by j. The fractions y j should add up to 1 i.e., Â j y j = 1. The vibrational entropy change is assumed to be zero due to the approximate cancellation of the vibrational entropy between the alloy and the parent compounds. A typical synthesis temperature of 823 K is used to compute DG mix in Figure 3. The trends remain unchanged at higher synthesis temperatures.

Effective Band Structures
We use the effective band structure (EBS) 28 to interpret the conduction band electronic structure of the alloys. The EBS is obtained by unfolding the SQS supercell bands to the primitive cell Brillouin zone using spectral decomposition. Each eigenstate in the primitive cell is assigned a spectral weight that represents how well the Bloch character of that eigenstate has been preserved when unfolded from the supercell wavefunctions. The spectral weight is defined by: where |y SC where d is the Delta function. Finally, the EBS is obtained by integrating the spectral function (A(k, E)) in an infinitesimal energy window dE: where d N(k i , E j ) gives the number of SC bands crossing PC wave vector k i and energy level E j , which is shown as the scale bar in Figure 5. EBS is calculated along the special k-point paths in the Brillouin zone of Ba 2 CdP 2 . We use 40 interpolation points between each pair of high-symmetry k points of the PC. We also calculate the EBS with 20 interpolation points (see comparison in Figure S3) to confirm that there are no interpolation artifacts. The EBSs are calculated using a modified version of the BandUP code. 50

Defect Energetics
We adapted the standard supercell approach 29 to calculate the defect energetics. First, we build a 40-atom SQS structure as the "unit cell" for the Ba 0.75 Sr 1.25 CdP 2 (x = 0.625) alloy. Next, we create a 2⇥2⇥1 supercell of the 40-atom SQS such that the supercell contains 160 atoms. We perform convergence test with 120-atom SQS to show that the 40-atom SQS reproduces the EBS and related electronic properties of the larger SQS supercell (Table S8). The defect formation energy (DE D,q ) is calculated from DFT total energies as follows: where DE D,q is the formation energy of a defect D in charge state q, E D,q is the total energy of the supercell containing defect D in charge state q, and E host is the total energy of the defect-free, neutral host supercell. E F is the Fermi energy, which referenced to the valence band maximum (VBM) i.e. E F = 0 at the VBM. E corr term contains the finite-size corrections, which are calculated following the methodology of Lany and Zunger 29 . The finite-size corrections comprises: (1) image charge correction for charged defects, (2) potential alignment correction for charged defects, and (3) band-filling correction for shallow defects. The pyladadefects software 51 is used to generate the defect supercells and to calculate the finite-size corrections.
The chemical potential of elemental species i is denoted by µ i and n i is the number of atoms of species i added (n i < 0) or removed (n i > 0) from the host supercell to form the defect D. µ i is expressed as µ i = µ 0 + Dµ, where µ 0 is the reference elemental chemical potential and Dµ i is the deviation from the reference value. µ 0 are fitted to a set of experimentally-measured formation enthalpy of compounds in the quaternary Ba-Sr-Cd-P chemical space. The fitted µ 0 values are listed in Table S10. For a material, the accessible range of Dµ i is constrained by the condition of its thermodynamic phase stability and is calculated through a convex hull analysis. Table S9 lists the Dµ i values at the corners of the polyhedra that defines the region of phase stability of Ba 0.75 Sr 1.25 CdP 2 . The positions of the ions in the defect supercells are relaxed with DFT using a plane-wave energy cutoff of 340 eV. The Brillouin zone of the supercell is sampled with a G-centered 2⇥2⇥2 k-point grid. We consider native vacancies (V Ba , V Sr , V Cd , V P ) and anti-site defects (Cd P , P Cd ). The interstitial defects Ba i and Sr i are not considered in this work because they are expected to have high formation energy due to their large ionic radii. For each defect type, we consider all unique Wyckoff positions in the structure and multiple charge states q = -3, -2, -1, 0, 1, 2, 3. The underestimation of the band gap in DFT is rectified by applying individual valence and conduction band edge shifts (relative to the DFT-computed band edges) as determined from GW quasiparticle energies. 29,52 Because of the direct band gap at G, the GW quasi-particle energies are calculated for the 40-atom SQS using a G-point only k-point grid.

Material Synthesis and Characterization
All manipulations and procedures were performed either in an inert argon environment or under a vacuum. The elements used in the synthesis of the title compounds (Ba, Cd, P, Pb, and Sr) were sourced from either Sigma-Aldrich or Alfa Aesar and carry a stated purity of, at least, 99 wt%. Due to the high purity of the elements, they were used as received. The elements were loaded, in desired stoichiometric ratios, into alumina crucibles, 1-11 | 9 which were subsequently flame sealed in silica jackets. The vessels were then placed into a programmable muffle furnace for heat treatment. The compounds synthesized, as a result of this inquiry into the Ba-Sr-Cd-P phase diagram, were made using a standard Pb flux technique wherein the elements were loaded with an excess of Pb intended as a solvent for the high temperature reaction. The heat treatment closely followed the optimized synthetic procedures used to obtain the parent Yb 2 CdSb 2 structure:1 the elements were rapidly heated to a peak temperature of 1233 K, allowed to equilibrate at this temperature for 24 hrs, and then slowly cooled to 773 K at a rate of 5 K/hr. Once the temperature of 773 K was reached, the vessels were removed from the furnace and centrifuged at high speeds to remove the excess Pb flux. Due to the experience with air sensitivity of compounds, the vessel was then transferred into an argon-filled glovebox for further manipulations. The crystals presented in this work were extracted from the matrix as black, needle-like crystals that and were used for further structural studies.

Single Crystal X-ray Diffraction
Single crystal X-ray diffraction data were collected for all of the compounds using the following procedure: Crystals, of suitable quality, were identified, while immersed in Paratone-N oil, under a high-powered microscope and then subsequently cut to size. The freshly cut crystals were then scooped up by a MiTeGen plastic loop and placed into the goniometer of a Bruker Apex-II diffractometer, which utilizes a sealed source Mo Ka radiation (0.71073 Å ) for diffraction. A constant temperature of 200 K was maintained by cold nitrogen vapor that was obtained from boiling liquid N 2 . This cold nitrogen stream also served to preserve an inert environment by hardening the Paratone-N oil, forming a hard shell around the crystal such that the decomposition of the material was prevented.
Data were collected, subsequently integrated, and corrected for absorption using the Bruker-supplied software 53 .