Electronic transport computation in thermoelectric materials: From ab initio scattering rates to nanostructures

Over the last two decades a plethora of new thermoelectric materials, their alloys, and their nanostructures were synthesized. The ZT figure of merit, which quantifies the thermoelectric efficiency of these materials increased from values of unity to values consistently beyond two across material families. At the same time, the ability to identify and optimize such materials, has stressed the need for advanced numerical tools for computing electronic transport in materials with arbitrary bandstructure complexity, multiple scattering mechanisms, and a large degree of nanostructuring. Many computational methods have been developed, the majority of which utilize the Boltzmann transport equation (BTE) formalism, spanning from fully ab initio to empirical treatment, with varying degree of computational expense and accuracy. In this paper we describe a suitable computational process that we have recently developed specifically for thermoelectric materials. The method consists of three independent software packages that we have developed and: 1) begins from ab initio calculation of the electron-phonon scattering rates, 2) to then be used within a Boltzmann transport simulator, and 3) calculated quantities from BTE are then passed on to a Monte Carlo simulator to examine electronic transport in highly nanostructured material configurations. The method we describe is computationally significantly advantageous compared to current fully ab initio and existing Monte Carlo methods, but with a similar degree of accuracy, thus making it truly enabling in understanding and assessing thermoelectric transport in complex band, nanostructured materials.


I. Introduction
Thermoelectric generators (TEGs) are solid-state devices able to convert the heat flow arising from temperature gradients directly into electricity. They have the potential to offer a sustainable path for power harvesting from a variety of industrial sectors at power levels from microwatts to tens/hundreds kW, and even MW. Their impact could be widespread across many applications including medical, wearable electronics, building monitoring, the internet of things, refrigeration, thermal management, space missions, transportation, and various industrial sectors [1]. However, high prices, toxicity, scarcity, and low efficiencies of the prominent thermoelectric (TE) materials are currently hampering their large-scale exploitation. On the other hand, progress on these materials has been rapidly expanding over the last two decades. Novel concepts and improved understanding of materials synthesis have provided new opportunities for the enhancement of the thermoelectric conversion efficiency across many materials [1]. The thermoelectric figure of merit ZT, which quantifies the ability of a material to convert heat into electricity, has more than doubled compared to traditional values of ZT~1, reaching values of ZT > 2 in several instances across materials and temperature ranges, and even approaching 3 in some cases [1,2,3,4,5].
The TE performance is quantified by the ZT figure of merit as ZT = σS 2

