Phase Stability and Defect Behavior in Complex Thermoelectric Zinc-Antimonides Citation Pomrehn, Gregory Schoelerman (2013) Phase Stability and Defect Behavior in Complex Thermoelectric Zinc-Antimonides. Dissertation (Ph.D.), California Institute of Technology. doi:10.7907/0AWH-4P69. https://resolver.caltech.edu/CaltechTHESIS:01252013-135311552 Abstract The Zn-Sb binary phase system has been of interest for many years in the search for efficient and low-cost thermoelectric materials. Of primary interest has been the Zn 4 Sb 3 phase which exhibits a thermoelectric figure of merit, zT , in excess of 1 in an intermediate temperature range. In this study, Zn 4 Sb 3 is shown to be entropically stabilized with respect to decomposition to Zn and ZnSb through the effects of configurational disorder and phonon free energy. Single-phase stability is predicted for a range of compositions and temperatures. Retrograde solubility of Zn is predicted on the two-phase boundary region between Zn 4 Sb 3 and Zn. The complex temperature-dependent solubility can be used to explain the variety of nanoparticle formation observed in the system: formation of ZnSb on the Sb-rich side, Zn on the far Zn-rich side, and nano-void formation due to Zn precipitates being reabsorbed at lower temperatures. A new binary compound, Zn 8 Sb 7 , known only in nanoparticulate form, is also studied using density functional calculations. The free energies of formation, including effects from vibrations and configurational disorder, are calculated to compare with the relevant phases ZnSb, Zn, and Zn 4 Sb 3 , yielding insight into the phase stability of Zn 8 Sb 7 . Band structure calculations predict Zn 8 Sb 7 , much like ZnSb and Zn 4 Sb 3 , to be an intermetallic semiconductor with similar thermoelectric properties. If sufficient entropy or surface energy exists to stabilize the bulk material, it would be stable in a limited temperature window at high temperature. In the AZn 2 Sb 2 series of materials—A = Ca, Sr, Yb, and Eu—I show that a large concentration of thermodynamically stable cation vacancies leads to high extrinsic carrier concentrations. The stable defect level depends on the choice of A, and is consistent with experimentally observed carrier concentrations in these materials. These results demonstrate that point defects are the primary mechanism by which the covalency of the cation bond can influence carrier concentration in nominally valence-precise AZn 2 Sb 2 compounds. This mechanism may be generally applicable to other Zintl phases, perhaps explaining similar trends seen in A 14 MSb 11 , A 2 MSb 2 (A=2+ cation, M = 2+ or 3+ metal),and similar materials. Item Type: Thesis (Dissertation (Ph.D.)) Subject Keywords: Thermoelectrics; density functional theory; thermodynamics; computational materials science; cluster expansion Degree Grantor: California Institute of Technology Division: Engineering and Applied Science Major Option: Materials Science Thesis Availability: Public (worldwide access) Research Advisor(s): van de Walle, Axel (co-advisor) Snyder, G. Jeffrey (co-advisor) Thesis Committee: Fultz, Brent T. (chair) van de Walle, Axel Snyder, G. Jeffrey Johnson, William Lewis Goddard, William A., III Defense Date: 7 January 2013 Record Number: CaltechTHESIS:01252013-135311552 Persistent URL: https://resolver.caltech.edu/CaltechTHESIS:01252013-135311552 DOI: 10.7907/0AWH-4P69 Related URLs: URL URL Type Description http://dx.doi.org/10.1103/PhysRevB.83.094106 DOI UNSPECIFIED http://dx.doi.org/10.1021/ja202458n DOI UNSPECIFIED Default Usage Policy: No commercial reproduction, distribution, display or performance rights in this work are provided. ID Code: 7439 Collection: CaltechTHESIS Deposited By: Gregory Pomrehn Deposited On: 14 Feb 2013 22:03 Last Modified: 03 Oct 2019 23:58 Thesis Files Preview PDF (Complete Thesis) - Final Version See Usage Policy. 14MB Repository Staff Only: item control page

