Two-Layer High-Throughput: Effective Mass Calculations Including Warping

Andrew Supk, Nichols A. Mecholsky, Mrco Buongiorno Nrdelli, Stefno Curtrolo,Mrco Fornri,d,*

a Department of Physics & Science of Advanced Materials Program, Central Michigan University, Mount Pleasant, MI 48859, USA

b Department of Physics & Vitreous State Laboratory, The Catholic University of America, Washington, DC 20064, USA

c Department of Physics & Department of Chemistry, University of North Texas, Denton, TX 76203, USA

d Center for Autonomous Materials Design, Duke University, Durham, NC 27708, USA

e Department of Mechanical Engineering and Materials Science, Duke University, Durham, NC 27708, USA

Keywords:High-throughput Electronic structure Band warping Effective mass

ABSTRACT In this paper, we perform two-layer high-throughput calculations. In the first layer, which involves changing the crystal structure and/or chemical composition, we analyze selected III–V semiconductors,filled and unfilled skutterudites,as well as rock salt and layered chalcogenides.The second layer searches the full Brillouin zone(BZ)for critical points within 1.5 eV(1 eV=1.602176×10-19 J)of the Fermi level and characterizes those points by computing the effective masses.We introduce several methods to compute the effective masses from first principles and compare them to each other. Our approach also includes the calculation of the density-of-states effective masses for warped critical points, where traditional approaches fail to give consistent results due to an underlying non-analytic behavior of the critical point. We demonstrate the need to consider the band structure in its full complexity and the value of complementary approaches to compute the effective masses. We also provide computational evidence that warping occurs only in the presence of degeneracies.

1. Introduction

The electronic band structure plays a fundamental role in our understanding of the origins of the physical properties of materials and in assessing paths for optimization and chemical substitutions.The dispersion relation of the solutions of the many-body electronic Schr?dinger equation provides quantitative information that is essential for understanding most of the functionalities that a material may exhibit.The band structure,En(k ),is a mapping from R3→ RN(where N is the number of relevant bands,n is the band index,and k is a vector indicating the crystalline momentum)and is usually represented by considering only two-dimensional (2D)band plots. Because any given band is a function whose domain is the three-dimensional (3D) Brillouin zone (BZ), k ∈ BZ ? R3,the graph of a given band is embedded in four dimensions. Band structure plots along high-symmetry lines are a useful tool for evaluating a material at a glance. However, the 2D representation hides the full complexity of the electronic spectrum by ignoring large sections of the BZ.Conventions,such as those in Ref.[1],provide common ground for band plots, but 2D band representations are intrinsically limited.

Critical points, where ?En(k )/?k = 0 (with ?En(k )/?k indicating the gradient of En(k ) with respect to k), are an important feature of the electronic band structure. At these points, the density of states (DOS; also represented by the function of energy E,D(E)) is large (or diverges; see Ref. [2]):

with the unit cell volume V and an infinitesimal element of the constant energy surface dS.Critical points are key in evaluating a material’s physical properties and fully characterizing the material. For example,for a 2D material,it is well known that saddle points lead to logarithmic singularities [3], and that maxima and minima in three dimensions lead to square root dependencies of the DOS and,in general,to van Hove singularities or non-smooth points[4].Thus, the identification of all critical points is an important goal.

A related concept associated with the local properties of the electronic band structure involves the effective mass tensor M*,which is—assuming Taylor expansibility near k0—a second-rank tensor in 3D with components:

where n is the band index, ˉh is the reduced Planck constant, meis the mass of electron, and the subscripts i, j are used to label the Cartesian components of the tensor M*or of the vector k.The reciprocal of the effective mass tensor is associated with the curvature of the energy dispersion, En(k ), and is a critical descriptor when discussing electronic transport and optical properties; its evaluation must be done by considering the full mathematical complexity of the band structure. In addition, the possible presence of nonanalytic points—such as points where the Hessian is not a symmetric matrix—leads to warped critical points,which play a role in several situations[5,6] but are difficult to identify with 2D band plots.By determining high-fidelity effective mass tensors at critical points, it is possible to formulate analytic models of band structures. Such models are important in many considerations in solidstate physics and electronic engineering; for example, they can be used as a starting point for Monte Carlo transport simulations [7]or for the multi-scale modeling [8] of electronic devices, batteries,and thermoelectric energy converters. Analytic band structures are used in applications such as modeling scattering rates where derivatives and integrals of the bands are necessary [7]. From an experimental standpoint,the DOS effective mass is a common property of the band structure that may be measured through cyclotron resonance [9–18], the four coefficients method [19], Shubnikov–de Haas oscillations [16,20–23], magnetophonon resonance[15,16,24,25], time-of-flight drift velocity [7,26–29], optical transmission and reflection[30–32], and infrared reflection and Faraday rotation [33–35]. Being able to compare the different measurements of the effective mass obtained with those indirect methods and reconcile such results with electronic structure calculations is a major goal.

