Relaxation time approximations in PAOFLOW 2.0.

Regardless of its success, the constant relaxation time approximation has limited validity. Temperature and energy dependent effects are important to match experimental trends even in simple situations. We present the implementation of relaxation time approximation models in the calculation of Boltzmann transport in PAOFLOW 2.0 and apply those to model band-structures. In addition, using a self-consistent fitting of the model parameters to experimental conductivity data, we provide a flexible tool to extract scattering rates with high accuracy. We illustrate the approximations using simple models and then apply the method to GaAs, Si, [Formula: see text], and [Formula: see text].

Introduction

Over the years, first-principles calculations have become a complementary tool for the experimental research aiming to discover and understand materials with relevant electronic properties. This has greatly improved also the understanding of the electronic transport properties which are crucial in applications ranging from electronics to energy conversion. With electronic transport coefficients, we indicate the response function to an applied external electric field and/or temperature gradient. The quantities of interest are the electric conductivity, , the electronic therma conductivity, , and the Seebeck coefficient, S, which couples the heat and charge fluxes. These quantities are, in the linear regime, constant tensors and include dissipation phenomena.

The most standard approach in calculating the transport coefficient is to use the semi-classical Boltzmann theory within the constant relaxation time approximation (CRTA). In most cases, this has been used without much analysis. Though it works for certain systems, recent research has shown that the CRTA have resulted in wrong predictions, missing vital information to understand the transport properties[1-3]. Moreover, the introduction of an arbitrary constant relaxation time severely limits the predictive capabilities of first-principles band-structures. In this work, we investigate how the CRTA affects the electric conductivity of well known materials by considering different scattering models that include energy and temperature dependence. Though methods exist for ab initio calculation of electron–phonon relaxation times[4,5], the calculation of the electron–phonon matrix requires extremely dense and point meshes. Therefore, prohibitively high computational costs make these techniques of limited practicality, especially when aiming to data driven high-throughput approaches. We chose to combine accurately interpolated band-structures and simplified mathematical models of the scattering phenomena in order to explore the consequences on the transport coefficients beyond CRTA and parabolic bands. Extrinsic scattering mechanisms (impurities, grain boundaries, alloy disorder) contribute significantly to the transport properties in a system. Indeed, these extrinsic scattering mechanisms may be often tuned during the synthesis of the system[6] and often beyond the power of predictive computational methods. The contribution of intrinsic scattering mechanism in a system, on the other hand, can be more easily addressed to obtain optimized performances. In addition to the insight offered by relaxation time and transport data studies offered by other methods[7,27], in this work, we introduce a self-consistent fitting of transport properties to experimental data which will give us an understanding of the temperature dependence of various scattering mechanisms for specific experimental samples. These relaxation time approximation (RTA) models are coded in our recently released PAOFLOW package[8,9]. In this paper, we will discuss in detail the theory and implementation of the RTA models in the newest release of the software, and illustrate the automated workflow while, at the same time, documenting the influence of scattering phenomena in the band-structure of cubium, graphene, and selected materials: Si, GaAs, MgSb, and CoSb.

Methods

PAOFLOW

PAOFLOW is a software tool to efficiently post-process standard first-principles electronic structure plane-wave pseudopotential calculations in order to promptly compute, from interpolated band-structures and density of states, several quantities that provide insight on transport, optical, magnetic and topological properties such as anomalous and spin Hall conductivity, magnetic circular dichroism, spin circular dichroism, and topological invariants. The methodology is based on the projection on pseudo-atomic orbitals (PAO) discussed in detail in Refs.[10-12].

Accurate PAO Hamiltonian matrices can be built from the direct projection of the Kohn–Sham (KS) Bloch states onto a chosen basis set of fixed localized functions. The Hamiltonian for a specific material, , is computed in real space using atomic orbitals or pseudo atomic orbitals from the pseudopotential of any given element. The key, in this procedure, is in the mapping of the ab initio electronic structure (solved on a well converged and large plane waves basis set) into tight-binding (TB) formalism that precisely reproduces a selected number of bands of interest. The crucial quantities that measure the accuracy of the basis set are the projectabilities ( is the operator that projects onto the space of the PAO basis set, as defined in Ref.[11]) which indicate the representability of a Bloch state on the chosen PAO set. Maximum projectability, , indicates that the particular Bloch state can be perfectly represented in the chosen PAO set; contrarily, indicates that the PAO set is insufficient and should be augmented. Once the Bloch states with good projectabilities have been identified, the PAO Hamiltonian is constructed asHere E is the diagonal matrix of KS eigenenergies and A is the matrix of coefficients obtained by projecting the Bloch wavefunctions onto the PAO set. Since the filtering procedure introduces a null space, the parameter is used to shift all the unphysical solutions outside a given energy range of interest. The procedure in Eq. (1) is recommended for most cases.