T/(κe+κl),
where σ is the electrical conductivity, S is the Seebeck coefficient, T is the absolute temperature, and κe and κl are the electronic and lattice parts of the thermal conductivity, respectively. The product σS 2 in the numerator of ZT is called the power factor (PF). The recent improvements in ZT are mostly attributed to drastic reductions of the lattice thermal conductivity in nanostructured materials and nanocomposites, which has reached amorphous limit values at κl = 1-2 W/mK and below [1,6,7,8,9]. One of the most successful strategies to reduce thermal conductivity is to hierarchically nanostructure the materials, by introducing different types and sizes of nano-features, which would scatter more effectively different groups of phonon mean-free-paths. Nanostructuring has brought tremendous success, but it is gradually running out of steam in providing further reductions to thermal conductivity. It is becoming increasingly clear that any further benefits to ZT must now come from power factor improvements. Current research efforts in improving the power factor aim towards identifying materials with favorable bandstructure features, such as resonant states and low-dimensional 'like' features within bulk materials [10,11], or bandstructure engineering such as bandconvergence strategies [11,12,13,14]. Indeed, most promising TE materials have complex bandstructures, with multiple valleys and bands that extend in the entire Brillouin Zone (BZ).
Due to the vast material possibilities, the exploration for PF enhancing band features needs to be led by advanced computational studies, which elucidate their role in electronic and TE transport. The most common computational approach is to extract the electronic bandstructures using ab initio methods, and afterwards use the Boltzmann Transport Equation (BTE) to extract the TE coefficients. For the scattering rates used within the BTE the most common approach is the constant relaxation time (CRT) or constant mean-free-path (CMFP) approximations, such as using BoltzTrap code [15], where a CRT around 10 fs, or a CMFP of around 5 nm, are routinely employed [15,16].
These approximations have the advantage of being computationally efficient, but have the disadvantage of providing uncertain and rather arbitrary outcomes, both quantitatively and qualitatively, e.g. with respect to materials ranking, temperature trends, and optimal carrier density [17,18]. However, we know that in these materials, the electronic scattering processes have complex energy, momentum and band dependencies [18].
More elaborate computational methods for the treatment of electron-phonon (e-ph) scattering vary, depending on the required accuracy and computational complexities, i.e. by only considering the energy dependence of acoustic phonons analytically [19]; by considered elastic scattering by acoustic phonons in a full band approach using deformation potentials [20]; by extracting deformation potentials from band shifts after applying stress along specific directions and using those to form scattering rates [21]; by using a full band approach and a numerical scheme for electron-phonon scattering based on calculating the ab initio electron-phonon coupling matrix elements, but with constant optical phonon frequencies [22,23]; and by full ab initio EPW+Wannier e-ph scattering rate calculations using the EPW code [24]. Other than for 2D materials, it is computationally challenging to utilize the entire phonon and electron dispersions for full transport calculations, thus, it is common to employ reasonable approximations [23]. Overall, however, the majority of current TE simulation works (justifiably) tend to heavily sacrifice accuracy for computational speed.
The need for theory to lead experiment, the search for new materials, and detailed understanding of the transport physics that would lead to performance optimization, requires electronic transport tools that provide confidence in accuracy, as well as computational robustness. We have recently developed such tools, which fill the gap of higher accuracy, in the inevitable expense of moderate computational cost, but still much lower compared to fully ab initio simulations. Our method is based on ab initio calculations of a selected few electron-phonon matrix elements, out of which we form deformation potentials and scattering rates based on the deformation potential theory. We then use those within our newly developed full-band BTE code, ElecTra, which provides the electronic properties of the bare material. This information is then passed on to a third simulator, a Monte Carlo real space ray-tracing code, specifically designed to address challenges for simulating TE transport in nanostructures. These three main sections are presented separately in the detailed works in Refs. [25,26,27], but in this paper we provide an overview of the entire process.
The paper is organized as follows. In Section II we describe the BTE code. In Section III we describe the extraction of the deformation potentials from ab initio calculations. In Section IV we present how we use the BTE and ab initio information in the MC code and present an example for transport simulations in a nanostructured material.
Finally, in Section V we conclude.

II. Boltzmann transport using the ElecTra code
To evaluate the transport properties of the base material, we use our newly developed simulator, ElecTra, that considers all appropriate scattering mechanisms (acoustic phonon, optical phonon, ionized impurity scattering), and the full energy, momentum, and band dependence of the relaxation times. Here we provide a few details on the process flow of the code, whereas more details can be found in the dedicated ElecTra paper and the user's manual [26,28].

Linearized Boltzmann Transport Equation (BTE) formalism
The TE transport coefficients are extracted within the Linearized Boltzmann Transport Equation (LBTE) formalism under the relaxation time approximation as [29,30]: where ( ) is the Transport Distribution Function (TDF) defined below in Eq. (2), EF, T, q0, kB, and f0, are the Fermi level, absolute temperature, electronic charge, Boltzmann constant, and equilibrium Fermi distribution, respectively.
The TDF is expressed as a surface integral over the constant energy surfaces, , for each band, and then summed over the bands, as [29,30,31]: where is a state on the surface and is its corresponding surface area element, computed as explained below below. ( , ) is the i-component of the band velocity of the transport state, ( , ) is its momentum relaxation time, |⃗ ( , ) | is its density-of-states (DOS), and s is the spin degeneracy.
The relaxation times for each individual scattering mechanism are combined following using Matthiessen's rule for each (k, ) state, to compute the comprehensive TE coefficients. Also, the overall energy-dependent relaxation time  and mean-free-path (mfp)  are returned, both per-band, per scattering mechanism, and overall for all mechanisms. These are computed as: These are meaningful full-band transport quantities, to be further exploited in subsequent simulations (such as Monte-Carlo).