This paper introduces a two-layer high-throughput (HT)methodology with the aim of partially addressing the misalignement between theory and experiment in the current literature.The first layer is a conventional chemical substitution/structural variation(e.g.,see Refs.[36–42]),and the second is a careful exploration of En(k )to identify the nature of the critical points and compute the effective mass tensors. The conventional HT considers an array of prototypical materials: III–V semiconductors, filled and unfilled skutterudites, as well as rock salt and layered chalcogenides. For each material, we search the full BZ for critical points and determine the effective mass at those points using several different methodologies for verification [43]. These methodologies also characterize the warping of bands and correctly compute the DOS effective mass at those warped points [6]. To our knowledge,this is the first HT calculation of effective masses in a large portion of the BZ.

In Section 2, we illustrate the details of the computations.Section 3 presents selected results(a large part of the data is included in Appendix A),and Section 4 discusses the impact of this work.

2. Methods

The prototypical materials selected for this work are as follows:III–V semiconductors(AlSb and AlP with the zincblende structure);rock salt (PbTe, GeTe, SnTe, PbS, GeS, SnS, PbSe, GeSe, and SnSe)and layered chalcogenides (Bi2Te3, Bi2Te2Se, Bi2Se2Te, Bi2Se3,Bi2Te2S, and Bi2Se2S); and pristine and fully filled cobalt antimonide (specifically CoSb3, CaCo4Sb12, and BaCo4Sb12). For each material, the workflow starts by generating a projected atomic orbital tight-binding (PAO-TB) Hamiltonian, which is exploited to interpolate the band structure efficiently and precisely [44]. We use Quantum Espresso (Quantum ESPRESSO Foundation, UK)[45,46] to calculate the electronic structure in the AFLOWπ HT computational framework [47].

The AFLOWπ’s workflow (Fig. 1) drives the calculation of the Hubbard U correction within the ACBN0 [48–51] scheme (see Tables S1,S4,S14,and S18 in Appendix A),optimizes the structure of the unit cell,and generates the PAO-TB Hamiltonian.The wavefunction and charge kinetic energy cutoffs were 150 and 600 Ry(1 Ry = 2.179872 × 10-18J), respectively, and a Monkhorst–Pack k-point grid with a density of about 0.01 ?-1has been used. The choice of pseudopotentials was driven by a need to maximize the number of well-projected bands in the PAO-TB model;for this purpose, Perdew–Burke–Ernzerhof (PBE)-projector augmented wave(PAW) pseudopotentials generated from the PSlibrary [52] were modified to have an extended basis [44]. The computations included spin–orbit coupling. For comparison and testing of warping, calculations without spin–orbit were also performed.PAOFLOW (University of North Texas, USA) [53] was used to project the computationally expensive Hamiltonian from the plane wave basis into the more efficient PAO-TB basis. The PAO-TB Hamiltonian allows the exploitation of the Fourier interpolation to obtain a smooth version of the band structure. The full BZ was divided into a 12×12×12 grid, and each voxel was searched for an isolated critical point. It was assumed that a voxel would contain at most one critical point. We chose a 12×12×12 grid to strike a balance between accuracy and computational efficiency.The locations of the critical points in k-space were then compared using operations of the symmetry group to identify unique critical points.This was done for ease of computation and for uniformity in analysis, given different crystal structures.

Fig. 1. The AFLOWπ workflow script used in this study. More details on AFLOWπ are available in Ref. [47].

In order to identify critical points,we identified k-points,where the band velocities vnare very small.

It should be noted that, in the definition of m*DOS, only one among the possible equivalent k-points is considered.The validity of the Fourier approach (which is very efficient computationally)holds for all non-warped bands but fails in the case of warped bands where the Hessian is not symmetric. In order to treat all the critical points in the same way,we calculated the inverse effective mass surface(IEMS)and used it to determine the three diagonal components of the effective mass tensor, the DOS effective mass accounting for band warping effects [6], and the band warping parameter (w) [5].