Band-structure interpolation on arbitrary Monkhorst and Pack (MP) k-meshes for the integration in the Brillouin zone (BZ) are at the very core of the ability of PAOFLOW to provide high-precision electronic structure data. Indeed, the TB Hamiltonian can be Fourier transformed from real space representation to the k-space and interpolated using an efficient procedure based on a zero-padding algorithm and fast Fourier transform routines.

The same accuracy defined by the projectabilities is conserved in this process. The expectation values of the momentum operator, which is the main quantity in the definition of the transport coefficients, is given bywith being the real space PAO matrix and the Bloch’s functions[13].

Boltzmann transport

In PAOFLOW the electrical conductivity is evaluated by solving the semi-classical Boltzmann equation (BTE) that describes the evolution of the distribution function f of an electron gas under external electric field and in presence of scattering mechanisms[14-16]. In the scattering-time approximation, the conductivity tensor can be expressed as an integral over the first BZ:where is the relaxation-time, is the i-th component of the electron velocity corresponding to the n-th band for each k-point in the BZ ( is derived by the diagonal of the momentum matrix element, Eq. (2)), is the equilibrium distribution function, and E is the electron energy.

Generalizing Eq. (4) it is also possible to define analogue expressions for the Seebeck coefficient S and the electron contribution to thermal conductivity . Following the notation of Ref.[17], we introduce the generating tensors ():where indicates the dyadic product, the band-structure, and is the chemical potential. The coefficients , S and can be expressed as follows:where T is the temperature. Our formalism based on PAO-TB performs the computations of the band-velocities and avoids issues with possible band-crossing. In addition, from Eqs. (4–5), it is evident that the evaluation of the transport properties requires an accurate integration over a fine grid of k-point in the BZ which becomes a trivial task using the TB representation from the PAO projections and Eq. (3).

Relaxation time approximation

The most common implementations of the Boltzmann transport equations assume the scattering time to be a constant (CRTA). A constant factors out of Eq. (5) and, thus, the method returns the quantities and (in the Seebeck coefficient cancels out). Clearly this is a severe approximation for a quantity that is expected to be significantly temperature-dependent and energy-dependent (typically via a power law). A direct estimate of the dependence of on energy and temperature is an important complement to any transport study, and would provide, even if at a phenomenological level, important insight into the relevant scattering mechanism present in any given system. Moreover, a direct comparison with existing experimental data would provide an extra layer of characterization for real world applications.

The approach implemented in PAOFLOW is based on the work of Jacoboni et al.[18] and recently included in the BoltzTrap[19] framework by the group of Fiorentini[20]. They employed analytical energy-dependent expressions for the relaxation time, which were developed on the basis of known semiclassical theories and include the most important mechanisms of electron scattering by acoustic phonons, polar-optical phonons, and charged impurities.

Acoustic phonon scattering is treated within the elastic deformation potential approach in the long-wavelength acoustic-phonon limit,where E is the band structure energy considered from the top of the valence bands and T is the temperature. All other parameters are defined in Table 1.

Table 1

Symbols and units for the scattering parameters required in various scattering models.

Parameter	Symbol	Units
Mass density	ρ	kg/m^3
Lattice constant	a	m
Low freq. dielectric constant	ϵ_0	–
High freq. dielectric constant	ϵ_∞	–
Acoustic velocity	v	m/s
Effective mass ratio	m^∗	–
Acoustic deformation potential	D_ac	eV
Optical deformation potential	D_op	eV
Optical phonon energy	ħω_op	eV
Number of impurities	n_I	cm^-3
Charge on impurity	Z_I	–
Piezoelectric constant	p	C/m^2