Electronic structure quantities
The first step in the simulation approach is to obtain the electronic structure using Density Functional Theory (DFT). ElecTra's interface can take as input a '.bxsf' format file [32,33] enabling to interface with any DFT code that provides the electronic structure data in this format. This is the format used by 'XCrysDen' [33], which is compatible with a number of DFT codes. For example, QE, CASTEP and VASP among others, have the fs, c2x, and vasp2x_fs routines, respectively, to save the computed DFT bandstructure in this format.
Then ElecTra builds the constant energy surfaces to grasp the 3D details of the bandstructure. For this, we map the function E(k) into a function k(E). The 3D mesh in the k-space is scanned, the mesh elements crossed by the constant energy surfaces are identified, and the coordinates of the points on these surfaces are computed. For every element of the k(E) mesh, the E(k) is linearly interpolated from the original mesh, and the three components of the band velocity ( , , ) are computed with the contragradient method [34]. ElecTra offers two schemes to perform this step: i) A local triangulation on the k-space mesh elements which are crossed by the surface of the constant energy of interest. Here the elemental surface area, ( , , ) , which will define the DOS of that specific state, is computed using Heron's formula. ii) An easier approximate method which samples the nearest-neighbour points on the k-mesh (NN-scan). Here the code scans all the k-points on the energy surface, checks for their nearest neighbours, then uses the distances between these k-points to approximate the dAk surface element areas. In the latter, ElecTra does not resolve constant energy surface elements, but acquires a collection of points on the energy surface of interest. Although this is an approximation, as it detects only the points along the edges of the k-mesh elements, it is around 15 to 30 times faster, but without noticeable penalty in the results compared to triangulation, either for isotropic or anisotropic bands. For further details about these methods we refer the interested reader to the original Electra paper and the user's manual [26,28].
The k-state-dependent DOS is then defined as Note that we have used energy surface elements that are extracted after we construct the constant energy surfaces, which is an alternative way compared to codes such as BoltzTrap which uses volume integral and delta-functions represented by smeared Gaussians [15].
The method has been systematically validated with details reported in the manual [28] for various example cases: i) comparisons with parabolic and non-parabolic bands for which the DOS, velocity and TDF solutions for isotropic scattering mechanisms are analytical, ii) comparisons for the mobility of materials cases for which the mobility is experimentally well-known, i.e. Si, Ge, SiGe, and GaAs, and iii) comparisons with quantities (i.e. DOS) extracted from existing codes. The latter is depicted in Fig. 1 for TiCoSb and Mg3Sb2. Espresso and ElecTra, the latter from the two methods implemented (nn-sampling and triangulation). The inset shows constant energy surfaces of one of the valence bands at E = -0.12 eV below the valence band edge. (b) DOS for the conduction band of Mg3Sb2, and comparison between the results obtained from BoltzTraP [15] and ElecTra [26]. The inset depicts the constant energy surface for the Mg3Sb2 conduction band at E = 0.1 eV above the band edge.