Phase Stability and Defect Behavior in Complex Thermoelectric Zinc-Antimonides Thesis by Gregory Pomrehn In Partial Fulfillment of the Requirements for the Degree of Doctor of Philosophy California Institute of Technology Pasadena, California 2013 (Defended January 7, 2013) ii © 2013 Gregory Pomrehn All Rights Reserved iii This thesis is dedicated to my loving wife, Charissa, and delight-filled sons, Noah and Nico—all of whom have kept me connected to the ‘real world’ in the midst of graduate student life and have been a constant reminder of what is truly valuable in this world. iv Acknowledgments None of this work would have been possible without the help of: Advisers: Axel van de Walle and Jeff Snyder Collaborators: Eric Toberer and Alex Zevelkink Colleagues: Ljubomir Miljacic, Pratyush Tiwary, Steve Demers, Qijun Hong, and Chirranjeevi Balaji Gopal Computers: ThinkBig, Eniac, Lonestar, Ranger, Trestles, and AppleII Church: Vision Christian Fellowship Psalm 16:5–6,11 The Lord is my chosen portion and my cup; you hold my lot. The [phase boundary] lines have fallen for me in pleasant places; indeed, I have a beautiful inheritance. You make known to me the path of life; in your presence there is fullness of joy; at your right hand are pleasures forevermore. v Abstract The Zn-Sb binary phase system has been of interest for many years in the search for efficient and low-cost thermoelectric materials. Of primary interest has been the Zn4 Sb3 phase which exhibits a thermoelectric figure of merit, zT, in excess of 1 in an intermediate temperature range. In this study, Zn4 Sb3 is shown to be entropically stabilized with respect to decomposition to Zn and ZnSb through the effects of configurational disorder and phonon free energy. Single-phase stability is predicted for a range of compositions and temperatures. Retrograde solubility of Zn is predicted on the two-phase boundary region between Zn4 Sb3 and Zn. The complex temperaturedependent solubility can be used to explain the variety of nanoparticle formation observed in the system: formation of ZnSb on the Sb-rich side, Zn on the far Znrich side, and nano-void formation due to Zn precipitates being reabsorbed at lower temperatures. A new binary compound, Zn8 Sb7 , known only in nanoparticulate form, is also studied using density functional calculations. The free energies of formation, including effects from vibrations and configurational disorder, are calculated to compare with the relevant phases ZnSb, Zn, and Zn4 Sb3 , yielding insight into the phase stability of Zn8 Sb7 . Band structure calculations predict Zn8 Sb7 , much like ZnSb and Zn4 Sb3 , to be an intermetallic semiconductor with similar thermoelectric properties. If sufficient entropy or surface energy exists to stabilize the bulk material, it would be stable in a limited temperature window at high temperature. In the AZn2 Sb2 series of materials—A = Ca, Sr, Yb, and Eu—I show that a large concentration of thermodynamically stable cation vacancies leads to high extrinsic carrier concentrations. The stable defect level depends on the choice of A, and is vi consistent with experimentally observed carrier concentrations in these materials. These results demonstrate that point defects are the primary mechanism by which the covalency of the cation bond can influence carrier concentration in nominally valence-precise AZn2 Sb2 compounds. This mechanism may be generally applicable to other Zintl phases, perhaps explaining similar trends seen in A14 MSb11 , A2 MSb2 (A=2+ cation, M = 2+ or 3+ metal), and similar materials. vii Contents Acknowledgments iv Abstract v List of Figures ix List of Tables xiv 1 Introduction 1 2 Retrograde Solubility and Entropic Stabilization in Zn4 Sb3 6 2.1 Introduction . . . . . . . . . . . . . . . . . . . . . . . . . . . . . . . . 6 2.2 Computational Methods . . . . . . . . . . . . . . . . . . . . . . . . . 9 2.3 Results and Discussion . . . . . . . . . . . . . . . . . . . . . . . . . . 11 2.4 Conclusions . . . . . . . . . . . . . . . . . . . . . . . . . . . . . . . . 19 3 Predicted Thermodynamic and Electronic Properties of a Newly Discovered Binary Zn8 Sb7 Compound 20 3.1 Introduction . . . . . . . . . . . . . . . . . . . . . . . . . . . . . . . . 20 3.2 Computational Methods . . . . . . . . . . . . . . . . . . . . . . . . . 22 3.3 Results and Discussion . . . . . . . . . . . . . . . . . . . . . . . . . . 25 3.4 Conclusions . . . . . . . . . . . . . . . . . . . . . . . . . . . . . . . . 38 4 Defect-Driven Properties in AZn2 Sb2 39 4.1 Introduction . . . . . . . . . . . . . . . . . . . . . . . . . . . . . . . . 39 4.2 Thermodynamics . . . . . . . . . . . . . . . . . . . . . . . . . . . . . 41 viii 4.3 Results and Discussion . . . . . . . . . . . . . . . . . . . . . . . . . . 43 4.4 Conclusion 48 . . . . . . . . . . . . . . . . . . . . . . . . . . . . . . . . A Cluster Expansion of III–V Alloys 50 A.1 Scalar Cluster Expansion for Energy . . . . . . . . . . . . . . . . . . 50 A.2 Tensorial Cluster Expansion of Strain and Elastic Constants . . . . . 56 Bibliography 61 ix List of Figures 1.1 Summary of significant results for Zn4 Sb3 , Zn8 Sb7 , and AZn2 Sb2 phases. (Left) Single-phase stability region for Zn4 Sb3 , with a predicted region of retrograde Zn solubility [1]. (Center) Free energy curves for Zn4 Sb3 and Zn8 Sb7 with respect to Zn and ZnSb phases [2]. Zn8 Sb7 is only predicted to be stable in a small high-temperature window. (Right) Stable cation vacancy concentration in AZn2 Sb2 for A = Ca, Sr, Yb, and Eu. Boundary lines correspond to the equilibrium condition with ZnSb. . . . . . . . . . . . . . . . . . . . . . . . . . . . . . . . . . . . . 2.1 4 Conventional hexagonal unit cell of Zn4 Sb3 . The Sb sub-lattice (orange) is composed of isolated Sb and Sb dimers. The Zn ‘A’ sub-lattice (blue) is 90% occupied. The ‘B’ (green), ‘C’ (brown), and ‘D’ (pink) Zn sublattices are 5% occupied. The primitive rhombohedral cell is 1/3 the volume of the conventional cell. . . . . . . . . . . . . . . . . . . . . . . 2.2 7 Formation enthalpy at 0 K for Zn4 Sb3 configurations with respect to Zn and ZnSb. The solid black line represents the convex hull of structures confined to the Zn4 Sb3 lattice. Supercell calculations in green were done using the conventional unit cell, each composed of three primitive unit cells. The dashed line shows the predicted energy for supercells following the independent cells approximation. . . . . . . . . . . . . . . . . . . 12 x 2.3 Formation free energy due to phonons with respect to Zn and ZnSb of the lowest energy Zn4 Sb3 configurations at each composition. There is a more favorable contribution to the formation energy for the more disordered Zn-rich configurations, evidenced by a steeper slope with temperature. Two configurations of Zn13 Sb10 are in good agreement to 1 meV/atom. . . . . . . . . . . . . . . . . . . . . . . . . . . . . . . 2.4 13 Single-phase stability region for Zn4 Sb3 . Solid black line corresponds to the charge balanced composition. Stable compositions to the left of the line will result in a p-type semiconductor. Inset shows a close-up of the region of retrograde Zn solubility. The orange band corresponds to compositions estimated from charge carrier measurements [3]. . . . . . 2.5 15 Possible phase stability regions of Zn4 Sb3 and Zn8 Sb7 under conditions favorable to the formation of Zn8 Sb7 . The dotted line represents the phase boundary of Zn4 Sb3 in the absence of Zn8 Sb7 . . . . . . . . . . . 3.1 18 Structures of the known room temperature zinc antimonide phases, 3− shown with zinc in blue, Sb4− in orange. The 2 dimers and isolated Sb anionic structure of ZnSb is exclusively Sb4− 2 dimers, while Zn4 Sb3 , and Zn8 Sb7 contain both Sb3− and Sb4− 2 moieties. . . . . . . . . . . . . . . 3.2 21 Electron density difference plots showing networked structures of diamond moieties of ZnSb, Zn4 Sb3 and Zn8 Sb7 . Yellow isosurfaces show regions of electron accumulation with respect to unperturbed overlapping atomic wave functions. In all structures, Sb atoms share a local maximum “banana” with two closely spaced Zn atoms, forming diamondshaped moieties. The network of diamonds becomes increasingly complex in Zn4 Sb3 and Zn8 Sb7 , providing a rich source of entropy. . . . . 28 xi 3.3 Top: Zn8 Sb7 lattice sites after applying all 12 symmetry operations of the P 6/m space group. The orange Sb sites map relatively closely to the original sites. Many of the blue Zn sites map to new sites, while a few map to nearly symmetrically equivalent sites. Bottom: Idealized lattice sites after removing nearly equivalent sites. The original 28 Sb sites and 32 Zn sites are in orange and blue, respectively. An additional 66 unique Zn interstitial sites are shown in pink. 3.4 . . . . . . . . . . . . 30 (a) Calculated density of states for of ZnSb, Zn4 Sb3 and Zn8 Sb7 , (b) comparison near the Fermi level and, (c) atomic orbital projected density of states for Zn8 Sb7 . . . . . . . . . . . . . . . . . . . . . . . . . . . . . 3.5 31 Calculated Seebeck coefficient as a function of (hole) carrier concentration assuming a rigid band model for ZnSb, Zn4 Sb3 , and Zn8 Sb7 . (Dashed lines correspond to calculations from a single parabolic band model with effective carrier mass m∗ under the constant relaxation time approximation.) 3.6 . . . . . . . . . . . . . . . . . . . . . . . . . . . . . . 32 Calculated Seebeck coefficient as a function of temperature at constant carrier concentration (assuming extrinsic doping regime). Carrier concentration is chosen to maximize the power factors of the respective phases under the constant relaxation time approximation. . . . . . . . 3.7 33 Formation enthalpy with respect to Zn and ZnSb at 0 K for Zn8 Sb7 (red) and Zn4 Sb3 (blue) configurations. The lines connect the respective lowest energy configurations at each composition. . . . . . . . . . . . 34 xii 3.8 Calculated formation free energy curves for Zn4 Sb3 and Zn8 Sb7 with respect to ZnSb and Zn at three temperatures. A convex hull is drawn connecting the free energy curves of the thermodynamically stable phase at a given temperature. At 650 K, Zn and ZnSb are in equilibrium over the entire composition range shown. At 900 K, Zn4 Sb3 has broken the convex hull and has a pure phase stability window on a narrow composition range. At 1200 K, Zn8 Sb7 has broken the convex hull and shows a favorable formation energy with respect to ZnSb and Zn4 Sb3 (connected by the black dashed line). 4.1 . . . . . . . . . . . . . . . . . . 36 (left) Layered crystal structure of AZn2 Sb2 , space group P-3m1. The cations are shown in green, Zn in blue, and Sb in orange. (right) The experimental Hall carrier concentration of AZn2 Sb2 compounds increased with increasing cation electronegativity. [4, 5, 6, 7, 8, 9, 10, 11, 12, 13] 4.2 40 Electronic structure of AZn2 Sb2 compounds along symmetric directions in reciprocal space, shows similar valence band structure. An antiferromagnetic supercell was used for EuZn2 Sb2 so the Brillouin Zone is different. . . . . . . . . . . . . . . . . . . . . . . . . . . . . . . . . . 4.3 44 (top) Defect formation enthalpy for various point defects in CaZn2 Sb2 as a function of electronic Fermi level. The slope of each segment (+1, 0, −1 or −2) corresponds to the charged state of the defect. For clarity, line segments are only shown for the lowest energy charge state of each defect type, but in reality extend beyond each intersection. In all AZn2 Sb2 compounds, cation vacancies, VA, were found to be the lowest energy defect. (bottom) With increasing cation electronegativity, the formation enthalpy, ∆Hdef ect for cation vacancies in AZn2 Sb2 compounds decreases. The colored dashed lines correspond to estimated conduction band minima for each compound. . . . . . . . . . . . . . . 46 xiii 4.4 Range of phase widths predicted for A1−x Zn2 Sb2 . The boundary line corresponds to the maximum vacancy concentration for A1−x Zn2 Sb2 compounds in thermodynamic equilibrium with ZnSb. . . . . . . . . . . . . A.1 47 Zinc-blende (left) and wurtzite (right) crystal lattices. Blue cations can be occupied by Al, Ga, and In. Orange anions can be occupied by N, P, and As. . . . . . . . . . . . . . . . . . . . . . . . . . . . . . . . . . . . A.2 Selected pair (left) and triplet (right) cluster motifs on the zinc-blende lattice. . . . . . . . . . . . . . . . . . . . . . . . . . . . . . . . . . . . . A.3 51 52 Pseudo-binary convex hulls for all [Al,Ga,In][N,P,As] combinations on the zinc-blende lattice. Green points are predicted configurational energies based on a 2-component cluster expansion fit to red points. A.4 . . . 53 Cluster expansion fit error for 6-component III–V zinc-blende compounds. Formation energy is with respect to the 9 binary end-members. In the least-squares fit, points are weighted with a Boltzmann factor according to the formation energy of each configuration. The inset shows an enlarged view of the configurations with very low formation energy. A.5 . 55 Cluster expansion fit for pseudo-ternary III–N wurtzite compounds using a Bayesian fitting approach. ECIs are reported for 2000 pair, triplet, and quadruplet clusters. Lines correspond to the Bayesian prior assumption of ECI magnitude. . . . . . . . . . . . . . . . . . . . . . . . . . . . . . A.6 56 Coupled tensor-cluster interaction terms for lattice strain mapped onto the unit sphere. Lighter colors represent strong positive coupling between a particular cluster and the direction of strain. . . . . . . . . . . A.7 59 Coupled tensor-cluster interaction terms for elastic constants mapped onto the unit sphere. Lighter colors represent strong positive coupling between a given cluster and the directional elastic constants. . . . . . 59 xiv List of Tables 3.1 Bulk moduli, B, as fit to the equation of state E(V ) = 2VB0 (V − V0 )2 (for total energy, E, unit cell volume, V, and ground state volume, V0 ) and the room temperature unit cell volume (calculated and measured) for three Zn-Sb phases. . . . . . . . . . . . . . . . . . . . . . . . . . . . . . 4.1 37 Experimental Hall carrier concentration alongside calculated vacancy concentration, carrier concentration, and Hall carrier concentration for AZn2 Sb2 at 800 K. All values are in units of cm−3 . . . . . . . . . . . . A.1 Basis set for 2-dimensional strain tensor for the zinc-blende lattice up to 2nd -nearest neighbors. . . . . . . . . . . . . . . . . . . . . . . . . . . A.2 48 57 Basis set for 4-dimensional elastic tensor for the zinc-blende lattice up to 2nd -nearest neighbors. . . . . . . . . . . . . . . . . . . . . . . . . . . 58 1 Chapter 1 Introduction Materials Science can be described as the study of the relationship between material processing, structure and properties. Experiments can be designed to synthesize a material under certain conditions, observe its structure in various ways and measure its physical properties. The experimental results may validate or refute theories of how the processing affects the structure and how its structure determines the properties. In many applications, a set of desired properties may be known, but the optimal structure and processing to achieve the properties are unknown. This is the “inverse design” problem in materials science. Computational advances have enabled the theoretical study of materials to ever increasing accuracy and capability from the atomic to the continuum scale. Computational tools are now an integral component in understanding the processing-structure-property relationships necessary for the design of functional materials. In this thesis, I present a computational study of Zn-Sb-based compounds of interest for thermoelectric power generation applications. The leading character in this work is a thermodynamic function—the grand canonical potential (GCP). The GCP is capable of describing a physical system with a varying composition or number of particles. It is the natural potential to describe material systems from multi-component alloys to stoichiometric compounds that— even in their purest forms—contain defects and impurities. In the main chapters of this work, I utilize the GCP to understand the phase behavior of two complex binary Zn-Sb phases: Zn4 Sb3 (Chapter 2) and Zn8 Sb7 (Chapter 3). I also describe the intrinsic defect behavior in a series of ternary Zn-Sb-based compounds, AZn2 Sb2 2 (Chapter 4). Specific findings of these studies are highlighted below in Figure 1.1. In a broader computational effort, the Appendix contains a summary of a cluster expansion model that expands the capability of sampling large phase space alloys for use in the GCP. Additional results not included in this thesis can be found in the bibliography for thermoelectric compounds ZnSb [14], PbTe [15], Ca5 M2 Sb6 [16], Ca3 AlSb3 [16], Sr3 GaSb3 [17] and Cu2 Zn1−x Fex GeSb4 [18]. In a general form, the GCP can be expressed as a function of temperature (T) and chemical potential (µ): ! Φ (T, µ) = −kB T ln X e−(Es −µ·Ns )/kB T s with kB as Boltzmann’s constant. The partition function inside the logarithm is a sum over states with energy, Es , and vectors of the chemical potential, µ, and the number of atoms of each species, N s . For the purpose of calculating the GCP using atomistic simulations—especially with computationally intensive quantum mechanical calculations—it is useful to coarse-grain the partition function. For solid, crystalline materials, a lattice model is usually the first valid approximation. That is to say that any given atom in the system has an overwhelming likelihood to occupy a well-defined lattice site. For a phase with interstitial defects—as is the case with Zn4 Sb3 —the interstitial locations can be included as sites of the parent lattice. For a given lattice, the partition function can be further course-grained [19]:  1 Z= exp E (σ, ν, e) kB T σ∈L ν∈σ e∈ν XXX  where σ represents a configuration of atoms on the parents lattice, ν accounts for small displacements from the lattice sites (such as due to phonons), and e represents electronic excitations. A natural interpretation of this factorization comes from considering the respective time scales of each sum. An ensemble average is presumably achieved much quicker for each successive sum in the factorization of states. Specifically, for a given atomic configuration, phononic modes are readily excited and 3 sampled. Similarly—as in the Born–Oppenheimer approximation—for a given atomic perturbation, we consider the electronic states to be well-sampled. This assumption is not valid, of course, in materials that exhibit strong electron-phonon coupling. In the case that each configuration, σ, has a local potential energy minimum, the harmonic or quasi-harmonic approximation can be employed to calculate the configurational energy [19]: E(σ) −→ min {E(σ, ν, e)}} = fvib (σ) + felec (σ) ν,e The vibrational and electronic contributions are not always significant in determining thermodynamic stability, but in many cases it is an important effect [20]. In this work, I find the vibrational contribution to the free energy to have a significant stabilizing effect for Zn4 Sb3 . A further simplification to the GCP is the low-temperature expansion [21] where the system primarily occupies the ground state. This state is factored out of the partition function  ! X E0 kB T ∆Es − µ · ∆Ns φ (T, µ) = − µ · x0 − ln m0 + ms exp − N N kB T s with a ground state energy, E0 , and composition, x0 . A finite unit cell size, N , should be large enough to describe its full configurational space. Under the “independent cells” approximation, the configurations in the system are assumed to behave independently from neighboring cells. This assumption is validated by comparing the energetics of larger supercells with that of the constituent unit cells. The sum inside the logarithm includes excited states with respective excitation energies, ∆Es , and relative change in each species of atoms, ∆Ns . The multiplicative factor, ms , corresponds to the number of symmetrically equivalent states are in the cell. In this approximation, excited states are sampled independently from one another at finite temperature. In the dilute limit of point defects, there may only be a single likely defect state in the system. In this case, the defect can also produce localized electronic states. Electrons occupying these states are in equilibrium with the electronic 4 1200 10 1100 5 Sr Temperature (K) 1000 900 1000 800 900 700 800 600 700 500 12.985 400 12.4 12.99 12.995 12.5 12.6 13 −5 T=900K −10 12.8 X 10 12.9 13 13.1 A = Yb 700 600 500 ZnSb + A1-xZn2Sb2 A1-xZn2Sb2 Zn8Sb7 400 −20 Zn4Sb3 12.7 X in Zn Sb 800 T=650K 0 −15 Eu Ca Temperature (K) Formation Free Energy (meV/atom) 900 0.52 0.53 0.54 0.55 0.56 Zn Atomic Fraction 0.57 T=1200K 300 0.58 ZnSb 0.59 0.04 0.03 0.02 x in A1-xZn2Sb2 0.01 0 Figure 1.1: Summary of significant results for Zn4 Sb3 , Zn8 Sb7 , and AZn2 Sb2 phases. (Left) Single-phase stability region for Zn4 Sb3 , with a predicted region of retrograde Zn solubility [1]. (Center) Free energy curves for Zn4 Sb3 and Zn8 Sb7 with respect to Zn and ZnSb phases [2]. Zn8 Sb7 is only predicted to be stable in a small hightemperature window. (Right) Stable cation vacancy concentration in AZn2 Sb2 for A = Ca, Sr, Yb, and Eu. Boundary lines correspond to the equilibrium condition with ZnSb. chemical potential (or Fermi level), EF . In this limit, for a single defect state, D,  ! X ∆ED,q + qEF − µ · ∆N D kB T E0 mD exp − − µ · x0 − ln 1 + φ (T, µ) = N N kB T D,q where q is the local charge of the defect state. In this way—for a semiconducting material—the defect concentration is strongly coupled to the Fermi level and hole and electron carrier concentrations. In this work, I explain the defect-driven electronic trends in AZn2 Sb2 for A = Ca, Sr, Yb and Eu. When the low-temperature expansion or the dilute limit of the GCP is not valid, the lattice cluster expansion (CE) can be an efficient means to sample the large configurational space of an alloy [22]. The cluster expansion is a generalized Ising model [23, 24] where the site interaction terms are fit to a database of configurational energies—known, for instance, from density functional calculations. The computationally intractable total energy functional of atomic coordinates is simplified to a polynomial expansion of site occupation variables, σ: E(r1 , r2 ...rN ) −→ E(σ1 , σ2 ...σN ) = X Jα σα α where Jα is the effective cluster interaction (ECI) and σα is a product of occupation 5 variables included for a given cluster motif, α. An expansion can include cluster motifs of single sites, pairs, triplets and so on. The CE Hamiltonian can be utilized on systems of thousands of atoms for statistically accurate Monte Carlo simulations. To fit a CE, known configurational energies are typically added to the database while the ECIs are fit using an increasing number of cluster motifs. A fit may be considered converged if the strength of ECIs decays close to zero with longer range cluster motifs and the predictive power of the CE can quantified with cross-validation techniques. Achieving a converged CE is increasingly difficult for large multi-component alloy systems where the configurational phase space becomes exponentially larger. In Chapter A, I summarize work to fit a 6-component CE for a III–V alloy system (Al, Ga, In cations and N, P, As anions). An additional aspect of this work was to extend the CE beyond scalar quantities to include tensorial quantities such as lattice strain and elastic constants. In this way, the GCP can include various stress and strain conditions such as might be present in epitaxial materials and hetero-structures:   φ(T, µ, ) = − kB T X  exp − ln  N s Es − µ · Ns + P i,j σs,ij ij + 21 kT P i,j,k,l Cs,ijkl ij kl    where  is the 2-dimensional strain tensor, C, is the 4-dimensional elastic tensor, and σ is the 2-dimensional stress tensor. These tensor properties are fit to a cluster expansion using an additional tensor basis, β, coupled with each cluster motif [25, 26]: Q(σ) = X X Jα,β βσα β∈C(α) Appendix A includes work to fit CEs for the lattice strain and elastic tensors for the same III–V alloy system above. 6 Chapter 2 Retrograde Solubility and Entropic Stabilization in Zn4Sb3 Reproduced with permission from Phys. Rev. B 83, 094106 (2011). Copyright © 2011 American Physical Society. 2.1 Introduction The Zn-Sb binary phase system has been of interest for many years in the search for efficient and low-cost thermoelectric materials. Of primary interest has been the Zn4 Sb3 phase which exhibits a thermoelectric figure of merit, zT , in excess of 1 in intermediate temperature ranges [27]. This phase, being composed of environmentally benign and relatively earth abundant elements, continues to draw active research. Zn4 Sb3 exhibits exceptionally low lattice thermal conductivity [27], due in part to its high configurational disorder. The room temperature structure of Zn4 Sb3 (R3̄c, Figure 2.1), denoted ‘β-Zn4 Sb3 ’ here [27] (or sometimes ‘-Zn4 Sb3 ’ [28]), contains an anionic framework composed of 30 Sb, divided between 6 Sb4− 2 dimers and 18 isolated Sb3− . With 39 Zn2+ a charge balanced composition is obtained at Zn13 Sb10 . Of these Zn cations, at most, 36 can fit on Zn framework sites and 3 must be distributed among 3 crystallographically distinct interstitial sites [29]. For clarity, we will always refer to the entire phase as ‘Zn4 Sb3 ’ and write ‘ZnX Sb10 ’ only when referring to a specific configurational composition within ‘Zn4 Sb3 ’. Previous ab initio studies of Zn4 Sb3 have included the electronic density of states 7 Figure 2.1: Conventional hexagonal unit cell of Zn4 Sb3 . The Sb sub-lattice (orange) is composed of isolated Sb and Sb dimers. The Zn ‘A’ sub-lattice (blue) is 90% occupied. The ‘B’ (green), ‘C’ (brown), and ‘D’ (pink) Zn sub-lattices are 5% occupied. The primitive rhombohedral cell is 1/3 the volume of the conventional cell. 8 and formation energy of various atomic configurations ranging in composition from Zn12 Sb10 to Zn14 Sb10 [30, 31, 32, 33]. These results confirm the expectation from charge counting that Zn13 Sb10 has a Fermi level in the band gap, Zn-deficient Zn12 Sb10 has a Fermi level in the valence band and Zn-rich Zn14 Sb10 has a Fermi level in the conduction band. All reports of formation energy [31, 32, 33] also show agreement that all Zn4 Sb3 configurations have a positive formation enthalpy at 0K with respect to Zn and ZnSb. Yet Zn4 Sb3 is observed below 700K. Furthermore, at low temperature, a reversible transition to a low symmetry, meta-stable, ordered ‘α’ (and α0 ) phase is observed [34, 35]. Synthetically, Zn4 Sb3 is usually prepared with excess Zn: Zn13.3 Sb10 . As Zn has a high vapor pressure even at the relatively low synthesis and operation temperature, as well as the tendency to readily oxidize [36], it is not very surprising that excess Zn helps form the phase. Zn4 Sb3 is always a heavily doped p-type semiconductor with room temperature (Hall) carrier concentration between 6−9x1019 /cm3 . [3] According to the simple charge counting and ab initio calculations, this would be best explained by a slight Zn deficiency, Zn13−δ Sb10 , where δ should be between 0.016 and 0.024. If n/nHall is 2 as predicted by Singh [30], then δ should be between 0.008 and 0.012. It is perhaps surprising that with all the interstitial Zn sites, synthesis with excess Zn does not produce n-type material. In this work, we expand on previous ab initio results to consider a thorough thermodynamic investigation of Zn4 Sb3 . A thermodynamic ensemble approach is necessary because no single configuration (or even small number of configurations) can adequately describe the phase. By assembling a Grand Canonical Partition Function we are able to quantitatively prove entropic stabilization. We also predict a region of single-phase stability on a temperature versus composition diagram that exhibits retrograde solubility of Zn with δ > 0 that would always produce p-type Zn4 Sb3 in intermediate temperatures. 9 2.2 Computational Methods The natural thermodynamic function to account for a range of possible compositions is the grand canonical potential [21] (GCP): X kB T ln e−(Es −µ· Ns )/kB T φ (T, µ) = − N s ! (2.1) where k is Boltzmann’s constant, T is temperature, N is the total number of atoms in the system, while, for a given state s, Es is the total energy of state, Ns is a vector containing the number of atoms for each chemical species (with elements summing to N ) and µ is a vector containing the chemical potential of each species. This greatly simplifies for both Zn and ZnSb, which do not exhibit configurational disorder: φ (T, µ) = ε0 − µ · x0 n (2.2) where ε0 is the energy per unit cell, n is the number of atoms per unit cell, and x0 is a vector of the atomic fraction of each species. For Zn4 Sb3 , we consider the GCP within an “independent cells” approximation. We assume that each primitive unit cell of 23 ± 2 atoms is non-interacting with neighboring cells, in the sense that the defect configuration present in one cell does not affect the energies of defect configurations in a nearby cell. (All our ab initio calculations are nevertheless performed on infinite periodic systems with suitable k-point sampling.) This approach is useful in this system because the unit cell is rather large. This assumption is validated by computing the energy of supercells with different configurations in each primitive unit cell and comparing it to that predicted by summing the energies of the constituent primitive cells. Under the independent cells approximation: X kB T ε0 φ (T, µ) = −µ· x− ln 1 + mi e−(∆εi −µ·∆ ni )/kB T n n i>0 ! (2.3) where ε0 is the ground state energy per unit cell, n is the number of atoms per unit 10 cell and x0 is the ground state composition. For each configuration, i, mi is the symmetric multiplicity, ∆εi is the change in energy from the ground state and ∆ni is the change in the number of atoms from the ground state. The energy for each configuration is calculated under the generalized gradient approximation (GGA) using the projector augmented wave (PAW) method with Perdew-Burke-Ernzerhof (PBE) potentials as implemented in VASP 4.6, neglecting spin-orbit coupling. All unit cell parameters and atomic positions were allowed to relax to find the lowest energy configuration to within 10−4 eV. A final static calculation was performed for an accurate total energy. Defect configurations were systematically generated by enumerating defect combinations deviating from the undefected Zn12 Sb10 structure (‘A’ sub-lattice fully occupied and ‘B’,‘C and ‘D’ sub-lattices fully unoccupied). Defect clusters were composed of a combination of vacancies on the ‘A’ sub-lattice, and occupation of an interstitial site on the ‘B’, ‘C’ or ‘D’ sub-lattices. Clusters of up to 6 defects were enumerated, excluding some structures because the defects were too close (nearest neighbor A and C sites both occupied). After allowing the atomic configuration to relax, the resulting atomic positions were projected onto the closest unrelaxed configuration and any duplicate configurations were excluded in order to avoid over-counting states. Phonon density of states and vibrational free energies were calculated using the ‘supercell’ method as implemented in the Alloy Theoretic Automated Toolkit (ATAT) [37, 38, 39]. Since the computational resources needed to compute phonon modes for all Zn4 Sb3 configurations are prohibitive, representative configurations were selected at four compositions between Zn12 Sb10 and Zn15 Sb10 configurations. The lowest energy configuration for each composition was selected as well as a second Zn13 Sb10 configuration to assess the error in our approximation. For each of these configurations and the end-members Zn and ZnSb, phonon modes were calculated at 0%, 1% and 2% strain to account for the effects of thermal expansion, under the quasiharmonic approximation. The vibrational contribution to the free energy is incorporated into the GCP through a nested sum in the partition function [19]. For each distinct configura- 11 tion, phonon occupation accounts for small displacements around the local energy minimum, resulting in a temperature dependent free energy correction. GCPs were assembled for Zn, Zn4 Sb3 and ZnSb. Phase equilibrium is determined by equality of two respective GCPs. The equilibrium composition of each phase can be determined by: ∇µ φ (T, µ) = −x 2.3 (2.4) Results and Discussion Over 100 unique, stable configurations were enumerated in the β−Zn4 Sb3 primitive rhombohedral unit cell with between 21 and 25 atoms per cell (the number of Sb was held constant at 10). The 0K formation enthalpy of these configurations with respect to Zn and ZnSb is shown in Figure 2.2. Consistent with previous theoretical studies [31, 32, 33], all configurations have a positive formation enthalpy with respect to Zn and ZnSb. The “convex hull” is plotted in solid black and connects the ground state configurations (for Zn4 Sb3 phase) at each composition. The composition Zn14 Sb10 does not have a configuration touching the convex hull. This means that for an ensemble of atoms of composition Zn14 Sb10 constrained to remain in the Zn4 Sb3 lattice, it would be more energetically favorable to form a mixture of cells of composition Zn13 Sb10 and Zn15 Sb10 . If the system is allowed to adopt any lattice, then, at 0K, it would be even more energetically favorable to separate into Zn and ZnSb crystals. Figure 2.2 also allows us to test the independent cells approximation. The data points marked “supercells” are all made up of 3 primitive cells of different compositions to yield the intermediate compositions shown. Between Zn12 Sb10 and Zn13 Sb10 , we see that the supercells produce a slightly lower energy by about 3 meV/atom than predicted by the independent cells approximation. Between Zn13 Sb10 and Zn14 Sb10 , the energies from the supercells are either very close or slightly higher than predicted by the independent cells approximation (the dashed line). The small energy difference of a few meV/atom justifies the use of the independent cells approximation, although a consideration of a possible systematic error 12 100 90 80 E f (meV/atom) 70 − − − − − − − − 60 50 30 20 0 − − − − − − − − − − − − − − − − − − − − − − − − − 40 10 Supercells − − − − − − − − − − − Zn 11 Sb 10 Zn 12 Sb 10 Zn 13 Sb 10 Zn 14 Sb 10 Zn 15 Sb 10 0.524 0.546 0.565 0.583 0.600 Zn Atomic Fraction Figure 2.2: Formation enthalpy at 0 K for Zn4 Sb3 configurations with respect to Zn and ZnSb. The solid black line represents the convex hull of structures confined to the Zn4 Sb3 lattice. Supercell calculations in green were done using the conventional unit cell, each composed of three primitive unit cells. The dashed line shows the predicted energy for supercells following the independent cells approximation. 13 80 70 Zn Sb Ef + Fvib (meV/atom) 60 15 10 50 Zn 14 Sb 10 40 30 Zn 12 Sb 10 20 10 Zn 13 Sb 10 0 0 200 400 600 Temperature (K) 800 1000 Figure 2.3: Formation free energy due to phonons with respect to Zn and ZnSb of the lowest energy Zn4 Sb3 configurations at each composition. There is a more favorable contribution to the formation energy for the more disordered Zn-rich configurations, evidenced by a steeper slope with temperature. Two configurations of Zn13 Sb10 are in good agreement to 1 meV/atom. is discussed later. The formation free energy due to phonons (relative to Zn and ZnSb) is shown in Figure 2.3 for the 5 representative configurations (the lowest energy structure at each composition). In all cases, the phonon contribution to the free energy favors the Zn4 Sb3 phase over Zn and ZnSb as evidenced by a negative trend with temperature. This is likely a result of the softer phonon modes in the more complex and open crystal structure. Increasing Zn concentration and thereby increasing the Zn-disorder results in a more favorable contribution to the free energy. There is good agreement between the two configurations of composition Zn13 Sb10 , differing by 1 meV/atom at 1000K. It then seems reasonable to assume that the vibrational free energy of the representative structures may be applied to all the configurations at that composition 14 when we compute the GCP of Zn4 Sb3 . The computed region of single-phase stability for Zn4 Sb3 is shown in Figure 2.4. Including configurational and vibrational effects to the free energy, Zn4 Sb3 is found to stabilize at around 700K at a composition of Zn12.992 Sb10 . As the temperature increases, the range of stable single-phase compositions increases, more broadly on the Zn-deficient side. (Eventually, the solid melts, but our analysis focuses on the solid-state portion of the phase diagram.) Interestingly, on the Zn-rich side, we observe retrograde solubility of Zn. That is, the highest stable concentration of Zn in Zn4 Sb3 decreases with increasing temperature. This unusual finding is a result of the composition-dependent relationship between the enthalpy and entropy of mixing. Typically defects of all kinds entropically stabilize stoichiometric variations on both sides of a valence-precise composition. In the case of Zn4 Sb3 , the valence-precise composition is already ‘defected’ with interstitial Zn. A miscibility gap arises for compositions above Zn13 Sb10 due to the large jump in formation enthalpy compared to the Zn-deficient compositions. The effect is lessened by the counteractive effects of higher vibrational entropy for the Zn-rich compositions. If the vibrational effects were neglected, the Zn-rich boundary would not curve to the right until a much higher temperature. The region of single-phase stability shown in Figure 2.4 has readily observable consequences. The black vertical line at composition Zn13 Sb10 represents the valenceprecise structure with a Fermi level inside the band gap. At any Zn-deficient composition from Zn13 Sb10 , we expect a partially filled valence band resulting in a p-type semiconductor. (Surprisingly, there is a recent report of a Zn13 Sb10 configuration being slightly p-type [40]. Our electronic density of states results are in agreement with previous studies [30, 31, 32, 33].) All stable compositions below 1000K in Figure 2.4 result in a p-type semiconductor. Attempting to dope Zn4 Sb3 with excess Zn will not result in an n-type semiconductor but a two-phase equilibrium between Zn metal and p-type Zn4 Sb3 . Furthermore, with the consideration of retrograde solubility, if Zn-saturated single-phase Zn4 Sb3 is heated up, it passes through a region where it becomes thermodynamically more stable to precipitate Zn metal. With more heating, 15 1200 1100 Temperature (K) 1000 900 1000 800 900 700 800 600 700 500 12.985 400 12.4 12.99 12.995 12.5 12.6 13 12.7 12.8 X in Zn Sb X 12.9 13 13.1 10 Figure 2.4: Single-phase stability region for Zn4 Sb3 . Solid black line corresponds to the charge balanced composition. Stable compositions to the left of the line will result in a p-type semiconductor. Inset shows a close-up of the region of retrograde Zn solubility. The orange band corresponds to compositions estimated from charge carrier measurements [3]. 16 the Zn is reabsorbed as more Zn-rich compositions are stable. Upon cooling, the same precipitation and absorption should occur if held in thermodynamic equilibrium. On the Zn-deficient side of the stable phase region the usual temperature dependent solubility is observed. Here one would expect ZnSb to precipitate as the temperature is reduced. Because this would occur at low temperature in the solid state, small nanometer sized precipitates would be expected. Such particles are indeed observed in Sb-rich samples [3]. The high-temperature phase boundary on the Zn-rich side could explain the formation of Zn nanoparticles observed in some Zn4 Sb3 samples [41]. Cooling a Zn-rich composition above Zn13 Sb10 would precipitate Zn (possibly nano-particles if cooled fast enough [42]). Upon further cooling through the retrograde region, some of the Zn would be reabsorbed into the β−phase. The absorption of nanoparticle Zn may explain the nano-voids observed in some samples [3]. The precipitation of Zn from the Zn4 Sb3 structure and then subsequent re-absorption also explains the increase in volume (leading to breakage of the glass ampoule) in some samples, that at room temperature appear to be single-phase. Experimental carrier concentration measurements from Toberer [3] from singlephase β−Zn4 Sb3 samples range from 6 − 9x1019 /cm3 . This corresponds to a composition range of 0.002% atomic Zn between Zn12.988 Sb10 and Zn12.992 Sb10 assuming n/nHall = 2 [30]. This concentration range, corresponding to the orange band in the inset of Figure 2.4, is in close agreement with the predicted stable compositions near the stabilization temperature. This composition range is much smaller than the 0.2% atomic Zn observed by microprobe analysis. This discrepancy is not surprising since the composition range is basically the limit of the microprobe resolution. There is some expected uncertainty in our predicted temperature of stabilization of Zn4 Sb3 . Our predicted temperature of 700K is much higher than might be expected. There are several possible explanations for this. One simple possibility is the uncertainty inherent in first principles calculations using the GGA-PBE functional. There could also be a deficiency in our methodology to fully account for the sources of entropy in such a complex disordered structure. Our enumeration method may un- 17 derestimate the number of Zn4 Sb3 configurations by comparing the relaxed structures to the fixed lattice in Figure 2.1, which represents, in a sense, a measured average over many possible configurations. The vibrational contribution to the free energy has a significant effect in lowering the stabilization temperature. Assuming a representative structure for each composition might be an underestimate for the more disordered configurations and the effect could be even more pronounced. Lastly, we noted earlier that the supercells on the Zn-deficient side of Zn13 Sb10 had slightly lower formation energy than expected under the independent cells approximation. This neglected cellto-cell interaction ( 2 meV/atom) could lower the stabilization temperature by several hundred degrees. In either case, we expect the retrograde solubility to remain, and possibly intensify, because the supercells on the Zn-rich side of Zn13 Sb10 are in good agreement with the independent cells approximation. To ensure that our finding of a single-phase region for the β-Zn4 Sb3 lattice is not an artifact of the independent cell approximation, we have conducted separate Monte Carlo simulations based on a cluster expansion Hamiltonian fitted to our database of structural energies. These simulations confirm the presence of a single-phase region over the temperature range where β-Zn4 Sb3 is stable. Hence, the independent cell approximation was deemed reliable and was used throughout this work. It provides a convenient explicit expression for the free energy and is immune to the fitting errors inherent to the Hamiltonian construction procedure. Finally, we consider how Zn4 Sb3 might interact with other nearby phases, namely, α−Zn4 Sb3 and Zn8 Sb7 . First principles calculation of α−Zn4 Sb3 yields positive formation energy of 19 meV/atom with respect to Zn and ZnSb. This is 3 meV/atom below the lowest energy β−Zn4 Sb3 configuration of composition Zn13 Sb10 . Assuming no configurational disorder in α−Zn4 Sb3 we predict β−Zn4 Sb3 to become energetically favorable compared to α−Zn4 Sb3 at 300K (although both are still meta-stable with respect to Zn and ZnSb at this temperature). It seems possible that there is a temperature range in which meta-stable β−Zn4 Sb3 is observed before becoming thermodynamically stable with respect to Zn and ZnSb at some higher temperature. Another phase of recent interest is Zn8 Sb7 , which has been characterized experi- 18 1200 1100 Zn 8 Sb 7 + Zn 4 Sb 3 Temperature (K) 1000 900 800 700 Zn Sb 8 7 600 500 0.53 0.535 0.54 0.545 0.55 0.555 0.56 0.565 0.57 Zn Atomic Fraction Figure 2.5: Possible phase stability regions of Zn4 Sb3 and Zn8 Sb7 under conditions favorable to the formation of Zn8 Sb7 . The dotted line represents the phase boundary of Zn4 Sb3 in the absence of Zn8 Sb7 . mentally and studied through first principles calculations [43, 2]. These calculations reveal that, even though Zn8 Sb7 is unstable in bulk form, it could be stabilized in nanocrystalline form if surface stress or energy contributions result in a free energy decrease of 5 meV/atom (relative to ZnSb and Zn4 Sb3 ). Taking these assumptions into effect, the resulting phase diagram is shown in Figure 2.5. We see that Zn8 Sb7 would stabilize at a temperature higher than Zn4 Sb3 and a two-phase region results. This also reduces the stable Zn-deficient Zn4 Sb3 composition range predicted in the absence of Zn8 Sb7 (dashed line in Figure 2.5). Other unidentified phases as reported in some phase diagrams [44] could have a similar effect. 19 2.4 Conclusions From our first principles investigation, we have shown that Zn4 Sb3 is entropically stabilized with the help of configurational and vibrational entropy. Under the independent cells approximation of the grand canonical potential we predict a region of single-phase stability near Zn12.992 Sb10 , which results in a nominally p-type semiconductor. Additionally, we predict a temperature range with retrograde Zn solubility. The temperature dependent solubility can be used to explain the variety of nanoparticle formation observed in the system: formation of ZnSb on the Sb-rich side, Zn on the far Zn-rich side and nano-void formation due to Zn precipitates being reabsorbed at lower temperatures. 20 Chapter 3 Predicted Thermodynamic and Electronic Properties of a Newly Discovered Binary Zn8Sb7 Compound Reproduced with permission from J. Am. Chem. Soc., 2011, 133 (29), pp 11255–11261. Copyright © 2011 American Chemical Society. 3.1 Introduction Intermetallic compounds are fascinating due to the diversity of their crystal structures and resulting properties. For thermoelectric applications, ZnSb and Zn4 Sb3 are two intermetallic compounds which exhibit good thermoelectric efficiencies. In particular, Zn4 Sb3 has a thermoelectric figure of merit (zT ) in excess of unity [27, 3] due to its low thermal conductivity, which arises in part due to its complex, defect-ridden structure. Both of these materials have band gaps at or near the Fermi level which can be rationalized from valence charge counting [29]. Recently, a new binary Zn-Sb phase between ZnSb and Zn4 Sb3 was discovered and characterized at room temperature by Birkel et al. [43] Binary nanoparticles were synthesized in solution through the controlled reaction of elemental Zn and Sb nanoparticles. The small particle size of the reactants ensures minimum diffusion paths, low activation barriers, and low reaction temperatures, thereby eliminating solid-solid diffusion as the rate-limiting step as found in conventional bulk-scale solid- 21 ZnSb Zn8Sb7 [001] Zn4Sb3 Zn8Sb7 [100] Figure 3.1: Structures of the known room temperature zinc antimonide phases, shown 3− with zinc in blue, Sb4− in orange. The anionic structure of 2 dimers and isolated Sb 4− ZnSb is exclusively Sb2 dimers, while Zn4 Sb3 , and Zn8 Sb7 contain both Sb3− and Sb4− 2 moieties. state synthesis [43]. In the present work, the authors denote this new phase as Zn8 Sb7 , based on the idealized crystal structure, rather than Zn1+δ Sb (with δ = 0.068) as used in by Birkel et al. [43] Determination of the Zn8 Sb7 structure was a tour de force of automated electron diffraction tomography (ADT) combined with precession electron diffraction (PED) on individual nanoparticle grains. This work yielded the structure in Figure 3.1, shown along the [001] and [100] directions in the distorted hexagonal cell (P 1̄). Much 3− like Zn4 Sb3 , Zn8 Sb7 contains both Sb4− 2 dimers and isolated Sb , however, the long- range atomic arrangement is completely distinct. The Zn8 Sb7 phase enters into an interesting region of the binary phase diagram 22 [44] between ZnSb and Zn4 Sb3 . Since the three phases are so close in composition, the competition for thermodynamic stability will be fierce. Below ∼250 K, Zn4 Sb3 transforms to a variant (the α or α0 phase) where the interstitial atoms form an ordered arrangement [34, 35]. The disordered β phase is stable from ∼250 K (to 700 K) while no compounds have been characterized in the bulk at compositions between Zn4 Sb3 and ZnSb. Several phase diagrams have been reported that include a phase between Zn4 Sb3 and ZnSb stable only at high temperature [44, 45, 46], sometimes denoted as γ, but no crystal structures have been determined for these high temperature phases. Ab initio calculations have considered various interstitial configurations of Zn4 Sb3 at 0K and found that they should decompose to ZnSb and Zn when considering the enthalpy of formation alone [30, 31, 32, 33]. Additional theoretical work has recently shown Zn4 Sb3 to be entropically stabilized due to the high configurational disorder with the three interstitial Zn sites and the effects of lattice vibrations [1]. The narrow stable composition range is near Zn3.898 Sb3 as inferred from Hall carrier measurements and valence charge counting [3, 1], which is off from the valence balanced composition of Zn3.9 Sb3 . We will continue to use the conventional Zn4 Sb3 terminology. In this work, we explore the electronic and thermodynamic properties of a theoretical bulk Zn8 Sb7 phase. Electronic structure calculations probe the band structure and valence band edge character of Zn8 Sb7 to compare with known bulk thermoelectrics ZnSb and Zn4 Sb3 . Since the formation of Zn8 Sb7 is not readily observed in practice, a thermodynamic investigation determines if the phase might be thermodynamically stable or meta-stable with respect to competing Zn-Sb phases. A thorough thermodynamic treatment considers the bonding enthalpy as well as possible sources of entropy such as configurational disorder and atomic vibrations. The results will guide future experiments to synthesize bulk Zn8 Sb7 . 3.2 Computational Methods To account for the configurational disorder present in many Zn-Sb phases, a thermodynamic ensemble approach is employed using data from first principles calculations. 23 A grand canonical potential [21] (GCP) is formulated for each phase according to: X kB T ln e−(Es −µ · Ns )/kB T φ (T, µ) = − N s ! (3.1) where kB is Boltzmann’s constant, T is temperature, N is the total number of atoms in the system, while, for a given state s, Es is the total energy of state, Ns is a vector containing the number of atoms for each chemical species (with elements summing to N ) and µ is a vector containing the chemical potential of each species. Since ZnSb and Zn phases exhibit relatively low configurational disorder, the GCP greatly simplifies: φ (T, µ) = ε0 − µ · x0 n (3.2) where ε0 is the energy per unit cell, n is the number of atoms per unit cell, and x0 is a vector of the atomic fraction of each species. For Zn8 Sb7 and Zn4 Sb3 , we factor the GCP according to an “independent cells” approximation. We assume that each primitive unit cell is non-interacting with neighboring cells, in the sense that the defect configuration present in one cell does not affect the energies of defect configurations in a nearby cell. (All our ab initio calculations are nevertheless performed on infinite periodic systems with suitable k-point sampling.) This approach is useful in these systems because the unit cells is rather large. This assumption is validated by computing the energy of supercells with different configurations in each primitive unit cell and comparing it to that predicted by summing the energies of the constituent primitive cells. Under the independent cells approximation: X ε0 kB T φ (T, µ) = − µ · x0 − ln 1 + mi e−(∆εi −µ ·∆ ni )/kB T n n i>0 ! (3.3) where ε0 is the ground state energy per unit cell, n is the number of atoms per unit cell and x0 is the ground state composition. For each configuration, i, mi is the symmetric multiplicity, ∆εi is the change in energy from the ground state and ∆ ni 24 is the change in the number of atoms from the ground state. The total energy of each atomic configuration was calculated using Density Functional Theory (DFT) with the PBE exchange-correlation functional [47]. Calculations utilized the projector augmented wave (PAW) method [48] as implemented in VASP[49] neglecting spin-orbit coupling. All unit cell parameters and atomic positions were allowed to relax to find the lowest energy arrangement to within 10−4 eV per primitive unit cell. A final static calculation was performed for an accurate total energy. Energy convergence with respect to k-points was achieved with a Γ-centered grid with 14 k-points in the irreducible Brillouin zone. Configurational arrangements for Zn4 Sb3 were enumerated by considering all possible interstitial combinations on known Zn interstitial sites [29, 1]. Configurational arrangements of Zn8 Sb7 are much more difficult to enumerate appropriately without guidance from experimentally known defect lattice sites. The configurational disorder then was probed by a systematic enumeration of new possible configurations on the Zn sub-lattice. The specific enumeration process is described later in the discussion section. Configurations with Zn vacancies and interstitials deviating from the Zn8 Sb7 composition are also considered. DFT calculations were performed on select Zn8 Sb7 configurations since it would be computationally prohibitive to calculate all possibilities. Atomic positions were allowed to relax to a local minimum and were verified to not have settled into identical configurations. While this method is by no means exhaustive in gauging the configurational disorder in Zn8 Sb7 , it does provide a basis for making certain assumptions to be discussed in the results section to follow. Phonon density of states and vibrational free energies were calculated using the ‘supercell’ method as implemented in the Alloy Theoretic Automated Toolkit (ATAT) [37, 38, 39]. The computational resources needed to compute phonon modes for all Zn8 Sb7 configurations are prohibitive (360 perturbation structures per configuration per volume), so the phonon modes were only calculated for the nominal configuration. In all, phonon modes were calculated for one Zn8 Sb7 configuration, several Zn4 Sb3 configurations as well as end-members Zn and ZnSb at 0%, 1% and 2% strain to account for the effects of thermal expansion, under the quasi-harmonic approximation. 25 The vibrational contribution to the free energy is incorporated into the GCP through a nested sum in the partition function [19]. For each distinct configuration, phonon occupation accounts for small displacements around the local energy minimum, resulting in a temperature dependent free energy correction. GCPs were assembled for Zn, Zn4 Sb3 , Zn8 Sb7 and ZnSb. Phase equilibrium is determined by equality of two respective GCPs. The equilibrium composition of each phase can be determined by: ∇µ φ (T, µ ) = −x (3.4) and the free energy is calculated as: f (T, x) = φ (T, µx ) + µx · x (3.5) where µx is the chemical potential that stabilizes the phase at composition x . The Seebeck coefficient is calculated for the ground state structures using the Boltzmann transport equation within the constant relaxation time approximation as implemented in the BoltzTraP code [50]. A dense k-point mesh is required to evaluate the Fermi-integrals over the entire Brillouin zone. Convergence of the Seebeck coefficient with respect to carrier concentration is achieved in Zn8 Sb7 with a 14x14x14Γ-centered k-point grid. For Zn4 Sb3 and ZnSb structures, 16x16x16 and 40x32x32 k-point grids are used, respectively. 3.3 Results and Discussion Crystal Structure and Charge Counting The Zn-Sb structures in Figure 3.1 can all be described using the Zintl-Klemm formalism [51], in which Zn cations donate electrons to Sb anions which may then form covalent bonds with other Sb to complete their valence shells. A molecular orbital (MO) approach, which includes interactions between Zn and Sb can also successfully rationalize the electron count [33]. In ZnSb (P bca), all the respective Zn and Sb atoms 26 are symmetrically equivalent. The Zn cations each donate 2 electrons to an Sb anion 2+ and the Sb atoms form Sb4− 2 dimers. A valence-precise structure results with 8 Zn and 4 Sb4− 2 per unit cell. One then can rationalize why ZnSb behaves electronically as an intrinsic semiconductor with a band gap between the valence band (dominated by Sb p-states) and the conduction band (composed primarily of Zn s-states). The Zn4 Sb3 (R3̄c) framework structure is composed of 30 Sb (orange sites in Figure 3.1) and 36 Zn (blue sites) in the conventional unit cell. The anionic frame3− work is divided between 6 Sb4− 2 dimers and 18 isolated Sb . A fully occupied cation framework of 36 Zn2+ is insufficient to satisfy a valence charge balance. Three additional Zn2+ are needed to satisfy the charge balance at a composition of Zn13 Sb10 . These additional Zn atoms can occupy any of three symmetrically distinct regions for interstitial sites (pink, green and brown sites), introducing configurational disorder not present in ZnSb. Thermodynamically, Zn4 Sb3 is only stable in a narrow range of composition near Zn12.992 Sb10 , which is slightly valence deficient. This explains the persistent p-type doping observed experimentally in Zn4 Sb3 . Examining the Zn8 Sb7 (P 1̄) structure reveals some similarities to Zn4 Sb3 . The anionic Sb framework is only slightly distorted from a high symmetry structure that can be described with a P 6/m space group. There are 10 Sb4− 2 dimers and 8 isolated Sb3− , all sitting close to points of symmetry. The cation Zn framework exhibits low symmetry with 32 Zn2+ . This satisfies a precise valence condition, suggesting that Zn8 Sb7 , like ZnSb and Zn13 Sb10 , should behave as an intrinsic semiconductor. We further investigate the structure of Zn8 Sb7 with the help of electron density difference (EDD) plots. The EDD compares the calculated electron charge density with the charge density resulting from the overlap of non-interacting single atom electronic wave functions. EDD plots thus reveal where charge accumulation and depletion occur as atomic wave functions interact. Charge accumulation between atoms suggests the presence of a covalent bond. Figure 3.2 shows the EDD around a diamond moiety in ZnSb, which clearly reveals charge accumulation lobes in a tetrahedral arrangement around each Sb atom. In agreement with a previous molecular orbital study [33], one lobe of each Sb is 27 shared between two closely spaced Zn atoms, forming a repeated diamond unit. This same diamond moiety is also apparent in Zn4 Sb3 . Figure 3.2 shows the EDD for a configuration of composition Zn12 Sb10 (no interstitial Zn). In this case, each lone Sb atom shares a lobe with each of two pairs of Zn atoms, forming zigzag chains of diamond moieties [33]. The Sb-dimers form a separate linear chain from the diamonds. The presence of interstitial Zn, such as in Zn13 Sb10 configurations, breaks up the chain and creates local disorder in the diamond network. High configurational entropy results because there are many possible ways for the interstitial Zn to distort and rearrange the diamond network, all with similar energetic consequences. In Zn8 Sb7 , we again see the diamond moiety but with added complexity. Nonlinear chains of diamond moieties incorporate both Sb-dimers (as in ZnSb) as well as lone Sb (as in Zn4 Sb3 ). All the diamonds are distorted out of plane, as is shown for two of the diamonds in Figure 3.2. The crystal structure is frustrated in the sense that it cannot settle into a single preferred ordered state, such as in ZnSb, but like Zn4 Sb3 , there are many configurational possibilities that result in similar energetics. One way to enumerate the possible configurational states is to consider the mismatch of symmetry between the Sb and Zn sub-lattices. There are 6 symmetrically equivalent ways to combine the higher symmetry P 6/m Sb sub-lattice with the lower symmetry P 1̄ Zn sub-lattice. By combining these 6 arrangements into a single lattice, one can discover possible interstitial sites on the Zn sub-lattice. The combined lattice is shown in Figure 3.3. All the orange Sb sites map closely to existing sites as expected. However, the blue Zn sites reveal many new locations. Nearly equivalent sites are combined to form an idealized lattice with 66 new possible Zn sites in pink as shown at the bottom of Figure 3.3. These sites provide a basis to probe the configurational disorder. Configurational arrangements involving Zn vacancies on the highly occupied blue sites and Zn interstitials on partially occupied pink sites can be enumerated and the energies calculated from first principles calculations. This configurational disorder could likely be undetected by electron diffraction analysis due to the low occupation of certain sites. The experimentally predicted [43] partial occupancy results in a composition of Zn29.9 Sb28 , which could be explained as an 28 ZnSb Zn8Sb7 Zn12Sb10 Zn13Sb10 Figure 3.2: Electron density difference plots showing networked structures of diamond moieties of ZnSb, Zn4 Sb3 and Zn8 Sb7 . Yellow isosurfaces show regions of electron accumulation with respect to unperturbed overlapping atomic wave functions. In all structures, Sb atoms share a local maximum “banana” with two closely spaced Zn atoms, forming diamond-shaped moieties. The network of diamonds becomes increasingly complex in Zn4 Sb3 and Zn8 Sb7 , providing a rich source of entropy. 29 average of 2.1 Zn atoms per unit cell occupying various interstitial sites. The entropy associated with this disorder will be addressed in the phase stability section. Electronic Structure The calculated electronic density of states of Zn8 Sb7 is shown in Figure 3.4 along with ZnSb and Zn4 Sb3 . As expected from the Zintl-Klemm charge counting, Zn8 Sb7 is found to be an intrinsic semiconductor like ZnSb and Zn4 Sb3 . The calculated band gap is comparable to that of ZnSb but slightly smaller than that of Zn4 Sb3 . The similar shape and magnitude of the valence band edges of Zn8 Sb7 and Zn4 Sb3 suggest a similar hole carrier effective mass. This result alone suggests that Zn8 Sb7 might have similar electronic properties to that of Zn4 Sb3 . The valence band edge is composed primarily of electronic states localized around the non-bonded Sb atoms. The projected density of states in Figure 3.4 shows that the valence band is composed primarily of Sb p-orbital character, which is similar to previous studies of Zn4 Sb3 and ZnSb [33]. The calculated Seebeck coefficient, S, as a function of (hole) carrier concentration is shown in Figure 3.5 for 700K. This assumes a rigid band model with a carrier concentration controlled by shifting the electronic chemical potential away from the intrinsic Fermi energy. For reference, dashed lines show trends derived from a single parabolic band approximation for an effective carrier mass, m∗ . The Seebeck coefficient of Zn8 Sb7 approaches that of Zn4 Sb3 at higher carrier concentration for similar doping levels. An optimum carrier concentration can be chosen to maximize the power factor, σS 2 , with electrical conductivity, σ, calculated under the constant relaxation time approximation. An approximate relaxation time is found by fitting to experimental resistivity data for ZnSb [14] and Zn4 Sb3 [3] and assuming a similar value for Zn8 Sb7 . The constant carrier concentration Seebeck as a function of temperature is shown in Figure 3.6. Both Zn4 Sb3 and Zn8 Sb7 have a maximum power factor for a hole concentration of around 5 × 1020 cm−3 and ZnSb for a concentration of 7.5 × 1019 cm−3 . It should be noted that in ZnSb it has proven difficult to 30 Figure 3.3: Top: Zn8 Sb7 lattice sites after applying all 12 symmetry operations of the P 6/m space group. The orange Sb sites map relatively closely to the original sites. Many of the blue Zn sites map to new sites, while a few map to nearly symmetrically equivalent sites. Bottom: Idealized lattice sites after removing nearly equivalent sites. The original 28 Sb sites and 32 Zn sites are in orange and blue, respectively. An additional 66 unique Zn interstitial sites are shown in pink. 31 a) ZnSb 1 States/Atom eV 0.5 0 Zn4Sb 3 1 0.5 0 Zn8Sb 7 1 0.5 0 −10 States/Atom eV b) States/Atom eV Energy (eV) 0 5 1 0.8 0.6 0.4 0.2 0 −2 c) −5 −1.5 −1 −0.5 Energy (eV) 0 0.5 0.5 1 Sb−p Zn−p Zn−d Zn−s Sb−s Sb−d 0.4 0.3 0.2 0.1 0 −2 −1.5 −1 −0.5 Energy (eV) 0 0.5 1 Figure 3.4: (a) Calculated density of states for of ZnSb, Zn4 Sb3 and Zn8 Sb7 , (b) comparison near the Fermi level and, (c) atomic orbital projected density of states for Zn8 Sb7 . 32 300 700K 250 Zn4 Sb3 S (µV/K) 200 m*=1.5 Zn8 Sb 7 150 ZnSb m*=1.0 100 50 0 m*=0.5 0.5 1 1.5 2 2.5 3 hole carrier concentration (cm−3 ) 3.5 4 x10 20 Figure 3.5: Calculated Seebeck coefficient as a function of (hole) carrier concentration assuming a rigid band model for ZnSb, Zn4 Sb3 , and Zn8 Sb7 . (Dashed lines correspond to calculations from a single parabolic band model with effective carrier mass m∗ under the constant relaxation time approximation.) achieve hole concentrations above 1019 cm−3 . In Zn4 Sb3 single-phase material, the Hall carrier concentration ranges between 6 − 9 × 1019 cm−3 . Phase Stability Finally, we consider the thermodynamic phase stability of this new Zn8 Sb7 phase. We note that for Zn8 Sb7 to be thermodynamically stable, the free energy must fall below the tie-line between ZnSb and Zn4 Sb3 . This will prove difficult for Zn8 Sb7 as Zn4 Sb3 has high entropy even at room temperature. We start by considering the formation enthalpy at 0K and then include entropic contributions at finite temperature. Figure 3.7 shows the calculated formation enthalpy of Zn8 Sb7 configurations (red) alongside those of Zn4 Sb3 (blue). We readily observe that Zn8 Sb7 configurations have positive formation enthalpy with respect to decomposition into ZnSb and Zn. This is perhaps not unexpected due to similarities with Zn4 Sb3 , which also has positive formation enthalpy at 0K. However, Zn4 Sb3 is known to be thermodynamically stable at higher temperatures [1] due to both configurational and vibrational entropic 33 200 Zn4 Sb3 ( ND = 4.9x1020 ) S (µV/K) 150 100 ZnSb ( ND = 7.5x1019 ) 50 Zn8 Sb 7 ( ND = 4.7x1020 ) 0 100 200 300 400 500 Temperature (K) 600 700 Figure 3.6: Calculated Seebeck coefficient as a function of temperature at constant carrier concentration (assuming extrinsic doping regime). Carrier concentration is chosen to maximize the power factors of the respective phases under the constant relaxation time approximation. 34 100 Formation Enthalpy (meV/atom) 90 80 − − − 70 60 50 40 30 − − − − − − − − − − − − − − − − − − − − − − − − − − − − − − − − − − Zn Sb 12 20 Zn 8 Sb 7 10 0 − − − − − − − − Zn8 Sb7 configurations 0.5 0.52 0.54 − − − − − − − − − − − − − − − − − 10 Zn 13 Sb 10 Zn4 Sb3 configurations 0.56 Zn atomic fraction 0.58 0.6 Figure 3.7: Formation enthalpy with respect to Zn and ZnSb at 0 K for Zn8 Sb7 (red) and Zn4 Sb3 (blue) configurations. The lines connect the respective lowest energy configurations at each composition. contributions to the free energy. These sources might have significant contributions for Zn8 Sb7 as well. Entropy arises in an ensemble when many configurations are energetically accessible. It is apparent from Figure 3.7 that we will expect some configurational entropy in Zn8 Sb7 at finite temperature from the several configurations with similar formation energy. Due to computational limitations on the large unit cell, we have only calculated 30 different configurations at the Zn32 Sb28 composition as well as 30 offstoichiometric configurations. These known configurations represent only a small subset of the possible arrangements of atoms that can be enumerated based on the idealized lattice in Figure 3.3. Each of our known configurations is likely representative of many other configurations of similar energy. To estimate this, we consider the ‘mean field’ entropy based on the average fractional occupancy, xi of each site, i: 35 Smf [kB /cell] = −kB X [xi ln(xi ) + (1 − xi )ln(1 − xi )] (3.6) i We can correct this value by assuming that Zn8 Sb7 deviates from the mean field value similarly to Zn4 Sb3 . For Zn4 Sb3 , considering 90% occupancy on the framework sites and 5% occupancy on each of 3 interstitial sites gives an average mean field entropy of 0.48 kB /atom. This is almost twice as high as the configurational entropy calculated from the grand canonical ensemble, 0.21 kB /atom. For Zn8 Sb7 , we consider a Zn framework site occupancy of 93% based on the experimentally determined partial occupancy. The remaining Zn is interspersed over 66 possible interstitial sites. This results in a mean field entropy of 0.28 kB /atom, which we reduce by the correction factor for Zn4 Sb3 to estimate the configurational entropy for Zn8 Sb7 as 0.12 kB /atom. The grand canonical potential (Equation 3.3) for Zn8 Sb7 is tuned by increasing the ∂φ multiplicative constants, mi , so that, in the high-temperature limit, S = − ∂T = 0.12 kB /atom (the adjusted mean field entropy). In essence, this procedure combines enthalpy obtained from the ab initio calculations with an entropy estimate based on structural similarities with the Zn4 Sb3 phase. The vibrational contributions to the free energy have also been accounted for in the quasi-harmonic approximation. Both Zn4 Sb3 and Zn8 Sb7 have a favorable contribution to the free energy from phonons with respect to Zn and ZnSb. In the high-temperature limit above 300 K, Zn8 Sb7 has an approximate vibrational entropy of 0.19 kB /atom, which is slightly less than that of Zn4 Sb3 with 0.22 kB /atom. The grand canonical potential for each phase was calculated according to Equation 3.3 and the boundary lines of phase equilibria were determined by equality of two respective GCPs. Formation free energy curves, as a function of composition, are shown for Zn8 Sb7 and Zn4 Sb3 at three different temperatures in Figure 3.8. The set up blue curves at 650K correspond to the region of phase decomposition into ZnSb and Zn. At the intermediate temperature in orange, the formation free energy of Zn4 Sb3 has dropped below that of Zn and ZnSb. The convex hull connects the lowest energy phase at each composition and includes the respective common tangent construction 36 10 Formation Free Energy (meV/atom) 5 T=650K 0 −5 T=900K −10 −15 Zn8Sb7 −20 Zn4Sb3 0.52 0.53 0.54 0.55 0.56 Zn Atomic Fraction 0.57 T=1200K 0.58 0.59 Figure 3.8: Calculated formation free energy curves for Zn4 Sb3 and Zn8 Sb7 with respect to ZnSb and Zn at three temperatures. A convex hull is drawn connecting the free energy curves of the thermodynamically stable phase at a given temperature. At 650 K, Zn and ZnSb are in equilibrium over the entire composition range shown. At 900 K, Zn4 Sb3 has broken the convex hull and has a pure phase stability window on a narrow composition range. At 1200 K, Zn8 Sb7 has broken the convex hull and shows a favorable formation energy with respect to ZnSb and Zn4 Sb3 (connected by the black dashed line). between Zn4 Sb3 , ZnSb and Zn. At this temperature, there is no thermodynamic stability for Zn8 Sb7 even though it is more favorable than decomposition to ZnSb and Zn. At the highest temperature in red, we see that the Zn8 Sb7 free energy curve touches the convex hull indicating a region of phase stability. The black dashed line indicates the convex hull without the Zn8 Sb7 phase. We see that Zn8 Sb7 could have a temperature window of phase stability at high temperature, if there is sufficient entropy as in our estimation. We note that the calculated temperature scale does not correspond exactly to the physical Kelvin scale. Our temperatures likely over-predict the stabilization temperatures for both Zn4 Sb3 and Zn8 Sb7 . 37 With our plausible assumptions for configurational entropy we show that Zn8 Sb7 can be thermodynamic stable at high temperatures. We also consider two possible reasons as to why Zn8 Sb7 has only been observed experimentally in nanoparticulate form—lattice strain and surface energy. Ab initio calculations under the generalized gradient approximation typically overestimate the experimentally observed lattice parameter. This is indeed the case with ZnSb and Zn4 Sb3 , for which calculations predict, respectively, a 1.6% and 2.1% lattice parameter overestimation at 300K (including the calculated thermal expansion due to lattice vibrations). We would expect to see a linear trend with composition between ZnSb and Zn4 Sb3 . [52] Contrary to expectation, the calculated lattice parameter of Zn8 Sb7 under predicts the experimentally characterized lattice parameter by 1.0% at 300K as is shown in Table 3.1. This error is likely within the limits of accuracy for the experimentally characterized lattice constant, but the presence of lattice strain is possible, either due to a surface stress or other synthetic factors. Table 3.1: Bulk moduli, B, as fit to the equation of state E(V ) = 2VB0 (V − V0 )2 (for total energy, E, unit cell volume, V, and ground state volume, V0 ) and the room temperature unit cell volume (calculated and measured) for three Zn-Sb phases. Bulk Modulus Zn4 Sb3 Zn8 Sb7 ZnSb eV/Å3 0.267 0.240 0.248 GPa 42.7 38.5 39.7 Unit Cell Volume (Å3 ) DFT (GGA) 571 1577 410 Experiment 536 1625 391 Error 2.1% -1.0% 1.6% Surface effects could explain the more favorable formation of Zn8 Sb7 in nanoparticulate form. If the surface energy of Zn8 Sb7 was lower than the average surface energy of ZnSb and Zn4 Sb3 , then for a fixed number of atoms in a nanoparticle, we would expect Zn8 Sb7 to be more favored for small particle sizes. Calculating the surface energy for each of these structures would be computationally prohibitive, but we can consider values that might be plausible. As an example, if Zn8 Sb7 has a lower surface energy than ZnSb and Zn4 Sb3 by 0.05 eV/Å2 (0.78 J/m2 , which is on the order of differences reported in some elemental surface energies [53]) then for a 50nm par- 38 ticle, Zn8 Sb7 would have an energy benefit of 6meV/atom with respect to ZnSb and Zn4 Sb3 . At a particle size of 90nm, the surface energy difference would be negligible. Even if we assume that the surface energy is constant between all phases, we still observe that Zn8 Sb7 has a lower bulk modulus than both ZnSb and Zn4 Sb3 (Table 3.1). In a stressed surface condition, we expect the softer phase, Zn8 Sb7 , to have a more significant energy lowering to the surface contribution according to P 2 /2B, where P is pressure and B is the bulk modulus. Thus, Zn8 Sb7 would be slightly favored over ZnSb and Zn4 Sb3 . At the pressure required to adjust the calculated Zn8 Sb7 lattice parameter to the experimental value (-3.4 GPa), the energy correction is 3meV/atom in favor of Zn8 Sb7 . We suggest that through some combination of surface energy and lattice strain, Zn8 Sb7 might be more favorable under certain conditions. 3.4 Conclusions From our first principles investigation of Zn8 Sb7 we expect the theoretical bulk Zn8 Sb7 phase to exhibit good thermoelectric properties. We predict Zn8 Sb7 to be entropically stabilized at high temperature if our assumption of configurational entropy holds true. There are also other mechanisms, such as lattice strain or surface energy, that would favor the stability of Zn8 Sb7 with respect to ZnSb and Zn4 Sb3 . The stability of Zn8 Sb7 at high temperatures could partially explain the complexity observed in this part of the Zn-Sb phase diagram. 39 Chapter 4 Defect-Driven Properties in AZn2Sb2 4.1 Introduction Zintl compounds have enjoyed a rich history of scientific study since their initial investigation by Edward Zintl in the 1930s [51]. Intermetallic Zintl compounds are characterized by electropositive cations donating electrons to support the formation of anionic polyhedral units. Anion to anion bonds form to fully satisfy the s2 p6 configuration of valence electrons. For example, in ZnSb, each Zn atom donates two electrons to the valence band which is insufficient to satisfy the valence of each independent Sb atom. Instead, Sb atoms form bonded dimer units which satisfy the condition for valence balance. Likewise, in Zn4 Sb3 , Sb dimers as well as lone Sb are present to satisfy complete valence for a composition of Zn13 Sb10 , Zintl compounds have recently emerged as excellent thermoelectric materials due mainly to their intrinsically low lattice thermal conductivity. Good thermoelectric materials are generally degenerate semiconductors, with electronic properties optimized between an insulator (high Seebeck coefficient) and a metal (high electrical conductivity). Achieving the optimal balance between these properties requires controlling the electron (or hole) carrier concentration. Generally, the optimal carrier concentration, n, is on the order of 1019 to 1020 carriers/cm3 . By definition, Zintl compounds are valence precise, where—if all sites are assumed to be fully occupied—the Zintl condition (filled octets) is met. As such, in their pure form, Zintl compounds would be expected to have very few free charge carriers. How- 40 c a b Figure 4.1: (left) Layered crystal structure of AZn2 Sb2 , space group P-3m1. The cations are shown in green, Zn in blue, and Sb in orange. (right) The experimental Hall carrier concentration of AZn2 Sb2 compounds increased with increasing cation electronegativity. [4, 5, 6, 7, 8, 9, 10, 11, 12, 13] ever, there are many examples of Zintl compounds, which, despite being nominally valence precise, are inevitably found to be electron deficient, exhibiting large extrinsic carrier concentrations. This is particularly apparent in AZn2 Sb2 Zintl compounds (A=Ca, Sr, Eu, Yb) with the CaAl2 Si2 structure type (Figure 4.1 - left). These compounds exhibit persistently high extrinsic p-type carrier concentrations, as illustrated by the right panel of Figure 4.1. Across multiple studies, including both polycrystalline and single crystal samples, there is a clear trend between the electronegativity of A and the carrier concentration. Although generally less strong, similar trends are seen in other Zintl systems. For thermoelectric applications, high n in AZn2 Sb2 compounds is extremely useful, as they form naturally with nearly optimized electronic properties. Isoelectronic substitutions are sufficient to fine-tune the carrier concentration and electronic properties to optimize thermoelectric efficiency. In contrast, most classic Zintl compounds (Ca5 M2 Sb6 , Ca3 AlSb3 , Sr3 GaSb3 , Yb14 AlSb11 , Ba8 Ga16 Ge30 ), are intrinsic semiconducting materials with low carrier concentration (∼ 1018 cm−3 ), requiring intentional 41 doping with alio-valent elements to increase the carrier concentration to 1019 –1020 carriers/cm3 . Despite the advantage conferred by high n in AZn2 Sb2 compounds, the origin of high extrinsic carrier concentration in valence-precise semiconductors, and the dependence of n on electronegativity has not been satisfactorily explained. Here we show that the persistent p-type behavior observed in AZn2 Sb2 Zintl compounds and the dependence on electronegativity can be explained by considering defect concentrations. Using density functional calculations, we describe the thermodynamic behaviour of point defects in AZn2 Sb2 . The respective point defect concentrations contribute to a finite phase width and slightly off-stoichiometric materials. We show that the predicted deviation from the nominal stoichiometry can be used to explain the experimentally observed carrier concentrations in these materials. 4.2 Thermodynamics The formation and growth of a given crystal phase is governed by thermodynamics and kinetics, both of which are dependent on temperature and chemical composition. A phase may be thermodynamically stable over a range of composition, and the range varies with temperature. The limits of this composition range are determined by the energetics of configurational disorder and the equilibrium relationships with other nearby phases. Unless a sample is cooled extremely slowly during synthesis, the phase at room temperature may also exist outside of its thermodynamic range of stability, due to kinetic limitations. One form of configurational disorder that is present in all crystal phases is the point defect—a change in the crystal lattice occurring only at one lattice site. Even crystals with very precise stoichiometry will invariably contain point defects due to the disordering tendency of entropy. Point defects can significantly alter the electronic properties of a material, depending on the ionization of the defect and the defect energy relative to the valence and conduction bands. For a crystal phase that has a small concentration of defects, we can consider treat them in the dilute limit, under 42 the assumption that the defects are non-interacting and the crystal band structure is largely unaffected [54, 55]. In this limit, we can calculate the concentration of defects as a function of the single defect energy (∆H = ED –E0 ), the Fermi level (EF ) and the chemical potential of each atomic species (µ). In the case of a single cation vacancy, the expression becomes     ∆HD,q ED,q − E0 + qEF + µA xD,q ≈ x0 exp − = xA0 exp − kB T kB T (4.1) Intuitively, we would expect ED –E0 to correlate with the relative strength of bonds that may be broken or perturbed through the formation of a defect. The chemical potential range for the phase is constrained thermodynamically by equilibrium with nearby phases. In AZn2 Sb2 samples, ZnSb precipitates are frequently observed. The chemical potential, µ, can then be chosen for the condition of equilibrium between AZn2 Sb2 , ZnSb and Zn. The limits of defect concentrations occur at the limit of phase stability with these surrounding phases. A charge neutral condition must also be satisfied between charged defects and charged carriers: nh − ne + X qxD,q = 0 (4.2) D,q where nh and ne are the respective free carrier concentrations of holes and electrons: Z εV D(ε) [1 − f (ε; εF , T )] dε nh = −∞ Z ∞ ne = D(ε)f (ε; εF , T )dε εC where D(ε) is the density of states and f (ε; εF , T ) is the Fermi distribution function centered at chemical potential, εF . From these relationships, the concentrations of defects, holes and electrons can be calculated self-consistently at a given tempera- 43 ture, presuming complete thermodynamic equilibrium. However, we do not expect AZn2 Sb2 to be in equilibrium through the complete solid-state temperature range. At some temperature, defect formation and mobility become kinetically limited and the meta-stable concentrations at lower temperatures reflect the thermodynamic equilibrium of a higher temperature. 4.3 Results and Discussion The calculated electronic structures for AZn2 Sb2 compounds are shown along lines of high symmetry in Figure 4.2. The structure of the valence band is similar across the series—each consisting of three bands near the valence band maximum. These can be accurately described, analytically, as ellipsoidal pockets—each with anisotropic effective carrier mass. These bands are primarily of Sb p-orbital character, similar to the valence bands of other Sb-derived Zintl phases such as ZnSb [2, 14], Zn4 Sb3 [1], Ca5 M2 Sb6 [16] and Ca3 AlSb3 [16]. At zero temperature, the valence band states are fully occupied, (Fermi energy residing in the band gap) which would result in an intrinsic semiconducting material. This suggests that, for a perfect AZn2 Sb2 crystal, thermally activated intrinsic carriers should dominate transport. In contrast, experimental transport measurements show large concentrations of extrinsic p-type carriers, suggesting that in the measured samples, the Fermi level is in the valence band. The presence of band gaps in AZn2 Sb2 compounds is consistent with experimental results, which suggest band gaps ranging from 0.2 to 0.5 eV. The variation in the calculated bandgaps is within the range of experimental estimations, but we do not place a high value on their precision. However, we will show that the magnitude of the gap does not greatly influence the subsequent results, due to the highly degenerate p-type behavior. To understand the defect behavior of AZn2 Sb2 , we must understand the phase relationships with the other compounds in the A-Zn-Sb ternary phase diagram. Based on total energy calculations, AZn2 Sb2 is indeed stable at zero temperature with respect to its elemental constituents as well as to reported surrounding binary compounds. 44 0.5 CaZn2Sb2 0 −0.5 0.5 SrZn2Sb2 E − EV (eV) 0 −0.5 0.5 YbZn2Sb2 0 −0.5 G M 0.5 K G A L M G A L M EuZn2Sb2 0 −0.5 G M K Figure 4.2: Electronic structure of AZn2 Sb2 compounds along symmetric directions in reciprocal space, shows similar valence band structure. An anti-ferromagnetic supercell was used for EuZn2 Sb2 so the Brillouin Zone is different. 45 Here, we focus on the boundary between AZn2 Sb2 and ZnSb, since ZnSb precipitates are often observed after synthesis. The total energy of a defect in AZn2 Sb2 was calculated after removing or replacing a single atom from a large supercell. These defect energies (ED in Equation 4.1) were calculated for vacancies, as well as substitutional, and interstitial point defects. These include vacancies on the cation, Zn and Sb site (VA , VZn , VSb ) substitutions of Zn or Sb on the cation site (ZnA , SbA ), or cations on the Zn site (AZn ), as well as cation interstitials (Ai ). The defect formation enthalpy [54, 55], ∆H = ED –E0 + qEF + µ, is often reported as a function of EF . Figure 4.3 (top) shows ∆H for all calculated defects in CaZn2 Sb2 . The slope of the line corresponds to the charged state of the defect (for clarity, we only show the lowest energy charge segment for each defect). It is apparent from these results that the cation vacancy (VCa ) is the most energetically favorable point defect. This was found to be true for all AZn2 Sb2 compounds—the cation vacancy is always the predominant defect. It is thus only necessary to consider cation vacancies when comparing the defect energies across the AZn2 Sb2 series. The bottom panel of Figure 4.3 illustrates the dependence of vacancy energy on our choice of cation. As expected, the vacancies of the most electronegative cation, Yb, has the smallest energy cost. Vacancy formation becomes more costly as electronegativity decreases. An important implication of these results is that strong n-type behavior in AZn2 Sb2 compounds is thermodynamically unattainable. This arises from the fact that the negatively charged vacancy has a negative formation enthalpy for Fermi levels still in the middle of the band gap. The formation of a negatively charged defect effectively reduces the Fermi level to maintain charge balance. Such low formation energies of these cation vacancies pins the Fermi level closer to the valence band. This explains the persistent p-type character of AZn2 Sb2 compounds. Attempts to increase the Fermi level through n-type chemical doping would only encourage the formation of more cation vacancies. This is most dramatically the case in YbZn2 Sb2 , where the Fermi level is pinned very closely to the valence band maximum. At higher temperatures, the formation enthalpies of the vacancies are more easily overcome and a higher 46 VSb Cai 1.5 HDefect (eV) CaZn SbCa 1 ZnCa 0.5 VZn VCa 0 -0.2 -0.1 0 0.1 0.2 0.3 EF - EV (eV) 0.4 0.5 0.6 0.6 ∆HDefect (eV) 0.4 Sr 0.2 Eu 0 Ca −0.2 Yb −0.4 −0.2 0 0.2 EF − EV (eV) 0.4 0.6 Figure 4.3: (top) Defect formation enthalpy for various point defects in CaZn2 Sb2 as a function of electronic Fermi level. The slope of each segment (+1, 0, −1 or −2) corresponds to the charged state of the defect. For clarity, line segments are only shown for the lowest energy charge state of each defect type, but in reality extend beyond each intersection. In all AZn2 Sb2 compounds, cation vacancies, VA, were found to be the lowest energy defect. (bottom) With increasing cation electronegativity, the formation enthalpy, ∆Hdef ect for cation vacancies in AZn2 Sb2 compounds decreases. The colored dashed lines correspond to estimated conduction band minima for each compound. 47 Sr 900 Eu Ca Temperature (K) 800 A = Yb 700 600 500 ZnSb + A1-xZn2Sb2 A1-xZn2Sb2 400 300 ZnSb 0.04 0.03 0.02 x in A1-xZn2Sb2 0.01 0 Figure 4.4: Range of phase widths predicted for A1−x Zn2 Sb2 . The boundary line corresponds to the maximum vacancy concentration for A1−x Zn2 Sb2 compounds in thermodynamic equilibrium with ZnSb. vacancy concentration becomes stable—further reducing the Fermi level. The thermodynamically stable defect concentrations in AZn2 Sb2 determine the compositional phase width around the pure stoichiometry. This is illustrated for each AZn2 Sb2 compound in Figure 4.4 in the form of a pseudo-binary phase diagram between AZn2 Sb2 and ZnSb. The edge of each respective shaded region corresponds to the maximum concentration of vacancies in single-phase AZn2 Sb2 , assuming full equilibrium with ZnSb at a given temperature. Any composition and temperature range outside of the respective shaded region would result in ZnSb precipitates. As discussed above, it is unlikely that samples at room temperature are in thermodynamic equilibrium. Below a threshold temperature, diffusion kinetics become too slow to equilibrate the system. A reasonable guess for this temperature might be 800 K, but it likely varies depending on the cation species. Below this temperature, the vacancy concentration becomes effectively fixed and ZnSb precipitates—if any—would 48 have already been formed. The hole and electron carrier concentrations reflect the charge balance of the frozen vacancy concentration according to Equation 4.2. This results in the behavior of an extrinsically doped semiconductor due to intrinsically formed cation vacancies. Assuming that the defect concentration “freeze” at 800 K during cooling, the expected defect and hole carrier concentrations can be calculated as described above, and are listed in Table 4.1. These tabulated results for Hall carrier concentration show both qualitative and quantitative agreement with all experimental reports in these compounds. Table 4.1: Experimental Hall carrier concentration alongside calculated vacancy concentration, carrier concentration, and Hall carrier concentration for AZn2 Sb2 at 800 K. All values are in units of cm−3 . nh,exp [VA ] ncalc nhall,calc CaZn2 Sb2 SrZn2 Sb2 YbZn2 Sb2 EuZn2 Sb2 4.0 × 1019 1.5 × 1019 1.2 × 1020 2.5 × 1019 6.3 × 1019 3.6 × 1018 2.2 × 1020 2.0 × 1019 6.0 × 1019 5.2 × 1018 1.6 × 1020 2.5 × 1019 6.2 × 1019 1.0 × 1019 1.6 × 1020 2.9 × 1019 There is a clear trend in these results with the electronegativity of the cation in Figure 4.1. We have found that the energy of formation for cation defects in AZn2 Sb2 compounds scales inversely with electronegativity. This leads to highest vacancy concentrations in the Yb-based compounds. The trend in defect formation extends to carrier concentration to maintain charge balance. 4.4 Conclusion In materials such as AZn2 Sb2 , we show that a large concentration of thermodynamically stable cation vacancies leads to high extrinsic carrier concentrations. The stable defect level depends on our choice of A, and is consistent with experimentally observed carrier concentrations in these materials. These results demonstrate that point defects are the primary mechanism by which the covalency of the cation bond can influence carrier concentration in nominally valence-precise AZn2 Sb2 compounds. This 49 mechanism may be generally applicable, perhaps explaining similar trends seen in A14 MSb11 , and A2 MSb2 (A=2+ cation, M = 2+ or 3+ metal), and similar materials. 50 Appendix A Cluster Expansion of III–V Alloys The goal here is to create a lattice cluster expansion (CE) Hamiltonian that is capable of describing the entire phase space of a large alloy system. The CE would include scalar quantities, like formation energy, as well as tensorial properties, like lattice strain and elastic constants. The III–V semiconducting system was selected as a model system due to its order-dependent properties and epitaxial nature of many applications. The III–V alloys form in one of two primary lattices. The ground state lattice of the Phosphides and Arsenides—[Al,Ga,In][P,As]—is Zinc-blende [56]. The Nitride compounds—[Al,Ga,In]N—form a ground state in the wurtzite structure [56]. Many of the pure binary compounds form complete solid solutions with each other by substitution on the cation or anion site [57]. There are also many ordered heterostructures that form in the intermediate composition range due to the large lattice mismatch in some alloys. It is this problem that the cluster expansion of strain and elastic properties can help understand. A.1 Scalar Cluster Expansion for Energy The energy of a given configuration, E(σ), using the scalar cluster expansion, is represented as polynomial of occupation variables, σi . Here, the angled brackets represent the average product of all distinct clusters of a given motif, α: * E(σ) = X α mα Jα + Y i∈α0 σi 51 Figure A.1: Zinc-blende (left) and wurtzite (right) crystal lattices. Blue cations can be occupied by Al, Ga, and In. Orange anions can be occupied by N, P, and As. For example, in a large supercell, there could be hundreds of nearest neighbor interactions. The average is taken over the products of each distinct pair of occupation variables. For a given cluster motif, α, mα is the symmetric multiplicity and Jα is the effective cluster interaction (ECI). Select pair and triplet cluster motifs are shown in Figure A.2. A 2-component cluster expansion was fit or each of 18 binary alloy systems in [Al,Ga,In][N,P,As]. The convex hull is shown for each in the Zinc-blende lattice in Figure A.3. Red points are calculated using density functional (DFT) calculations and green points are predicted using the fitted CE. All but [Al,Ga]N show only the pure binary end-members as ground state configurations. In the case of [Al,Ga]N, the wurtzite structures have slightly lower energy, so the intermediate configuration is not an absolute ground state. From inspection of these convex hull diagrams, it becomes clear that some alloys—such as Al[N,As]—will have a strong tendency to phase separate, while others—such as [Al,Ga]P—have a much smaller energy cost for form a mixed configurations. Adding to these binary alloy data points, ternary and quaternary configurations were generated and fit with multi-component CEs—finally building up to the 6-component CE. 52 Figure A.2: Selected pair (left) and triplet (right) cluster motifs on the zinc-blende lattice. The multi-component CE is more complicated than its binary counterpart because the occupation variables cannot simply be ±1. The extra degrees of freedom from multiple site occupation require more than one cluster function for each cluster motif. To simplify the interpretation of the occupation variables, we used the following cluster functions: • Function 1— A: σi = − 12 , B: σi = 1, C: σi = − 21 – This can be interpreted as the effect of a B occupation relative to an random occupation of [A,B,C] atoms. • Function 2— A: σi = − 12 , B: σi = − 12 , C: σi = 1 – This can be interpreted as the effect of a C occupation relative to an random occupation of [A,B,C] atoms. • Redundant Function— A: σi = 1, B: σi = − 12 , C: σi = − 21 – Can be calculated from functions 1 and 2 for a symmetric interpretation of the average occupation variables. 53 Al(N,P) 0.7 Al(P,As) 0.02 predicted known str known gs 0.9 0.8 0.016 0.7 0.014 0.4 0.6 0.012 0.5 energy 0.3 0.4 0.3 0.006 0.1 0.2 0.004 -0.1 0.1 0 0.2 0.4 0.6 concentration 0.8 1 0.002 0 0 -0.1 -0.002 0 0.2 Ga(N,P) 0.6 0.4 0.6 concentration 0.8 1 0.8 0.02 energy energy energy 0.3 0.2 0.4 0.6 concentration 0.8 1 0.002 0 0.2 Ga(N,P) 0.35 0.2 0.15 0.8 -0.05 1 0.2 0 0 0.05 0.8 0 1 0.002 0.0015 0.2 0.0005 0.01 0 0 0.8 -0.01 1 0.003 0.02 0.001 0.01 0.4 0.6 concentration 0.08 0.8 1 0 1 predicted known str known gs 0.07 0.04 0.03 0.02 0.01 0 0.2 0.4 0.6 concentration 0.8 0 1 0 0.2 0.4 0.6 concentration 0.8 1 (Ga,In)As 0.07 predicted known str known gs predicted known str known gs 0.06 0.05 0.04 0.03 0.002 0.8 (Ga,In)P energy 0.04 energy 0.05 0.4 0.6 concentration 0.05 0.06 0.004 0.2 0.2 (Al,In)As 0.005 0 0 0.06 0.07 predicted known str known gs 0.006 -0.02 1 predicted known str known gs (Al,Ga)As 0.007 0.8 0.03 0.02 0.4 0.6 concentration 0.4 0.6 concentration 0.04 0.001 0.2 0.06 0 0 energy energy 0.05 0 0.08 0.02 0.07 0.06 -0.0005 predicted known str known gs 0.04 0.08 0.003 1 0.1 0.09 0.0025 0.8 (Ga,In)N (Al,In)P predicted known str known gs 0.004 0.0035 0.4 0.6 concentration 0.12 (Al,Ga)P 0.0045 0.2 0.14 energy energy 0.1 0.4 0.6 concentration 0 0.16 0.2 0.2 0.006 1 0.15 0.005 0 0.8 predicted known str known gs 0.25 0.01 0 0.008 (Al,In)N 0.015 -0.005 0.4 0.6 concentration 0.3 predicted known str known gs 0.02 predicted known str known gs 0.002 0 (Al,Ga)N 0.025 1 0.004 0 0.4 0.6 concentration 0.8 0.01 0.25 0.05 0 0.4 0.6 concentration 0.012 0.1 0.05 0.2 0.2 In(P,As) energy energy energy 0.15 0 0 0.014 0.3 0.1 energy 0 1 predicted known str known gs 0.35 0.2 energy 0.8 Ga(N,As) 0.4 0.25 -0.05 0.4 0.6 concentration 0.45 predicted known str known gs 0.3 0.01 0.008 0.004 0 -0.1 0.012 0.006 0.1 0 1 0.014 0.4 0.2 0.1 0.8 predicted known str known gs 0.016 0.5 0.2 0.4 0.6 concentration 0.018 0.6 0.3 0 0.2 Ga(P,As) predicted known str known gs 0.7 0.4 -0.1 0 Ga(N,As) predicted known str known gs 0.5 energy 0.01 0.008 0.2 0 predicted known str known gs 0.018 0.5 energy energy Al(N,As) 1 predicted known str known gs 0.6 0.03 0.02 0.01 0 0 0.2 0.4 0.6 concentration 0.8 1 -0.01 0 0.2 0.4 0.6 concentration 0.8 1 Figure A.3: Pseudo-binary convex hulls for all [Al,Ga,In][N,P,As] combinations on the zinc-blende lattice. Green points are predicted configurational energies based on a 2-component cluster expansion fit to red points. 54 The 6-component scalar CE because increasingly difficult to fit. To improve the fit, two schemes were utilized. First, a weighted least-squares fit was used to reduce the residual error of the lowest energy structures. The weight of each data point was
constructed from a Boltzmann factor—exp − EkTF —according to the formation energy of the configuration with respect to the pure binary end-members. The residual error of the fit is represented in Figure A.4. There is generally good agreement and a cross-validation (CV) score of 25 meV/atom. The inset shows an enlarged view of the configurations with very low formation energy. While the ground state configurations are correctly represented, there is still significant scatter for the low energy configurations that could adversely affect the Monte Carlo simulations. As a potential improvement, a Bayesian fitting procedure was used [58]. This approach has a much different view of the data being fit. It assumes that the configurational energy data is ’precisely known‘, rather than a frequentist view that the data has a random sampling error. This is appropriate with the use of DFT, which may have a systematic error for certain quantities, but should be very precise in predicting the total energy. With a Bayesian approach, there can be more fitting terms than data points. Using a prior assumption, the most likely set of fitting parameters can be determined. The prior used for fitting the cluster ECIs is a Gaussian function about zero: 0 p (J) = Y α ωα ωα = Abnα 1 √   Jα2 exp − 2 2ωα 2π Y  rij −p a {i,j}⊂α where the Gaussian width, ωα , of each ECI is decaying function of cluster interaction distance, rij and the number of sites in the cluster, nα . The adjustable parameters, b and p can be selected according to physical knowledge of the system and to optimize the predictive power according to the CV score. This fitting method was successfully applied to the ternary [Al,Ga,In]N wurtzite system. The known DFT energies are nearly perfectly matched and the fit has a CV score of 0.4 meV/atom. The ECIs are shown in Figure A.5 on a logarithmic 55 Predicted Formation Energy vs DFT Formation Energy 0.01 0.7 Cluster Expansion Formation Energy (eV) 0.0075 0.005 0.6 0.0025 0.5 0 0 0.0025 0.005 0.0075 0.01 0.4 0.3 0.2 0.1 0 0 0.1 0.2 0.3 0.4 0.5 0.6 0.7 DFT Formation Energy (eV) Figure A.4: Cluster expansion fit error for 6-component III–V zinc-blende compounds. Formation energy is with respect to the 9 binary end-members. In the least-squares fit, points are weighted with a Boltzmann factor according to the formation energy of each configuration. The inset shows an enlarged view of the configurations with very low formation energy. 56 3 12 3 12 3 8 3 5 3 5 0 10 −2 ECI ECI (eV) (eV) 10 −4 10 −4 -10 −2 -10 0 -10 3 Pairs (Distance,Ang) Triplets 8 Quadruplets Figure A.5: Cluster expansion fit for pseudo-ternary III–N wurtzite compounds using a Bayesian fitting approach. ECIs are reported for 2000 pair, triplet, and quadruplet clusters. Lines correspond to the Bayesian prior assumption of ECI magnitude. scale, along with lines showing the decay of the Bayesian prior with effective cluster distance. The prior parameters, b = 0.25 and p = −4, yielded the best CV score, though many other parameter combinations gave similar fitting results. Future study of this fitting method should include a test of robustness of the fitted ECIs to the variation in prior parameters. A.2 Tensorial Cluster Expansion of Strain and Elastic Constants To extend the scalar cluster expansion to tensorial properties, a tensor basis set is included in the expansion [25]: 57 Table A.1: Basis set for 2-dimensional strain tensor for the zinc-blende lattice up to 2nd -nearest neighbors. Cluster Symmetrically Allowed Tensor Basis   1 0 0 Null & Point Cluster β = √13 0 1 0   0 01  1 0 0 0 1 1 1st NN β = √13 0 1 0 , √16 1 0 −1  0 0 1  1 −1 0    1 0 0 0 0 0 0 0 0 0 1 1 2nd NN β = 0 0 0 , √12 0 1 0 , √12 0 0 1 , 21 1 0 0 0 0 0 0 0 1 0 1 0 1 0 0 * Q(σ) = X X α β∈C(α) mαβ Jαβ β0 + Y i∈α0 σi (α,β) In this case, each cluster, α, is coupled to a set of tensor bases, β, restricted by the symmetry of the cluster motif. Table A.1 lists the set of allowed 2-dimensional tensor bases for the first few cluster motifs on the zinc-blende lattice. Table A.2 shows the same for 4-dimensional tensors. Significant work remains to determine the best way to fit the tensorial CE to the DFT data. A simple least-squares fit was performed for strain and elastic constants for the full 6-component III–V system, using only pair clusters up to 2nd nearest neighbors. Figure A.6 shows the fitted ECI results for strain—visualizing the directional tensor coupling of each cluster mapped onto the unit sphere. Figure A.7 shows the coupling interactions for the elastic constants. Further work is needed to gauge the quality and predictive ability of these CE fits. As is, these results can still provide valuable insight into the coupling of order and directional properties. One of the ultimate goals of the tensorial cluster expansion is to be able to perform lattice Monte Carlo simulations under stressed and strained conditions. The grand canonical ensemble would include additional terms to account for these applied conditions: 58 Table A.2: Basis set for 4-dimensional elastic tensor for the zinc-blende lattice up to 2nd -nearest neighbors. Cluster Symmetrically Allowed Tensor Basis 1 0 0 0 0 0 0 0 0 Null & Point Cluster 1st NN 2nd NN 0 0 0 0 0 0 0 0 0 0 0 0 0 1 0 0 0 1 000000000 010000000 000100000  00 00 00 00 00 00 00 00 00   00 00 10 01 00 00 00 00 00   00 01 00 00 00 00 10 00 00        1 1 1 β = √ 0 0 0 0 1 0 0 0 0, √ 0 0 0 0 0 0 0 0 0, √ 1 0 0 0 0 0 0 0 1 3 0 0 0 0 0 0 0 0 0 6 0 0 0 0 0 1 0 0 0 12  0 0 0 0 0 0 0 1 0  000000000 000000100 001000000 000000000 000000010 000001000 000000001 0 0 0 00 0 0 0 0 100010000 0 0 0 0 0 0 0 0 0 1 0 0 0 0 0 0 0 0 0 0 1 0 0 0 0 −1 0 000000000 0 1 0 0 0 −1 0 0 0   00 00 00 00 00 00 00 00 00   0 0 0 0 0 −1 1 0 0    0 0 0 0 0 0 0 0 0, β = √1  0 0 0 0 1 0 0 0 0  , √1  3 0 0 0 0 0 0 0 0 0 12  0 0 −1 −1 0 0 0 0 0   0 0 0 1 0 0 0 −1 0  000000000 000000000 0 −1 0 0 0 0 −1 0 0 000000001  0 0 0 00 00 00 01 0 0 0 1 0 0 0 0 0 0 0 0 0 0 0 0 0 0 0 0 0 −1 1 0 −1 010000000  0 0 0 1 −1 0 0 −1 0   00 00 10 01 00 00 00 00 00   0 0 1 0 0 0 0 −1 −1    1 1  0 0 −1 0 0 0 −1 0 0  , √ 0 0 0 0 0 0 0 0 0, √ 6 0 0 0 0 0 1 0 0 0 24  1 −1 0 0 0 0 −1 0 0   0 1 0 0 −1 −1 0 0 0  000000100 000000010 1 0 −1 −1 0 0 0 0 0 000000000 0 −1 0 0 0 0 0  0 10 1−1 0 0 0 0 1 0 0 0 1 1 0 0 1 0 0  1000 1 0 0 0 0 000100000 1 0 0 0 0 0 0 0 1   00 01 00 00 00 00 10 00 00  1 0 0 0 1 0 0 0 0    √1  1 0 0 0 0 0 0 0 1  , √1  0 1 0 1 0 −1 0 −1 0  12  0 0 0 0 0 0 0 1 0  24  0 0 0 0 −1 0 0 0 −1  1 0 0 0 0 0 0 0 1  001000000 0 0 0 0 −1 0 0 0 −1 000001000 100010000 1 −1 0 1 0 0 0 0 0 0 0 0 0 00 01 00 0 00 0−1 0 0 0 000000000 000000000  00 00 00 00 00 00 00 00 00   00 00 00 00 00 00 00 00 00      β =  0 0 0 0 0 0 0 0 0  , 21  0 0 0 0 0 0 0 0 1  , 0 0 0 0 0 0 0 0 0 0 0 0 0 0 0 0 1 0 000000000 000000000 000000000 000001000  0 00 00 00 00 00 00 00 00 0  0 00 00 00 00 10 00 00 00 0 010000000 001000000 0 0 1 0 0 0 0 0 0 0 1 0 0 0 0 0 0 0 1 0 0 0 1 0 0 0 0 0 1 0 0 0 0 0 0 1 0 0 0 0 0 0 0 0 0 0 0, 0 0 0 0 0 0 0 0 0, 2  2 0 0 0 0 0 0 0 0 0 0 0 0 0 0 0 0 0 0 000000100 000100000 000000000 000000000 000000000   00 00 00 00 00 00 00 00 00   000000000 000000000 000000000 0 0 0 0 0 0 0 0 0 0 0 0 0 0 0 0 0 0  00 00 00 00 00 00 00 00 00      1 1 √ 0 0 0 0 1 0 0 0 0, √ 0 0 0 0 0 0 0 0 0, 2 0 0 0 0 0 0 0 0 0 2 0 0 0 0 0 1 0 0 0 000000000 000000000 000000000 000000010  00 00 00 00 00 00 00 00 10   00 00 00 00 01 00 00 00 01  000000000 000100000  00 00 00 00 00 00 00 00 00   00 01 00 00 00 00 10 00 00      √1  0 0 0 0 0 1 0 1 0  , √1  1 0 0 0 0 0 0 0 0  , 8 0 0 0 0 1 0 0 0 1 8 0 0 0 0 0 0 0 0 0 000000000 001000000 000010001 000000000  00 00 00 00 00 11 00 11 00   10 01 01 01 00 00 01 00 00  000000100 100000000  00 00 01 10 00 00 00 00 00   11 00 00 00 00 00 00 00 00      1 1 √ 0 0 0 0 0 0 0 0 0, √ 0 0 0 0 0 0 0 0 0, 8 1 0 0 0 0 0 0 0 0 8 0 0 0 0 0 0 0 0 0 010000000 100000000 100000000 000000000 000000000 0 0 0 0 00000 0 0 0 0 0 0 0 0 0 0 0 0 0 0 0 0 0 0 0 0 0 0 0 0 0 0 0 0 0 0 0 0 1 0 −1 1 000010000 000001021 0 0 0 0 1 −1 0 1 0    00 00 00 00 01 00 00 00 10   00 00 00 00 10 22 00 11 01   0 0 0 0 0 −1 0 1 1      1 1 1 0 0 1 0 0 0 1 0 0 √ 0 1 0 1 0 0 0 0 0, √ 0 0 1 0 0 0 1 0 0, √  8 0 0 0 0 0 0 0 0 0 32  0 1 2 2 0 0 1 0 0  24  0 1 −1 −1 0 0 1 0 0   0 0 0 0 1 1 0 −1 0  000000001 000011020 000000000 021100100 0 −1 1 1 0 0 −1 0 0 001000100 010100000 0 1 0 1 0 0 0 0 0 59 Figure A.6: Coupled tensor-cluster interaction terms for lattice strain mapped onto the unit sphere. Lighter colors represent strong positive coupling between a particular cluster and the direction of strain. Figure A.7: Coupled tensor-cluster interaction terms for elastic constants mapped onto the unit sphere. Lighter colors represent strong positive coupling between a given cluster and the directional elastic constants. 60    P P Cs,ijkl ij kl Es − µ · Ns + σs,ij ij + 21 i,j i,j,k,l  X  Φ(µ, T, ) = −kT ln  exp −  kT s where  is the 2-dimensional strain tensor, C, is the 4-dimensional elastic tensor, and σ is the 2-dimensional stress tensor—all of which could be calculated using the tensor CE. In this way, the large phase space of an alloy system can be explored in various epitaxial conditions or accounting for strain fields. This effort still holds promise to be an effective tool, but additional work is required to refine the fitting approaches to find an accurate and robust CE fit. 61 Bibliography [1] G. S. Pomrehn, E. S. Toberer, G. J. Snyder, A. van de Walle, Entropic stabilization and retrograde solubility in Zn4 Sb3 , Physical Review B 83, 094106 (2011). [2] G. S. Pomrehn, E. S. Toberer, G. J. Snyder, A. van de Walle, Predicted electronic and thermodynamic properties of a newly discovered Zn8 Sb7 phase, Journal of the American Chemical Society 133, 11255-11261 (2011). [3] E. S. Toberer, P. Rauwel, S. Gariel, J. Tafto, G. Jeffrey Snyder, Composition and the thermoelectric performance of β-Zn4 Sb3 , Journal of Materials Chemistry (2010). [4] E. S. Toberer, A. F. May, B. C. Melot, E. Flage-Larsen, G. J. Snyder, Electronic structure and transport in thermoelectric compounds AZn2 Sb2 (A = Sr, Ca, Yb, Eu), Dalton Transactions 39, 1046-1054 (2010). [5] A. F. May, et al., Properties of single crystalline AZn2 Sb2 (A = Ca, Eu, Yb), Journal of Applied Physics 111 (2012). [6] K. Guo, et al., Enhanced Thermoelectric Figure of Merit of Zintl Phase YbCd2−x Mnx Sb2 by Chemical Substitution, European Journal of Inorganic Chemistry pp. 4043–4048 (2011). [7] H. Zhang, et al., Thermoelectric properties of Ybx Eu1−x Cd2 Sb2 , Journal of Chemical Physics 133 (2010). [8] H. Zhang, et al., Thermoelectric Properties of Polycrystalline SrZn2 Sb2 Prepared by Spark Plasma Sintering, Journal of Electronic Materials 39, 1772-1776 (2010). 62 28th International Conference on Thermoelectrics/7th European Conference on Thermoelectrics, Freiburg, Germany, JUL 26-30, 2009. [9] Q.-G. Cao, et al., Zintl phase Yb1−x CaxCd2 Sb2 with tunable thermoelectric properties induced by cation substitution, Journal of Applied Physics 107 (2010). [10] H. Zhang, et al., Thermoelectric properties of EuZn2−x Cdx Sb2 , Dalton Transactions 39, 1101-1104 (2010). [11] H. Zhang, et al., Synthesis and properties of CaCd2 Sb2 and EuCd2 Sb2 , Intermetallics 18, 193-198 (2010). [12] X.-J. Wang, et al., Synthesis and high thermoelectric efficiency of Zintl phase YbCd2−x Znx Sb2 , Applied Physics Letters 94 (2009). [13] H. Zhang, et al., A new type of thermoelectric material, EuZn2 Sb2 , Journal of Chemical Physics 129 (2008). [14] P. H. M. Böttger, G. S. Pomrehn, G. J. Snyder, T. G. Finstad, Doping of ptype ZnSb: Single parabolic band model and impurity band conduction, Physica Status Solidi (a) 208, 2753–2759 (2011). [15] Y. Takagiwa, Y. Pei, G. Pomrehn, G. J. Snyder, Dopants effect on the band structure of PbTe thermoelectric material, Applied Physics Letters 101, 092102 (2012). [16] A. Zevalkink, et al., Influence of the Triel Elements (M = Al, Ga, In) on the Transport Properties of Ca5 M2 Sb6 Zintl Compounds, Chemistry of Materials 24, 2091-2098 (2012). [17] A. Zevalkink, et al., Thermoelectric properties of Sr3 GaSb3 — a chain-forming Zintl compound, Energy and Environmental Science pp. – (2012). [18] W. Zeier, et al., Phonon scattering through a local anisotropic structural disorder in the thermoelectric solid solution cu2zn1-xfexgese4., Journal of the American Chemical Society (2012). 63 [19] A. van de Walle, G. Ceder, The effect of lattice vibrations on substitutional alloy thermodynamics, Reviews of Modern Physics 74, 11 (2002). [20] B. Fultz, Vibrational thermodynamics of materials, Progress in Materials Science 55, 247 - 352 (2010). [21] A. van de Walle, M. Asta, Self-driven lattice-model monte carlo simulations of alloy thermodynamic properties and phase diagrams, Modelling and Simulation in Materials Science and Engineering 10, 521 (2002). [22] F. Ducastelle, Order and phase stability in alloys, Cohesion and structure (NorthHolland, 1991). [23] E. Ising, Beitrag zur theorie des ferromagnetismus, Zeitschrift für Physik 31, 253-258 (1925). [24] J. Sanchez, F. Ducastelle, D. Gratias, Generalized cluster description of multicomponent systems, Physica A: Statistical Mechanics and its Applications 128, 334 - 350 (1984). [25] A. van de Walle, A complete representation of structure-property relationships in crystals, Nature Materials 7, 455 (2008). [26] A. van de Walle, et al., Symmetry-adapted bases for the parametrization of anisotropic properties, arXiv preprint arXiv:1301.0168 (2013). [27] T. Caillat, J.-P. Fleurial, A. Borshchevsky, Preparation and thermoelectric properties of semiconducting Zn4 Sb3 , Journal of Physics and Chemistry of Solids 58, 1119 - 1125 (1997). [28] G. Nakamoto, et al., Homogeneity range of epsilon-phase Zn4 Sb3 , Journal of Alloys and Compounds 432, 116-121 (2007). [29] G. Snyder, M. Christensen, E. Nishibori, T. Caillat, B. Iversen, Disordered zinc in Zn4 Sb3 with phonon-glass and electron-crystal thermoelectric properties, Nature Materials 3, 458-463 (2004). 64 [30] S.-G. Kim, I. I. Mazin, D. J. Singh, First-principles study of zn-sb thermoelectrics, Physical Review B 57, 6199-6203 (1998). [31] F. Cargnoni, et al., Interstitial Zn atoms do the trick in thermoelectric zinc antimonide, Zn4 Sb3 : A combined maximum entropy method x-ray electron density and ab initio electronic structure study, Chemistry - A European Journal 10, 3861-3870 (2004). [32] E. S. Toberer, et al., Local structure of interstitial Zn in β-Zn4 Sb3 , Physica Status Solidi - Rapid Research Letters 1, 253-255 (2007). [33] A. S. Mikhaylushkin, J. Nylen, U. Haussermann, Structure and bonding of zinc antimonides: Complex frameworks and narrow band gaps, Chemistry - A European Journal 11, 4912-4920 (2005). [34] J. Nylen, M. Andersson, S. Lidin, U. Haussermann, The structure of α-Zn4 Sb3 : Ordering of the phonon-glass thermoelectric material β-Zn4 Sb3 , Journal of the American Chemical Society 126, 16306-16307 (2004). PMID: 15600316. [35] J. Nylen, et al., Low-temperature structural transitions in the phonon-glass thermoelectric material β-Zn4 Sb3 : Ordering of Zn interstitials and defects, Chemistry of Materials 19, 834-838 (2007). [36] H. Yin, et al., Thermal stability of thermoelectric Zn4 Sb3 , Journal of Electronic Materials 39, 1957-1959 (2010). 10.1007/s11664-009-1053-3. [37] A. van de Walle, M. Asta, G. Ceder, The Alloy Theoretic Automated Toolkit: A user guide, Calphad Journal 26, 539 (2002). [38] A. van de Walle, Multicomponent multisublattice alloys, nonconfigurational entropy and other additions to the Alloy Theoretic Automated Toolkit, Calphad Journal 33, 266 (2009). [39] A. van de Walle, G. Ceder, Automating first-principles phase diagram calculations, Journal of Phase Equilibria 23, 348 (2002). 65 [40] A. N. Qiu, L. T. Zhang, J. S. Wu, Crystal structure, electronic structure, and thermoelectric properties of β − zn4 sb3 from first principles, Physical Review B 81, 035203 (2010). [41] Ø. Prytz, et al., Nanoscale inclusions in the phonon glass thermoelectric material Zn4 Sb3 , Philosophical Magazine Letters 89, 362 - 369 (2009). [42] T. Ikeda, V. A. Ravi, G. J. Snyder, Formation of sb2te3 widmanstätten precipitates in thermoelectric pbte, Acta Materialia 57, 666 - 672 (2009). [43] C. S. Birkel, et al., Solution synthesis of a new thermoelectric zn1+xsb nanophase and its structure determination using automated electron diffraction tomography, Journal of the American Chemical Society 132, 9881-9889 (2010). [44] V. Izard, M. Record, J. Tedenac, Mechanical alloying of a new promising thermoelectric material, sb3zn4, Journal of Alloys and Compounds 345, 257-264 (2002). [45] J.-B. Li, M.-C. Record, J.-C. Tedenac, A thermodynamic assessment of the sb-zn system, Journal of Alloys and Compounds 438, 171 - 177 (2007). [46] X. Liu, C. Wang, I. Ohnuma, R. Kainuma, K. Ishida, Thermodynamic assessment of the phase diagrams of the cu-sb and sb-zn systems, Journal of Phase Equilibria 21, 432-442 (2000). 10.1361/105497100770339608. [47] J. P. Perdew, K. Burke, M. Ernzerhof, Generalized gradient approximation made simple, Physical Review Letters 77, 3865 (1996). [48] P. E. Blöchl, Projector augmented-wave method, Physical Review B 50, 17953 (1994). [49] G. Kresse, J. Furthmüller, Efficiency of ab-initio total energy calculations for metals and semiconductors using a plane-wave basis set, Computational Materials Science 6, 15 (1996). 66 [50] G. Madsen, D. Singh, Boltztrap. a code for calculating bandstructure dependent quantities, Computer Physics Communications 175 (2006). [51] S. M. Kauzlarich, Chemistry, Structure, and Bonding of Zintl Phases and Ions (Wiley-VCH, 1996). [52] A. van de Walle, G. Ceder, Correcting overbinding in LDA calculations, Physical Review B 59, 14992 (1999). [53] H. L. Skriver, N. M. Rosengaard, Surface energy and work function of elemental metals, Physical Review B 46, 7157–7168 (1992). [54] S. Lany, A. Zunger, Accurate prediction of defect properties in density functional supercell calculations, Modelling and Simulation in Materials Science and Engineering 17 (2009). [55] S. Lany, A. Zunger, Assessment of correction methods for the band-gap problem and for finite-size effects in supercell defect calculations: Case studies for ZnO and GaAs, Physical Review B 78 (2008). [56] C. Yeh, Z. Lu, S. Froyen, A. Zunger, Zinc-blende–Wurtzite Polytypism in Semiconductors, Physical Review B 46, 10086-10097 (1992). [57] I. Vurgaftman, J. Meyer, L. Ram-Mohan, Band parameters for III-V compound semiconductors and their alloys, Journal of Applied Physics 89, 5815-5875 (2001). [58] E. Cockayne, A. van de Walle, Building effective models from scarce but accurate data: Application to an alloy cluster expansion model, Physical Review B 81, 012104 (2010).