Similarly to the assumptions that were used in acoustic phonon scattering, we model optical phonon scattering with an elastic deformation potential ():The first term in the denominator of Eq. (89) represents the absorption of optical phonons and the second term represents the emission of the optical phonons. The probability of emission of a phonon when is zero since the electron does not have enough energy to emit the phonon and this is represented by the Heaviside step function included in the second term. represents the number of optical phonons. Polar optical scattering is modeled following Ridley:[21]where the sum is over all longitudinal-optical phonons, with energy ; the functions A, B, C, and Z are omitted for brevity and can be found in Appendix I of the Supplementary Information.

Note that the deformation potential theory used to model the electron phonon interaction assumes that the atomic displacements can be described by long wavelength atomic and optical waves[18,37]. This is a wide spread assumption commonly used in semiconductors since carriers are typically confined to narrow energy windows around the band extrema (therefore the scattering mechanism will be dominated by long wavelength phonons). However, this approximation should be used with caution in systems where short wavelength phonon scattering is known to be important. For impurity scattering we use the Brooks–Herring approach:[22]Finally, in semiconducting compound the strain induced by acoustic phonons creates a piezoelectric field. This piezoelectric scattering is modelled as in Ref.[18].where and the piezoelectric effect is captured by the piezoelectric constant, p.

The global relaxation time is then obtained using Matthiessen’s rule:

Simple models and the parabolic band approximation