Scattering rates and transport:
For each transport state (k,n,E) and each scattering mechanism ms, the corresponding where the sum runs over all the allowed final states k' of the same carrier spin [30,31].
|Sk,k'| is the transition rate between the initial k and final k' states, computed as detailed by Eq. 6 below for the different mechanisms. The (1 − ) term is an approximation for the momentum relaxation time, which is used to solve the BTE in the closed form, as commonly done in the literature when computing the transport coefficients [35]. ElecTra computes the scattering rates for the different scattering mechanisms as [30,31]: Optical Phonon' and describes the inelastic/anisotropic scattering of charge carriers with polar phonons, which here is treated as both intra-and inter-band [30]. ElecTra allows different phonon frequencies for all these inelastic processes separately, for example, for each non-polar and polar phonon branches, different frequencies can be used. IIS stands for 'Ionized Impurity Scattering' and describes the elastic scattering rate due to ionized dopants, for which the user can choose both intra-and/or inter-band transitions. "Alloy" represents the alloy scattering due to intrinsic disorder in alloys or solid solutions and is considered both intra-and inter-band [35]. k and k' are the wave vectors of the initial and final states. "-" and "+" in Eqs. 6b-7d indicate the phonon absorption and emission processes, respectively. These type of transition processes are illustrated in Fig. 2a. We allow for the directional dependence of the momentum scattering rates, because in an anisotropic band and /or under the influence of anisotropic scattering mechanisms, the rate at which the carrier's momentum relaxes, will depend on the carrier's initial momentum and the distribution of momenta of the final states.
The variables that appear in Eq. 6 are as follows: DADP, DODP, DIVS are the deformation potentials for the ADP, ODP, and IVS mechanisms.  is the mass density, vs the sound velocity,  the dominant frequency of optical phonons, considered as constant over the whole reciprocal unit cell, which has been validated to be a satisfactory approximation, [22] and N is the phonon Bose-Einstein statistical distribution; 0 the vacuum dielectric constant, ks and k∞ the static and high frequency relative permittivities, is the generalized screening length with EF being the Fermi level and the carrier density [29,30]. c is the volume of the primitive cell, x the fraction of one of the alloy elements, and EG the difference between the energy gap of the two constituent materials that form the alloy. The G function is the form factor of a hard sphere [35].
The scattering rates and the transport coefficients are computed along the orthogonal  EPW+Wannier interpolated e-ph calculations for this material. We have performed such comparison for Si, for both n-type and p-type, using a full set of deformation potentials that we have extracted, with very good agreementsee Fig. 5 below [25].
well as ADP and ODP, for a donor doping density n = of 1.1×10 20 cm -3 . Note that in this case we consider the effect of this IIS donor scattering on both the majority electrons and minority holes [29].