is the matrix of second partial derivatives of the energy dispersion and A is the Euler rotation matrix (see Eq. 4.46 in Ref. [57]). Here,

Fig.2. A 2D slice k = (i , j, 0.5)of the 3 components of the gradient of the energy for the bottom conduction band in the first BZ for SnTe.The purple and green represent k with large positive and negative values, respectively, for each gradient component.

Fig.3. Energy dispersion in several radial directions around the degenerate,warped Γ point k = (0 , 0, 0) near the Fermi level (EFermi) in CaCo4Sb12. The lattice parameter is indicated with a.

The sign of the integral depends on the sign of the IEMS. For saddle points, the value of the integrand in Eq. (14) approaches infinity when f2(φ, θ) approaches zero. In principle, one can integrate over the range of φ and θ where f2(φ, θ)is positive and separately where it is negative.In practice,however,it proves difficult to achieve convergence for the integral,and we do not include the calculation of the DOS effective mass for warped saddle points in our results. To the best of our knowledge,the problem of properly accounting for saddle points is still unsolved.

3. Results

There is a wide range of effects on the DOS and other material properties from the local band structure at critical points [4]. The HT procedure produces many such critical points for a given crystal, and the distribution and type of critical points ultimately account for the similarities and differences in material properties.Table 1 [58,59] summarizes some of our computations for several common materials in comparison with reference values.The effective mass was computed by all three methods discussed above:m*DOSby fitting Fourier derivatives (Eq. (6), m*F), m*DOSby fitting the IEMS(m*I),and m*DOSfrom the warped definition(m*wevaluated from Eqs. (13) and (14)). It should be noticed that, in Table 1, all calculations for masses agree for non-warped critical points,including the warped calculation of Eqs. (13) and (14). This is to be expected and is reflected in the data in general.Large discrepancies in the experimental effective masses of different materials make it difficult to directly compare with the calculated effective mass; however, our study compares favorably in most cases with the results from previous work (Table 1). All our results are included in Appendix A, specifically the band structure plots, the Hubbard U corrections, and the effective mass computed for critical points in the proximity of the Fermi level. It should be noted that experimental data are not available for critical points away from the Fermi level.

Table 1 Comparison of the IEMS fit(m*I),Fourier derivative(m*F),and warped(m*w)DOS effective masses calculated in our study with measured values(m*exp).The subscripts indicate the band-type (e: electron, h: hole, so: split-off), L and Γ are the conventional labels for two different critical points.

Fig.4. Critical points for FCC GeS in the first BZ.Each polyhedral shape relates to a different symmetry (star of the k-point) of the BZ, depending on the location of a representative point.

Consider Fig. 6 where the relative error of the direct fitting relative to the warped DOS effective mass is plotted versus the warping parameter. It can be seen from the figure that there is a general positive correlation,indicating that the warping parameter is a measure of how badly traditional methods(e.g.,direct fitting or the Fourier method) perform when evaluating effective masses;thus,the warping parameter is an indication of how badly a Taylor approximation will perform in approximating the surface of an IEMS at a warped point. This is similar to what is reported in Ref.[6]for a surface warped in Kittel’s form.Thus,the only correct value for the DOS effective mass is m*w. Fig. 7 shows the DOS of CaCo4Sb12. Critical points (van Hove singularities) are marked as red disks. The size of the disk indicates the amount of warping,as calculated in Ref. [5]. It can be seen that the DOS effective masses give rise to the shape of the non-smooth points.The particular shape can be expressed in terms of the DOS effective mass[4,6].With spin–orbit coupling included, the largest warping value w is about 1.94 in CaCo4Sb12at Γ. In general, without spin–orbit coupling, there are many more warped points (since spin–orbit coupling lifts degeneracy), and the warping is larger in that case.This study provides a computational answer to questions related to the origin and occurrence of warping effects [60]. Without exception,warping did not occur at a specific critical point without the presence of degeneracies.

Fig.5. An example of the Si band structure computed along the path suggested by Ref. [1] without spin–orbit coupling and with warping noted in red circles. The magnitude of the warping parameter w, is given as the size of the disk.

Table 2 Warping extent w and DOS effective mass at selected warped critical points calculated with 3 different methods: integral over the IEMS, ellipsoidal fit of the IEMS, and secondorder Fourier derivative.The table includes the following critical points:hole band m*e (H),heavy-hole band(m*hh),and light-hole band(m*lh)at Γ= (0 , 0, 0),the maximum of the valence band (m*max) at X = (1, 0, 0) or Γ, and selected saddle point (m*sad) at Γ.