In order to quantify the improvement of a richer RTA, and to gain a better understanding how varying various parameters affect the overall transport properties of a system, we start with two simple TB models: cubium and graphene. Cubium was chosen as representative of a 3D solid with quasi-parabolic bands (near the BZ center, and graphene for its 2D character and its linear dispersion at the Fermi level. The TB Hamiltonian for a system with two atoms per unit cell with contributions from a single orbital is given bywhere t is the first nearest-neighbor hopping parameter and defines the band gap of the band-structure. and runs over nearest neighbors. For the cubium, because of its simple structure, the vectors for the six nearest neighbours are , where a is the lattice constant, so that:Similarly for graphene, the vectors for the three nearest neighbors are , andSetting to 0 eV reproduces a graphene like band-structure where the bands show a linear dispersion at the Dirac point K in the BZ.

Since the energy and temperature dependence of the functional form of the scattering times were obtained using a parabolic band approximation, it is useful to examine the validity of such an approximation within electronic transport. The formulas for transport coefficients for a parabolic band within the CRTA, derived in detail in Ref.[17], have been compiled in Eqs. (14)–(19).

In the insulating regime for parabolic bands, when and in the metallic regime for parabolic bands, when where, represents either a band-edge minimum or maximum, and is +1 for a conduction-like band, and -1 for a valence-like band, and are the x, y and z components of the effective mass, is the chemical potential of interest, and n is the charge carrier concentration. We use the cubium model (Eq. 21) and compare the transport properties to those of the parabolic bands in Fig. 1:

Figure 1

Band-structure of a cubium with two bands (top left panel, dashed line) and a parabolic fit near point (solid line). In the top right, bottom left, and bottom right panels, the Seebeck coefficient, the conductivity, and the carrier concentration, respectively, are reported.

As apparent from the top left panel of Fig. 1, the cubium bands start deviating from the parabolic bands at  1.5 eV. Subsequently, all the transport properties for the parabolic model and the cubium band-structure are expected to match from 0 eV to 1.5 eV as the cubium band is a good approximation of a parabolic band in this range. This is reflected in the rest of the panels of Fig. 1. Since the parabolic band approximation is no longer valid in the cubium model after this limit, the transport properties of cubium start to deviate from those calculated using Eqs. (14)–(19). This deviation from parabolicity in transport properties is slightly enhanced at 4 eV due to contributions from the flat feature of the bands of cubium around the X point of the band-structure.

In Fig. 2, we consider a graphene like band-structure with a band gap of 0.5 eV. As expected, graphene bands do not follow a parabolic approximation, except very close to the Fermi surface. Therefore, the conductivity calculated using a parabolic band approximation is able to reproduce the conductivity calculated using BTE only extremely close to the valence and conduction band edges.

Figure 2

Band-structure of a graphene model and parabolic fit at K (left panel, dashed and solid line, respectively). The right panel shows the corresponding conductivities computed with Eqs. (14), (17), and PAOFLOW.

Simple models beyond CRTA

Simple models allow to investigate the effect of improving the CRTA with minimal computational overload and monitor the consequences in physically transparent scenarios. Our simplified approach is conducive to an exploration of the parameters’ space of the RTA models which can be then used with more realistic band-structures. We can control the parameters of the models to enhance one or the other scattering mechanisms by varying the values of the deformation potentials, or of the relevant optical frequencies or of the sound velocity, to name a few. This can give insight into the potential design of materials with optimal properties for any given application.

A first observation involves the comparison between the conductivity calculated within CRTA and other other RTA models. We present results for graphene and cubium, Fig. 3. Experimentally, the distinction between semiconducting and metallic behavior (including the case of heavy doping) is understood in term of the temperature dependence of the conductivity. Samples whose conductivity increases with temperature are semiconducting and samples whose conductivity decreases with temperature are metallic. In semiconductors, the increase in the number of charge carriers prevails over the reduction of the relaxation time; in metals, the reduction of induces the reduction of .

Figure 3

Electronic conductivity as a function of temperature at various chemical potentials (representing the doping level, referred to the middle of the gap). The left panel corresponds to graphene (with zero band gap) and the right panel corresponds to cubium with a band gap of 0.5 eV. The red lines denote the conductivity calculated using the CRTA whereas the black lines correspond to the conductivity calculated using the RTA; different markers correspond to different chemical potentials as in the legends.

In the case of graphene (a quasi-linear dispersion), the CRTA ( s) provides temperature-independent conductivity in the metallic case and monotonically increasing conductivity in the semiconducting regime (the position of the chemical potential wrt the band edge determine the regime from the electronic structure point of view): only the variation of the carrier density due to temperature is captured in the calculation (Fig. 3, left panel, red lines). When applied to a two-band cubium with a forbidden energy gap of 0.5 eV, the same phenomenology is obtained in the CRTA (Fig. 3, right panel, red lines). Let’s consider a specific RTA model constructed using and optical phonon and acoustic phonon scattering mechanisms (see section “Relaxation time approximation”). In graphene, the chosen RTA model introduce dissipation phenomena that shorten as the temperature increases: this induce a reduction in conductivity. In the cubium model, we recover the experimental evidence of decreasing conductivity as function of temperature in the heavily-doped (metallic case) and increasing conductivity as function of temperature in the semiconducting regime.

Relaxation time models and experimental conductivity

The scattering models described in section “Relaxation time approximation” have limited validity when extrinsic effects such as size of dopants or impurities, edge effects, etc. can significantly affect the scattering time. The electronic conductivity calculated using the relaxation times from the models, however, provide a framework for comparison to experiments and provide insight on the wide variations in experimental conditions and doping. We propose a modified Mathiessen’s rule as:In the above equation, , , , and are correcting functions that are fitted to reproduce the experimental conductivity. The fitting procedure uses the sequential least squares programming (SLSQP)[23] method, which allows for constrained non linear optimization of the fitting functions. A more detailed discussion of the fitting procedure is presented in Appendix 2 of the Supplementary Information. To demonstrate the effectiveness of this approach and its implementation in the PAOFLOW package, we present results for four prototypical systems: GaAs, Si, MgSb and CoSb.

Computational details and implementation in PAOFLOW

The calculation of the Boltzmann transport with the modified RTA models is implemented in PAOFLOW and follow a standard algorithmic flow. PAOFLOW requires a few basic calculations performed with the Quantum ESPRESSO (QE) package[24,25]. The first (self-consistent) run generates a converged electronic density and Kohn–Sham (KS) potential on an appropriate Monkhorst and Pack (MP) k-point mesh. The second (non self-consistent) one evaluates eigenvalues and eigenfunctions on a larger MP mesh and often for an increased number of bands. After these preliminary steps PAOFLOW’s most fundamental procedure is the construction of accurate PAO Hamiltonians following the theory outlined in section “PAOFLOW”.

The density functional theory (DFT) calculations for Si were performed using a local density approximation (LDA). A kinetic energy cut-off of 18 Ry (180 Ry cut-off for the charge density) and a Monkhorst–Pack k-point mesh were used for the non self-consistent calculation. This was further increased to a grid using PAOFLOW’s Fourier interpolation method in order to accurately integrate transport tensors.

The DFT calculations for GaAs were performed using the generalized gradient approximation (GGA) in the parametrization of Perdew, Burke, and Ernzerhof (PBE). Projector augmented wavefunctions (PAW) pseudopotentials were used to treat the ion-electron interactions. A kinetic energy cut-off of 60 Ry (600 Ry cut-off for the charge density) and a Monkhorst–Pack k-point mesh were used. This was further increased to a grid using PAOFLOW (The convergence of electronic conductivity as a function of k-grid size is shown in Appendix 4 of the Supplementary Information).

The DFT calculations for were performed using a GGA functional in the parametrization of PBE. PAW were used to treat the ion-electron interactions. A kinetic energy cut-off of 45 Ry (450 Ry cut-off for the charge density) and a Monkhorst–Pack k-point mesh were used. This was further increased to a grid using PAOFLOW.

The DFT calculations for were performed using a LDA functionals. A kinetic energy cut-off of 45 Ry (450 Ry cut-off for the charge density) and a Monkhorst–Pack k-point mesh were used. This was further increased to a grid using PAOFLOW. The pseudopotentials for all the atomic species were obtained from pslibrary1.0.0[26].

The experimental parameters used in the calculations of the relaxation times for different systems have been listed in Table 2. The energy, E is taken from the original DFT Hamiltonian processed by PAOFLOW. The calculations are done for user defined ranges of temperature T. An example workflow is discussed in Appendix 3 of the Supplementary Information. The relaxation time as a function of energy obtained using models have been compared to those obtained from first principles[27-29] in Appendix 4 of the Supplementary Information. A reasonable qualitative agreement was observed with the additional benefit of low computational cost.

Table 2

Symbols and units of the parameters to be input in calculation of scattering models. The values for Si and GaAs are obtained from Ref.[18], for from Ref.[20] and for from Refs.[30,31].

Symbol	Units	Mg_3Sb_2	GaAs	Si	CoSb_3
ρ	kg/m^3	3.9×10^3	5.3×10^3	2.3×10^3	7.8×10^3
a	m	8.7×10^-10	5.6×10^-10	5.4×10^-10	9.1×10^-10
ϵ_0	–	26.7	13.5	11.7	33.5
ϵ_∞	–	14.2	11.6	–	25.6
v	m/s	2.7×10^3	5.2×10^3	6.6×10^3	3.3×10^3
m^∗	–	0.3	0.7	0.29	3
D_ac	eV	6.5	7	9.5	5
D_op	eV	–	–	8×10^10	1×10^11
ħω_LO	eV	[0.0205,0.0248,0.031]	[0.03536]	–	[0.0264]
p	C/m^2	–	0.16	–	–

GaAs

The scattering models were implemented on n-type GaAs for two different carrier concentrations, and and was used to calculate conductivities.

The scattering rates as a function of temperature are depicted by the solid lines in Fig. 4. The dominant scattering mechanisms were determined from Refs.[32] and[33]. The calculated conductivities were then compared to the experimental values as shown in the inset of Fig. 4. The fitting procedure is performed as well and the scattering rates obtained as a result of the fitting procedure are shown by the dashed lines in Fig. 4. The fitting procedure produces negligible changes to the scattering rates, signifying that the original models themselves represent the scattering rates in GaAs well.

Figure 4

A comparison of the scattering rates in GaAs obtained using the original scattering models (labelled RTA) and those obtained using the fitting procedure (labelled RTA fitted) for samples of two different doping concentrations. The left panel corresponds to an n-type sample with a doping concentration of while the right panel corresponds to an n-type doping of . The inset shows the electrical conductivity calculated using both the original and the fitted scattering rates.The experimental data have been obtained from Ref.[32].

Reference[32], from which the experimental data have been obtained, uses analysis of their Seebeck and Hall coefficient data and notes that the relative weight of polar scattering increases with increasing temperature, whereas the contribution of impurity scattering decreases with increasing temperature. This is confirmed by our results as well.

Si

Electronic conductivities were calculated for intrinsic Si, and n-type Si with doping concentrations of cm and cm. Similar to GaAs, the solid lines in Fig. 5 represent the scattering rates calculated from the original models while the dashed lines represent the scattering rates obtained using the modified Mathiessen’s rule. It is clear from the results that the original scattering rates require significant modification for the calculated conductivities to match experiments.

Figure 5

Conductivities and scattering rates for the different Si samples in a low temperature regime. The left panel shows the data for a sample with n-type doping with a carrier concentration of  cm while the right panel shows the data for a sample with n-type doping with a carrier concentration of  cm. The inset shows the goodness of fit of theoretical electrical conductivity to experiments resulting from the fitting procedure as well as the electrical conductivity calculated using the original scattering models

Low temperature experimental data for transport properties in Si were obtained from Ref.[34]. There, the authors discuss the electrical conductivity of the heavily doped sample ( cm) being weakly temperature dependent due to the large number of impurity atoms forming an impurity band. They state that in contrast, the weakly doped  cm exhibits an exponential behavior of electrical conductivity due to the freeze out of impurities at low temperatures. This behavior is fully captured by our models as well. In panel 1 of Fig. 5, for  cm the fitting procedure produces negligible change to the scattering rate due to impurity scattering in order to match experimental data. However, in panel 2 of Fig. 5, for  cm we see the contribution of impurities to the overall scattering rate is underestimated by the original models, but is rectified by our fitting procedure which significantly increases the scattering rate due to impurities. The dominant scattering mechanisms are the optical and acoustic phonons scattering whose values require significant correction by the fitting procedure.

MgSb

, a well studied thermoelectric material was chosen as a test system in order to verify our implementation of the various scattering models as well as the fitting procedure. The effect of scattering rates on transport properties in and comparison to experiments have been extensively carried out in Ref.[20]. Their results show that transport properties calculated with scattering models are in good agreement with experiments. This is confirmed by our implementation of scattering models and the subsequent fitting procedure. As shown in Fig. 6, the fitting procedure produces negligible modifications to the original scattering models in order to match calculated electrical conductivity to experimental data.

Figure 6

A comparison of the scattering rates in MgSb obtained using the original scattering models (RTA) and those obtained using the fitting procedure (RTA). The left panel corresponds to an n-type sample with a doping concentration of while the right panel corresponds to an n-type doping of . The inset shows the electrical conductivity calculated using the respective scattering rates.The experimental data have been obtained from Ref.[20]

CoSb

The conductivity for the p-type samples of seem better represented by the scattering models than that of those for n-type samples. As seen in Fig. 7, the scattering rates for n-type sample require significant corrections at higher temperatures in order to reproduce experimental results.

Figure 7

A comparison of the scattering rates in CoSb obtained using the original scattering models (RTA) and those obtained using the fitting procedure (RTA) are shown. The top left panel corresponds to a p-type sample with a doping concentration of while the top right panel corresponds to a p-type doping of . The bottom left panel corresponds to an n-type sample with a doping concentration of while the bottom right panel corresponds to an n-type doping of . The inset shows the electrical resistivity calculated using the both the original and fitted scattering rates. The experimental data have been obtained from Ref.[35].

CoSb, a well known thermoelectric, is studied over a wide range of temperatures. Caillat et al.[35] analyzed their experimental mobility data for p-type samples and suggests that the dominant scattering mechanism, at least below 500 K is acoustic phonon scattering since the mobility followed a T behaviour. However, our models seem to suggest that the dominant scattering mechanism is optical phonon scattering for p-type samples and depending on the doping levels, maybe acoustic phonon scattering or optical phonon scattering for n-type samples. This inconsistency is also noted by Kajikawa[36] who carried out analysis of p-type CoSb within a two-valence and two-conduction band model.

Conclusion

We have implemented relaxation time models that allow calculation of conductivities beyond the constant relaxation time approximations and are able to provide reasonable agreement to experimental conductivities in various systems. Moreover, unlike the CRTA, it allows for a quantitative and qualitative description of the scattering mechanisms themselves. We introduce an automated self consistent fitting procedure that allows one to see how various the base scattering models need to be tuned in order to reproduce experimental conductivities. This is highly advantageous in determining sample specific scattering properties which is beyond the scope of the base models. Future work includes the implementation of mobility using the BTE and therefore the extension of PAOFLOW as a tool to analyze experimental mobility data to extract scattering time trends.

Supplementary Information.