III. Ab initio extraction of deformation potentials for scattering
The main scattering parameters required by ElecTra as described above, are the deformation potentials for electron-phonon scattering. Electron-phonon (e-ph) scattering is a vital part of simulations for materials properties. Deformation potentials are based on a theory developed by Bardeen and Shockley [37], and provide a convenient way to treat e-ph scattering using the analytical expressions in the previous section [30,38]. The deformation potential essentially describes the shift in the bands upon a change in the lattice caused by a perturbation from specific phonon modes, the ones that dominate the overall process.
Deformation potential theory is instrumental for the calculation of low-field mobility [35,39]. Traditionally, e-ph scattering is employed within transport methods such as the Boltzmann transport equation (BTE) [40,41], Monte Carlo [42,43], Landauer-Buttiker [44], etc.. It is also routinely used in semiconductor device transport simulators and high-field calculations in such devices, still with adequate accuracy [45,46,47]. The analytical scattering rates that result, are convenient when the e-ph scattering needs to be combined with other scattering mechanisms, such as for nanostructured materials [48], or highly doped materials and alloys for which ionized impurity scattering [18] and alloy scattering are important (such as transistor devices and thermoelectric materials [29]). The use of deformation potentials can allow for the flexibility and computational robustness that this process requires.
However, deformation potential values are only known for common semiconductor materials, and those are mostly extracted from experiment. For the arbitrary complex bandstructure materials, these are not known and need to be derived from first principles. Such approaches, however, are computationally extremely expensive, starting from their DFPT part. This can be typically performed on coarse meshes, but it is then interpolated using Wannier methods, because a dense electronic and phononic mesh discretization is required [56]. This leads to a large number of possible combinations of electronic and phononic states in the calculation of the e-ph interaction. For example, the use of interpolated 50×50×50 k-mesh and 25×25×25 q-mesh (typically half the size of the k-mesh), results in billions of possible matrix elements (the total number is given by the product of the two meshes). In fact, even larger meshes can be needed for mobility convergent calculations as shown in Ref. [58]. Typically, the symmetry of the unit cell is considered, which would reduce this number down to several million depending on the material specifics. Although these matrix elements are not calculated, but interpolated from the original coarse mesh onto the fine e-ph Wannier mesh [57], where they are subsequently used in transport calculations, this step is still computationally expensive. It can be even more expensive than the original DFPT/DFT on the coarse mesh (see Table 1 and discussions later on for a comparison example).
Here we describe a first-principles framework to extract deformation potentials in complex band materials based on density-functional theory (DFT) and density-functional perturbation theory (DFPT). We show that with the calculation of a reduced set of matrix elements and the formation of deformation potential scattering rate expressions rather than by using the full ab initio e-ph calculation, we can reduce the computational cost significantly, while not jeopardizing the accuracy [25]. We first compute the electronic band structures, phonon dispersion relations, and electron-phonon matrix elements to extract deformation potentials for acoustic and optical phonons for all possible processes.
The matrix elements clearly show the separation between intra-and inter-valley scattering and quantify the strength of the scattering events. The deformation potentials are extracted and then be used within the BTE as described in the previous section.
The process is as follows: Using ab initio self-consistent field DFT calculations we extract the electronic structures of the material of interest [58,59]. We employ density-functional perturbation theory (DFPT) [49,50]  The DADP (for LA and TA modes) computed in this way is comparable to the deformation potentials derived from band shifts after applying stress along specific directions, which is a good indication for an effective deformation potential, Deff, as described recently in Ref. [21]. That method extracts the corresponding scattering rates based on Deff, which is then used as a 'holistic' electron-phonon scattering process for all states and valleys in the bandstructure. This is indeed a very efficient way to extract first order transport properties at a very low computational cost, making it very effective form materials screening studies. However, it does not capture the details of the underlying physical mechanisms correctly, which can render it incomplete for TEs in certain cases.
For example, it does not capture optical phonon scattering (important for nanostructuring), or inter-valley processes, crucial for band-alignment optimization, which is the main strategy for power factor improvements in complex materials. Our method provides a step forward in terms of accuracy and physical relevance. We distinguish between all acoustic and optical phonon scattering processes, and all intra-and inter-valley processes, accounting for all selection rules automatically as well. The price we pay is the need for the initial coarse mesh expensive DFPT calculations, but overall our method is less expensive compared to full EPW+Wannier due to the reduced number of matrix elements required in the calculation of the transport properties.  Fig. 5b participate. This information allows for deeper understanding of transport. In a similar manner, to capture this for the rest of the scattering processes, the form factors, Fk,k', can also be computed, similarly to the deformation potentials, only for a few points along high symmetry lines and between valleys. The entire analysis for Si, as well as the result for the mobility calculations which match experiment very well, are all described in detail in Ref. [25], and here we have only provided a summary of this study.  [63]. (e-f) The Si electron and hole mobilities computed with our deformation potential method and the BTE code ElecTra, compared to experimental data. Data are from Ref. [64,65,66,67] for electrons, and Ref. [38,66,67] for holes. Blue dots are the simulation results and black triangles the experimental data.