Fig. 6. The error between the warped effective mass DOS calculation and the DOS effective mass calculated with the fit to the IEMS versus the dimensionless warping parameter (w), for the materials in this study.

Fig.7. DOS for CaCo4Sb12 including spin–orbit coupling.The van Hove singularities are identified by red disks, where the diameter of the disk is determined by the warping parameter of the critical point located at the center of the disk.The shapes of the IEMS at the warped critical points are included in the figure. Blue (red)portions of the IEMS represent negative (positive) band curvature.

In all of the cubic materials, all saddle points at Γ and at L (or other critical points) are warped. The reason for this can be seen in the requirements placed on the IEMS: The IEMS will have regions of positive and negative values due to the saddle behavior of the energy dispersion, and the cubic symmetry forces the principal directions to be the same. These two competing geometric requirements lead to warping at these points.

The HT list of critical points and effective masses has a variety of applications.For example,a full approximate analytic model of the BZ may be accomplished by approximating each critical point with the corresponding ellipsoid or saddle in the case of materials where warping is not an issue for critical points near the Fermi energy.This could be used for calculating scattering rates in Monte Carlo simulations[7]or for calculations of relaxation times[61]for different scattering processes.Moreover,a full analytic band structure would allow for possible band engineering,for the exploration of material properties,and even for the low-dimensionality reduction of thermoelectric properties [62]. A full description of the effective masses would also provide a model for the exploration of anisotropic transport considerations. As opposed to critical points that are not near high-symmetry points, lines, or planes,none of the effective masses from a quadratic expansion of the energy dispersion need to be the same. This gives rise to anisotropic transport that is all but ignored to avoid excessive complication. This procedure provides examples where anisotropy may be relevant for transport or other properties and provides a first step toward realistic analytic descriptions of band structures with these properties.

The Appendix A includes: ①the Hubbard U correction determined with the ACBN0 self-consistent protocol (Tables S1, S4,S14, and S18); ②characterization of the critical points found in the proximity of the Fermi level for AlSb, AlP, GeS, GeSe, GeTe,PbS, PbSe, PbTe, SnS, SnSe, SnTe, CoSb3, CaCo4Sb12, BaCo4Sb12, Bi2-Se2Te, Bi2Te3, Bi2Te2Se, Bi2Se3, Bi2Te2S, and Bi2Se2S (Tables S2, S3,S5, S7–S13, S15–S17, and S19–S24 in Appendix A); and ③band structure plots (Figs. S1–S20 in Appendix A).

4. Conclusions

We have designed and applied a two-level HT calculation with the aim of finding and characterizing critical points in the BZs of several materials.This search identified features of the band structures, En(k ), that are not visible in 2D plots but nonetheless contribute to electronic transport, DOS, and other material properties. We used three different methods to compute the DOS effective masses at these critical points:namely,the direct method,the Fourier derivative method,and the warped method.At analytic(ellipsoidal) critical points, all approaches agreed with each other and aligned well with literature values,where available.This study also identified many warped critical points and used Eqs. (13) and(14)to correctly calculate the DOS effective mass.At these warped points, we provided reliable values for m*DOSusing the theory reported in Ref.[5].All three methods disagreed,but the‘‘warped”method was correct. Standard 2D band structure plots are often inadequate to convey the complicated nature of the effective masses at important critical points. Our HT procedure overcomes this issue and facilitates comparison with experiments and the development of analytical models, as well as band structure engineering and the tuning of properties for technological applications.

Acknowledgments

The authors would like to thank Lorenzo Resca for helpful discussions regarding the symmetry properties of the Brillouin zones,along with other aspects concerning warping. This work was performed in part through computational resources and services provided by the Institute for Cyber-Enabled Research at Michigan State University.Nicholas A.Mecholsky would like to acknowledge financial support from the Vitreous State Laboratory.

Compliance with ethics guidelines

Andrew Supka, Nicholas A. Mecholsky, Marco Buongiorno Nardelli, Stefano Curtarolo, and Marco Fornari declare that they have no conflict of interest or financial conflicts to disclose.

Appendix A. Supplementary data

Supplementary data to this article can be found online at https://doi.org/10.1016/j.eng.2021.03.031.