Simulation steps in our method
The overall simulation steps for the electron-phonon interactions are as follows: Up to here, these are the same steps that fully ab initio methods (such as EPW [24]) perform (from here on, those calculations compute >billion matrix elements, usually 100k k-points and 10k q-points). We continue with the process of extracting deformation potentials using DFPT, which provides the matrix elements, one q-point at a time: 7) Find the coordinates of the valleys in the BZ (for CB/VB) and select the extrema points. These will form the initial/final k-points for the matrix elements. The fact that we compute a small amount of matrix elements, makes this method much more computationally affordable compared to fully ab initio Wannier methods, still with accuracy close to those methods (see Ref. [25] for Si). It is on the other hand, it is more expensive compared to other, simpler, deformation potential methods such as the one in Ref. [21], which is designed for reduced computational cost. Thus, our method is a step closer to full ab initio accuracy, with some associated increase in computational cost.
Currently it requires some hands-on time to identify the k-vectors and q-vectors for which the deformation potentials are computed, but we are in the process of automating this step.
In Table 1 below we show a computational cost comparison for Si between our deformation potential method (indicated as DP) and DFT+Wannier using the EPW code [24]. The full calculation is presented in our previous work [25]. The calculation starts with the DFPT step, which is common to both methods, typically on a 6×6×6 q-mesh (1 st row).
Our DP method then computes some matrix elements and then performs Boltzmann transport (2 nd and 3 rd rows). The corresponding EPW calculation that we have performed (4 th row) require significantly more CPU hrs. Note that the 64×64×64 interpolated mesh used is typical for EPW calculations [56], and it can even be an underestimation as the mobility might not be converged for that mesh size. Thus, we save consistently on the generation of the matrix elements and their subsequent use within BTE, which is a cost comparable typically even larger than the coarse-mesh DFPT calculation cost itself. The transport part requires less than 10% and the full calculation less than 20% CPU hours, compared to the fully ab initio method.  [24]. Step

(electrons and holes)
ElecTra It is worth noting that the unit cell of many promising TEs is large and their crystal symmetry low, for full ab initio calculations with practical computation expense. For example, Mg3Sb2 has five atoms in its unit cell and crystallizes in the trigonal space group, with lower symmetry than cubic. To accurately compute mobility using DFPT is extremely difficult for 3D materials (for 2D is also challenging but feasible [23]). General DFT calculations for Mg3Sb2 will need ~20 times more CPU hours than for Si. The bandstructure is more complex and many more phonon modes participate (see Fig. 6a-b).
Si would need 5k CPU hrs for standard 64×64×64 k-and q-meshes as shown in Table 1 above. First principle e-ph Wannier calculations would need ~100k CPU hrs for Mg3Sb2 (much more time for lower symmetry and larger unit cell). Rather than computing matrix elements throughout the Brillouin zone, the extraction of deformation potentials only needs a limited number of computations (for Mg3Sb2~50), and will require only a few thousand CPU hrs. An example of deformation potential extraction is shown in Fig. 6c. The CBM consists of six equivalent valleys, which results to three different transition types from a certain valley to the other 5 (inset of Fig. 6c). This inter-valley transitions are caused by all phonon modes. In Fig. 6c we compute the deformation potentials caused by each phonon mode individually and for each of the three transition types. Their strength can be quantified, and the dominant ones in each type identified (green box in Fig. 6c). All, or only the dominant ones, can make it to the BTE code for the computation of the transport properties.
Finally, for the rest of the scattering mechanisms, parameters such as Δ G can be computed by calculating the difference between the energy gap of the two constituent materials that form the alloy; ∞ and are standard outputs of general DFPT codes; and parameters like and can be calculated or found in many databases. Fk,k' can also be computed for a few cases within and from valley to valley.

IV. Transport in nanostructured materials using Monte Carlo
The Monte Carlo (MC) computational formalism we use is described in detail in Ref. [27]. It focuses on the effect of nanostructuring, and for this reason it uses many 'shortcut' deviations from usual MC algorithms to facilitate the challenges that simulating nanostructured geometries present. Specifically, we focus on geometries with the so-called  For computational simplicity and speed reasons, we follow two directions in using the real material properties and ab initio information extracted in the previous sections: 1) We consider parabolic bands and map the complexity of the real material bandstructure onto parabolic bands. For this we extract a conductivity and DOS effective masses (mC and mDOS, respectively) for the real bandstructure. We also use as input the deformation potentials to form the scattering rate equations. 2) In the second method, we extract the DOS(E), band velocity, v(E), mean-free-path λ(E), and scattering times, τ(E), for the material of interest using the ElecTra BTE code (which already uses the deformation potentials).
We begin below by describing the approach we follow to extract the effective masses to be used in MC, which is presented in [29,72], and then we describe the MC algorithm we employ, with specific modifications and peculiarities compared to current algorithms, designed specifically for the numerical complexities introduced by nanostructuring.

Effective mass extraction:
We compute the mDOS as the effective mass of an isotropic parabolic band that gives the same carrier density as the actual band structure. We evaluate mC as the effective mass of an isotropic parabolic band which maps the average velocity of the band states weighted by their contribution to transport. For this, we employ a simple ballistic field effect transistor model, extract the average injection velocity in the sub-threshold regime, and map that velocity to a parabolic band which provides the same injection velocity. The process is as follows (using the conduction band as an example): We consider the nondegenerate regime, in which the carrier density can be expressed as: where E0 is the energy of the band edge, EF is the Fermi level, kB is the Boltzmann constant, T is the temperature, and NC is the effective density of states calculated as: For a generic numerical band structure, the carrier density can be calculated as: where the sum is over all k-points and bands in the first Brillouin zone, ( , ) is the Fermi-Dirac distribution, and is the volume element in space, which usually depends only on the mesh. Then the mDOS can be obtained by combining Equations (7)(8)(9).
The conductivity effective mass in a specific direction, mC,i, is calculated from the injection velocity vinj of the carriers in the sub-threshold regime of a field effect transistor (FET): where the injection velocity is extracted as the FET current divided by the injected charge at the source contact: where EF,S is the source Fermi level, and IFET is the FET current density, which is computed from the difference between the source and drain currents. We assume that the drain current is negligible under high drain voltage, thus IFET can be given by: The vinj is identified and the conductivity effective mass is extracted from Equation (10) along the three Cartesian directions. Thus, three conductivity effective masses are extracted, and the final conductivity effective mass is calculated by averaging as:

Monte Carlo method for nanostructures
Monte Carlo simulations involve the ray-tracing of particle trajectories rather than the direct solution of partial differential equations. These particles are where it is ray-traced to either of the contacts, iii) does not require the application of a driving force. The method we present provides the same accuracy as common methods, but with a significantly reduced computational cost. It has many differences and peculiarities compared to standard MC methods, to tackle the issue of computational complexity. Details are provided in [27], and in this paper we provide the essentials and focus on the coupling and use of the various ab initio parameters within MC.
We use the incident flux (single-particle) approach, where the electrons are initialized at the domain boundaries one-by-one and propagate until they exit at either boundary (propagated to the other side or back-scattered as shown in Fig. 7d). Here we consider a two dimensional domain for numerical simplicity. Regarding the large computation time associated with the ray-tracing of low energy electrons, we consider a mean-free-path (mfp) approach, rather than the picking of random free-flight time and the self-scattering approach as is common practice [30,77]. We compute the total scattering rate of the particle using the deformation potential expressions described earlier, and using its bandstructure velocity we calculate its mean-free-path, λ(E). The particle propagates one mfp at a time (as in Fig. 7d), and then undergoes (enforced) scattering. In the case of acoustic phonon scattering for 3D parabolic bands, for example, the mean-free-path is where, tr(E) is the ToF for a single electron, and Nr(E) is the number of electrons that make it to the right end. We chose to keep the energy dependence since we only deal with elastic transport conditions for the example we present here. Then the averaged <ToF(E)> is used to calculate the flux per simulated electron at each energy as: Using the flux per electron at a given energy, we can form the overall flux in energy by multiplying by the density of states (DOS), g(E), either in its parabolic form, or its tabulated form from ElecTra for the real complex band material, which essentially is proportional to the transport distribution function (TDF) of the analytical BTE as: This is simply because the product of flux with DOS provides essentially the flow of charge, which is directly related to conductivity in the same way the TDF determines the conductivity. The proportionality constant C in the equation of the TDF accounts for the super-electron charge that is typically used in MC (the fact that we only simulate a finite number of electrons), and geometrical factors related to the simulation in a finite 2D domain rather than an infinite 3D domain (and connects the conductance to conductivity).
We use C to map the MC simulated TDF to the TDF that can be derived and extracted from the analytical solution of the BTE (or the one extracted from ElecTra), for the case of pristine material alone. To extract C, we calculate the electrical conductivity as a function of the Fermi energy both in the linearized BTE and with MC by integrating the TDF times the Fermi derivative in energy. We then find the mapping factor CEF for the conductivity at each Fermi energy. This comes to be almost constant for all EF, so we take the average of those values as the overall final C. After this, it also turns out that the TDFs from analytical BTE and MC are almost identical (see Fig. 8a).
After obtaining the transport distribution function numerically from MC and for the MC TDF, Ξ ( ) by calibrating to the linearized BTE one, we can substitute Ξ ( ) in the place of the analytical Ξ( ), which is given by ( ) 2 ( ) ( ) in the usual BTE formulas below. The electron conductivity is then calculated as: and the Seebeck coefficient as: Further, the TE power factor and the electronic thermal conductivity are evaluated, respectively, as: The above transport coefficients are computed using MC initially for the pristine material configuration, where the calibration of the constant C takes place. The electronic conductivities from the analytical BTE and the simulated MC in this way, are again very similar, almost identical (see Fig. 8c). The idea is that once this is calibrated, we can perform MC simulations for complex nanostructured domains without further calibration and benefit from the robustness of the single-flux injection method.
Thus, by using the MC extracted Ξ ( ), and by multiplying it by the derivative of the Fermi distribution function with respect to energy (i.e. = ) for the electrical conductivity (see Fig. 8b), and with respect to temperature (i.e. = − × ) for the Seebeck coefficient calculations, we essentially account for linear response. In this way we effectively eliminate: i) the need for two counter propagating simulation fluxes (now this is captured by the derivatives), and ii) the need for an application of a driving force, either a voltage difference or a temperature difference. Thus, we avoid the peculiar situation where a small enough potential difference window does not provide enough statistics, while a large enough could make the simulation deviate from linear response or from the range of voltages that TE materials utilize. Note also that the acquisition of adequate statistics is even more difficult in the case of a temperature gradient for the Seebeck coefficient in common bi-directional flux methods. We do not only need to differentiate between the right and left going fluxes, but also the ones which flow above and below the Fermi level.
We now consider the effect of nanostructuring. We simulate the geometries shown in Fig. 7 above and compute the electrical conductivity and power factor, shown in Fig. 9a and Fig. 9b, respectively. We show the conductivity as a function of the Fermi level, EF,

V. Conclusions
In this paper, we provided a review of a computational methodology which allows the extraction of the electronic and thermoelectric coefficients of complex nanostructured materials. The method merges three distinctive parts: i) a Boltzmann Transport Equation solver that we have recently developed, named ElecTra, which takes as input a complex bandstructure and scattering parameters, and provides transport coefficients and relevant quantities; ii) an ab initio methodology to extract the scattering parameters needed by ElecTra, and in particular deformation potentials; and finally iii) a Monte Carlo simulator which can utilize transport quantities from the ab initio calculations and the BTE transport above, with many features that deviate from common methods, specifically designed to provide robust computation for nanostructured simulation domains. In fact, each of the three methods we describe utilize unconventional routes that allow for faster and more robust calculations, without compromising accuracy at a significant degree. We believe this multi-physics framework, but also each of the simulators individually, are truly enabling, and can contribute significantly in exploring the electronic and thermoelectric properties of pristine and nanostructured materials in an efficient manner.