The elphbolt ab initio solver for the coupled electron-phonon Boltzmann transport equations

elphbolt is a modern Fortran (2018 standard) code for efficiently solving the coupled electron–phonon Boltzmann transport equations from first principles. Using results from density functional and density functional perturbation theory as inputs, it can calculate the effect of the non-equilibrium phonons on the electronic transport (phonon drag) and non-equilibrium electrons on the phononic transport (electron drag) in a fully self-consistent manner and obeying the constraints mandated by thermodynamics. It can calculate the lattice, charge, and thermoelectric transport coefficients for the temperature gradient and electric fields, and the effect of the mutual electron–phonon drag on these transport properties. The code fully exploits the symmetries of the crystal and the transport-active window to allow the sampling of extremely fine electron and phonon wave vector meshes required for accurately capturing the drag phenomena. The coarray feature of modern Fortran, which offers native and convenient support for parallelization, is utilized. The code is compact, readable, well-documented, and extensible by design.


Introduction
Ab initio computation of the transport properties of materials allows a deeper understanding and probing of the rich underlying physics.It is also crucial for the predictive designing of materials for industrial applications.The advances in density functional theory (DFT) [1,2] and density functional perturbation theory (DFPT) [3] have enabled accurate electron (e) and phonon (ph) band structure calculations.Furthermore, various freely available codes exist that allow fast computation of ph-ph [4] and e-ph [5,6] interactions.On the transport side, significant progress has been made that allows computations of the mode (band/branch and wave vector) resolved lifetimes of electrons and phonons.The phonon thermal and the electronic charge conductivity can now be calculated ab initio in the relaxation time approximation (RTA) or, going one step further, from a full solution of the corresponding single-species Boltzmann transport equation (BTE) [7,4,8,9,6].In an interacting e-ph gas, however, the transport of the two systems is intimately connected, and a unified, coupled e-ph transport may, instead, be needed for an accurate description of the underlying physics.
The idea of a complete separation of the e and ph BTEs dates back to the 1930's and is known as Bloch's Assumption [10].A paradoxical consequence of this assumption is that when one solves the e (ph) BTE, the phonon (electron) system is taken to remain in equilibrium.This framework was promptly questioned by Peierls [11], who argued that there must exist a momentum-mixing between the electrons and the phonons, which would, in turn, cause both systems to move under the influence of a driving field.A theory of the coupled e-ph transport was developed by Gurevich in 1946 [12], which allowed calculation of the effect of the non-equilibrium electrons (phonons) on the phonons (electrons).The first is called the electron drag and the latter, phonon drag.In an interacting electron-phonon gas, however, the two effects are inseparable and a mutual electron-phonon drag effect occurs.A decade later the first experimental evidence for the drag effect on the thermopower of germanium and silicon was found by Frederikse [13] and Geballe and Hull [14,15].Around the same time, an influential theory was developed by Conyers Herring [16] to explain the experimental findings.Herring's theory, which contains several free-parameters, partially decouples the e and the ph BTEs, retaining an approximate term to account for the drag effect.From a physical point of view, however, Herring's theory violates a deep, thermodynamic connection between the Seebeck and the Peltier effects knows as the Kelvin-Onsager relationship [17].Semi-analytical work on 2D systems was carried out in 1987 by Cantrell, Butcher, and coworkers [18,19].
More recently, in 2014, a semi-analytical model with ab initio fitted parameters and a partially decoupled framework was employed by Mahan, Lindsay, and Broido to calculate the drag effect on the thermopower of silicon [20].Soon after, fully ab initio drag calculations with partially decoupled solutions of the e and ph BTEs were devised by Zhou et. al. in 2015 for silicon [21].The code used from that work was released in 2020 [22].Ab initio solutions to partially decoupled e and ph BTEs were also developed by Fiorentini and Bonini in 2016 for silicon [23] and Macheda and Bonini [24] in 2018 for diamond.
Finally, in 2020, a fully coupled e-ph BTEs solution was devised and applied to gallium arsenide using model e-ph interactions by Protik and Broido [25].In the same year, that method was improved to include fully ab initio e-ph interactions to calculate the drag effect in silicon carbide by Protik and Kozinsky [26].
Here we present elphbolt (short for electron-phonon Boltzmann transport), a code that features major improvements over the methods given in Refs.[25,26].Moreover, we release the code as Free/Libre software under a GNU General Public License version 3, bringing the capabilities for calculating the e-ph drag physics via an ab initio solution of the fully coupled e-ph BTEs, almost a century after Peierls' conception of the idea, to the broader transport physics community.Our code is hosted on github [27].• mode resolved phonon thermal conductivity; • mode resolved phonon Peltier coefficient; • mode resolved electronic charge conductivity; • mode resolved electronic thermal conductivity; • mode resolved electronic Seebeck and Peltier coefficients; • and the effect of e-ph drag on all of the above quantities.This code is suitable for the study of 3d and 2d insulators, semiconductors, semimetals, and metals.
In the sections below, we present the theory behind elphbolt, the implementation, and outlook.

Results
In this section, we give results for the calculated thermopower, mobility, and thermal conductivity of n-doped silicon.First, the basic ingredients of the theory are described in detail.We present the elementary interactions considered in this work.This is followed by a discussion of the BTEs in the forms in which they are implemented, and the various types of solutions that elphbolt offers.Lastly, we present the transport coefficients that are obtained from the solution of the BTEs along with a brief discussion of the Kelvin-Onsager reciprocal relationship connecting the thermoelectric coefficients.

Electrons, phonons, and interactions
First, we need the electrons and phonons, and their interactions calculated on arbitrarily fine wave vector meshes.These are achieved by the Wannier interpolation techniques described in Refs.[28,29,30,5].We refer the readers to these seminal works for the details of the calculation of the Wannier functions and obtaining the Wannier representations of the various physical quantities starting from their Bloch representations.
In particular, the expressions for the Hamiltonian, dynamical matrix, and the e-ph matrix elements are given in Refs.[30,5].Below we simply provide the expressions for the Wannier to Bloch transformations that are computed within elphbolt.
The Hamiltonian in the Bloch representation can be obtained by the Fourier transformation: where m and n are band indices, k is an arbitrary wave vector, H(R e ) is the Hamiltonian in the Wannier representation, living on the real space spanned by {R e }, and N e is the number of real space unit cells.
Diagonalizing H(k), we obtain the electronic band energies mk .The unitary matrix U k digonalizing the Hamiltonian contains the eigenstates |mk .The band velocities are obtained from the Hellmann-Feynman theorem [31]: where is the reduced Planck constant.
Similarly, the dynamical matrix in the Bloch representation is obtained from its Wannier representation using where s and s denote the phonon branches, q is an arbitrary wave vector, D ss (R ph ) is the dynamical matrix in the Wannier representation, and R ph locates each unit cell in real space containing N ph cells.The first term is short-ranged, and the term D NAC ss is the long-range (the, so called, non-analytic) correction due to the dipole-dipole interaction given by [32] where e is the electronic charge, V is the primitive unit cell volume, τ and τ label the basis atoms with masses m τ and m τ , Z * is the Born effective charge tensor, ε 0 is the permittivity of free space, and ∞ is the high-frequency dielectric tensor.This additional term is only required for polar materials.
We obtain the phonon branch energies ω sq and the eigenstates |sq by diagonalizing D ss (q) with the unitary matrix u q .Here ω sq is the phonon angular frequency.The phonon group velocities are calculated using Equipped with these electronic and phononic quantities, we move on to the calculation of the various interactions in the electron-phonon system.
We start with the e-ph interactions.The vertices (matrix elements) in Bloch space are given by [30] g smn kq = 2m τ ω sq where g smn ReR ph are the e-ph matrix elements in Wannier space.The additional final term is the long-range correction due to the so-called Fröhlich interaction [33,34] required for polar materials.For these materials, the g ReR ph is constructed to be the short-range part of the full interaction.The Fröhlich term takes the following form: where ξ τ,sq = u τ,sq is the eigendisplacement of the basis atom τ due to the phonon mode sq, r τ is the position of the atom τ in the primitive unit cell, and the sum over G is performed using the Ewald sum technique by multiplying each term of the summation by where the parameter α is set to 1 Bohr −2 .
Recently, it has been shown that quadrupolar corrections may be important for an accurate description of the e-acoustic phonon interactions [35,36].Particularly, the electronic mobility can be strongly affected for those piezoelectric materials in which the transport active band extremum is at the BZ center.Currently, we do not have the capabilities for including the quadrupolar corrections.As such, care must be taken when dealing with strongly polar materials such as cubic GaAs and wurtzite GaN that feature zone-centered band extrema.We plan to include support for quadrupolar corrections in a future release of the code.
In terms of the e-ph matrix elements described above, the temperature dependent transition rates of the electrons due to the phonon absorption (+) and emission (-) processes are given by [37] where f 0 and n 0 is the equilibrium, i.e. the Fermi (Bose), distribution of the electron (phonon) gas, and the delta functions enforce the energy conservation in a scattering process.The notation k k |q above means that the initial electron wave vector is taken to be in the irreducible Brillouin zone (IBZ), the final electron wave vector is on the full first Brillouin zone (FBZ), the phonon wave vector mediating this transition is q = [k − k ], with the square brackets denoting a modulo operation with a reciprocal lattice vector.N k is the number of electronic wave vectors in the FBZ.
In elphbolt, we have the option of including electron-charged impurity (e-chimp) scattering to capture the effect of charged dopants on the transport properties.The e-chimp interaction is calculated within the first Born approximation, taking the impurity potential to have a static, Yukawa (screened Coulomb) form.The modulus squared of the interaction vertex is given by (generalized from Ref. [38]) where i = p, n denotes p-or n-type doping, n i is the concentration of the dopant, Z i is the ionization of the dopant, e is the absolute electronic charge, 0 is the zero frequency dielectric constant of the material, q is the wave vector magnitude measured from the nearest BZ center, and q TF is the Thomas-Fermi screening wave vector given by where d s is the spin degeneracy of the electronic state, and β ≡ (k B T ) −1 in terms of the Boltzmann constant k B and temperature T .
In terms of the matrix elements, the charged impurity mediated m k to nk transition rates are given by where we have multiplied the squared matrix elements with a transport factor corresponding to the in-scattering for elastic processes.Note that this term suppresses forward scattering.
At the present, we do not include e-e interactions.Typically, the e-e scattering does not strongly affect the transport in metals [37].However, in degenerate semiconductors, its effects might be strong [39].A rigorous treatment of the e-e interaction necessarily requires the calculation of the dynamical screening of the electron gas.Typically, this is done in the random phase approximation which introduces an additional Bosonic system -the plasmons.We wish to include e-e and the related e-plasmon scattering in a future release.
Next, we look at the interactions of the phonons.We start with the ph-e interaction.
The lowest-order process involves two electrons, and the transition probabilities are given by [37] Y ph-e s q |mknk = 2π N q g smn k q where N q is the number of phonon wave vectors in the FBZ and k = [k + q ].This describes the same + process given in Eq. ( 8).
The lowest-order ph-ph interaction vertices describing the q → q ± q processes are given by [4] V ±,ss s q q q = i jk αβγ where i, j, k label the atoms in the supercell, the primed sum indicates restriction to the central primitive unit cell, the Greek letters are Cartesian directions, e j,s q = exp(ir j • q )ξ j,s q is the eigendisplacement of atom j in the supercell due to the phonon mode s q , and Ψ ijk are the third order anharmonic force constants.
We will require only the modulus squared of the above expression.As such, we need only calculate the V − processes, since V +,ss s qq q In terms of the ph-ph vertices, the anharmonic phonon transition rates can be calculated as [4] W ± s q s q |s q = π 4N q V ±,ss s q q q 2 ω s q ω s q ω s q (14) We note that four-phonon interactions have been shown to be important for phonon transport in strongly anharmonic materials and those that feature anomalously weak three-phonon scattering rates for modes that dominate the transport at high temperatures; see Ref. [40] and the references therein.We do not currently have the functionality for calculating four-phonon interactions but plan to include it in a future release.
In elphbolt, we may also consider lowest-order ph-isotope and ph-substitution defect interactions within the first Born approximation and assuming the defect to be an on-site perturbation and low in concentration.This approximation can be described as the Tamura model [41].The transition rates for these two phonon processes are given by where x = {iso, subs} for isotope and substitution defect scattering, respectively, and is the mass variance parameter of the host atom τ , t denotes the type -isotope or substitution -of the guest atom, f tτ is the ratio of type-t guest to type-τ host, and M τ is the average on-site mass.
It is worth noting that the above description of the ph-defect interaction is simplistic and can fail for defects which cause extended bond perturbations.There exist advanced diagrammatic techniques that can handle such cases better [42,43,44].Such methods are planned for inclusion in a future release of elphbolt.

Coupled and decoupled BTEs
In elphbolt we may consider the electric (E) and temperature gradient (∇T ) fields as the drivers of the electron and the phonon currents.These applied fields cause the distribution functions of the electrons and the phonons, f mk and n sq , respectively, to deviate from their equilibrium forms.In the linear response regime, these are given by The deviation functions above of the electrons and the phonons, respectively, can be written as where I mk (F sq ) is the electron (phonon) response function that measures the deviation from equilibrium of the occupation of the electron (phonon) state mk (sq) due to applied ∇T field.Similarly, J mk (G sq ) is the electron (phonon) response to the E field.
The coupled electron and phonon BTEs have been given in numerous references e.g.[21,23,24,25,26,37].Here, we write them in the following form in terms of the response functions: ∇T field: In the equations above, the terms with the superscript 0 are the relaxation time approximation (RTA) terms.These involve a direct coupling to the applied field, which is why in the phonon equation to the E field, the RTA term is identically zero.Furthermore, the RTA only involves out-scattering processes.Truncating the BTEs up to the RTA term is equivalent to ignoring the in-scattering corrections and the drag effect.Next, the in-scattering corrections are given by the terms with the superscript S.These terms are functionals of the response function of the same species and, as such, they are called the self terms.The inclusion of these terms renders the BTEs (barring the one for G) into a set of decoupled equations which can be solved iteratively.Truncating the BTEs up to the self term, however, still ignores the drag effect.This is because, at this level of the approximation, the electron equations implicitly take phonons to remain in equilibrium, while the phonon equations take electrons to remain in equilibrium.Finally, the terms with the superscript D represent the drag terms.These are functionals of the other species.Inclusion of these terms allows an accurate description of the mutual drag effect in the interacting electron-phonon system.At this level, each pair of BTEs must be solved self-consistently.Furthermore, there exists a thermodynamic restriction known as the Kelvin-Onsager relationship [17] that reflects a deeper coupling of the two pairs of equations for the two fields.In elphbolt, we have devised a fast, iterative scheme that allows us to obtain the full solution of the Eqs.(19), while respecting the thermodynamic restrictions mandated by the Kelvin-Onsager relationship.This is discussed further in the Transport coefficients subsection.Now we provide the expressions for each term in the BTEs in Eq. ( 19).First we give the expressions for the electronic equations.The field coupling, RTA terms are given by where µ c is the chemical potential of the electronic system.The electron RTA scattering rates are given by The summation above of the independent scattering channels at the RTA level is known as the Matthiessen's Rule.
The self terms are given by Lastly, the drag terms are Similarly, for the phonons we have where the phonon RTA scattering rates are given by The self and drag terms are, respectively With these expressions, we solve the coupled e-ph BTEs, Eqs. ( 19), using an iterative procedure which is an improvement over the one used earlier in Refs.[25,26].Here we use a unified scheme, indexed by a single iterator, that of the ph BTE.In this approach, for a coupled BTEs solution, the e BTE is internally iterated to self-consistency for each iteration of the ph BTE.Thus, once the ph BTE has achieved self-consistency, the e BTE is guaranteed to do so also.This procedure is physically justified since the electron system is, in general, faster than the phonon system.Moreover, this scheme allows the strict enforcement of the Kelvin-Onsager relationship after each ph BTE iteration because the electron system is always brought to consistency with the current nonequilibrium phonon system.However, since the Kelvin-Onsager relationship mandates a thermodynamic constraint over the thermoelectric transport coefficients of both the electron and the phonon systems (see subsection Transport coefficients), in order to confirm that this relationship is satisfied, the coupled BTEs iteration must be performed for both the E and ∇T fields, simultaneously.Thus, all four BTEs in Eqs. ( 19) are tightly coupled.This approach is an improvement over the simple, dual iterator scheme used in Refs. [25] and [26], both conceptually and numerically.
We offer four different levels of solutions of the BTE for each species and field: RTA, partially decoupled, dragless full, and dragged full.Note that the RTA and dragless full solutions for the case of the phonons under the influence of the electric field are trivially zero.In our iterative approach, the computational time and space difference between the partially decoupled solution and the dragged full solution is not large, and both the RTA (iteration 0) and partially decoupled (iteration 1) solutions are provided en route to the dragged full solution.By default, elphbolt will carry out a set of decoupled, hence, dragless, ph and e BTEs immediately after the coupled BTEs are fully solved.This allows the users to compare the transport coefficients in all the above mentioned approximation levels after a single computational run.The users, however, also have the option of calculating only the decoupled e and ph BTEs (dragless full), if they wish to do so.

Transport coefficients
From the solutions of the response functions, we can obtain the following transport coefficient tensors.
From the electronic charge current, we get where σ is the electronic charge conductivity and S is the (Seebeck) thermopower.
From the electronic heat current, we get where α el is closely related to the electronic Peltier coefficient.The electronic component of the (Peltier) thermopower is given by Q el = α el (σT ) −1 .The tensor κ 0,el is the electronic thermal conductivity in the zero E field (closed circuit) condition.The electronic thermal conductivity in the open circuit condition, which is what can be measured in experiments, is given by κ el = κ 0,el − α el S.
Lastly, from the phonon heat current, we obtain where κ ph is the phonon thermal conductivity and α ph is related to the phonon Peltier coefficient.Specifically, the phonon component of the (Peltier) thermopower is given by The Kelvin-Onsager relationship unifies the Seebeck and the Peltier effects and mandates where α = α el + α ph [17].Note that this implies that the same thermopower Q ≡ S = Q el + Q ph is obtained by both the Seebeck and the Peltier pictures.
As such, the above relationship connects the ∇T field e BTE and the ∇T and E field ph BTEs.Since, for a given field, the BTE for one species is coupled to the one for the other via the drag term, the Kelvin-Onsager relationship effectively leads to a tight physical connection between all the four BTEs of the interacting e-ph system.However, during the course of the iterations, numerical issues may cause the system of equations to deviate from the Kelvin-Onsager relationship.To remedy this, corrective measures must be employed.Now, within the iterative scheme described in the previous subsection, it is straightforward to enforce the Kelvin-Onsager relationship.To do this, we split σS into a diffusion and a drag term and a draw parallel with the Peltier decomposition: Then, decomposing I as (band and wave vector indices dropped for brevity) and demanding σS diff [I diff ] = α el T −1 , we can arrive at the relation: The problem of enforcing the Kelvin-Onsager relationship then boils down to satisfying where λ is a small, corrective scalar which can be easily found using a bisection method.

Implementation
The coupled BTEs solver described above is implemented in Fortran 2018.This allows us to make use of the object-oriented programming (OOP) support and the built-in coarray functionality that provides concise, native syntax for parallelization.Specifically, we create the following 7 derived types: crystal, symmetry, numerics, electron, phonon, epw wannier, and bte dealing with the components of the problem that the names suggest.Each derived type contains its own data and procedures (functions and subroutines).Apart from these, there are separate modules for immutable parameters, We tried to strike a balance between the speed of development, execution, readability, and extensibility.Below we discuss the general workflow and the logical structure of the program.

Workflow and structure
In Fig. 2 we present the workflow of elphbolt.For a full calculation complete with the e-ph drag effect, we require the second-order interatomic force constants (IFC2s) from Quantum Espresso [45,46,47]; the third-order interatomic force constants (IFC3s) from thirdorder.py [4] which provides an interface with Quantum Espresso; and the Wannier space electronic Hamiltonian, dynamical matrix, e-ph matrix elements, and the real space cell maps and degeneracies from EPW [5,34,30].EPW internally interfaces with Quantum Espresso and Wannier90 [48].The generation of the IFC2s and IFC3s are also part of the ShengBTE workflow, and the users of that code will find that only one extra step -i.e. an EPW calculation -is required to generate all the input data for an eph coupled BTEs calculation with elphbolt.We have provided two modified EPW source files that allow the generation of some of the required Wannier space information.The user also needs to provide an input file called input.nml.The format of the input is described in details in the file README.org.
Next, we outline the logical flow of the program.The main program is in the file elphbolt.f90.
The program starts with creating objects of the seven derived types mentioned earlier.
Then, following a welcome message, it initializes the crystal object.This involves reading the information about the crystal and initializing the appropriate internal variables.
Following this, the numerics object is initialized which involves the reading in of the transport wave vector meshes, data output directory, BTE solution type, etc. Next, the symmetry object is initialized.At this step, the symmetries of the crystal and the BZ are calculated.Following this, the epw wannier object is initialized which involves reading in the Wannier space information.Next, we initialize the electron object.In this step, the information about the bands, transport active energy window, chemical potential, etc. are read in.The transport energy window restricted electronic IBZ is generated along with the bands, eigenstates, and velocities.Next we initialize the phonon object which involves the calculation of the phonon branches, eigenstates, and velocities.The next major step is the calculation of the interaction vertices.These zero temperature quantities are stored in the disk for later temperature and carrier concentration sweeps.
The temperature dependent transition rates are calculated and stored in the disk also for their reuse in the various types of BTE solutions that are available.The transition rates expressions include the energy conserving delta functions.We provide two methods for the evaluation of these delta functions: the triangular method [49,50] and analytic tetrahedron method [51].The first is the only option for 2d systems and is the default for 3d systems.Unlike the more commonly used Gaussian or Lorentzian methods, the triangular and tetrahedron methods do not have any additional smearing parameter, and, as such, the electron and phonon wave vector meshes are the only parameters to converge.
Finally, the bte object is used to solve the BTEs.The code produces a large amount of data for analysis, including the RTA scattering rates, band/branch resolved transport tensors, response functions at the RTA, partially decoupled, and fully iterated levels, among others.The users may also run the code in a post-processing mode to generate spectral transport coefficients, if required.Full descriptions of all the input options and the output data are provided in the README.orgfile.

Example: cubic silicon
As a demonstration of the code, we calculate the drag effect on the thermopower in cubic silicon.In addition, we include calculations of the mobility and thermal conductivity.
We consider first an n-type doped sample of low carrier density 2.75 × 10 14 cm −3 over a range of temperatures.This is close to the carrier concentration in the best sample (number 537) of Geballe and Hull's experimental work [15].We show in Fig. 3 a comparison of the calculated and the measured thermopower magnitudes.The experimental data is collected from Fig. 1 of Ref. [15].The black curve is the calculated total thermpower.In the Seebeck picture, this includes the phonon drag contribution.And in the Peltier picture, this is the sum of the electronic contribution and the phonon contribution, the latter being purely due to the electron drag effect.The Peltier picture also allows a clean separation of the electronic and the phonon contributions.This Peltier breakdown of the total thermopower is shown in Fig. 3.The solid blue curve is the electronic contribution to the thermopower.This contribution decreases slightly with decreasing temperature.The solid green curve is the phonon thermopower.This starts off as significantly lower than the electronic thermopower at 300 K, but quickly overtakes the latter around 175 K, before dominating by more than an order of magnitude at 50 K.The calculated total thermopower is in excellent agreement with the experimental measurements (red circles).Without the drag effect, the phonon contribution would be  trivially zero, and the electronic contribution would be an order of magnitude off from the experimental values at low temperatures.Furthermore, without drag, the temperature dependence of the thermopower would be spectacularly wrong.The explanation of the origin of the strong drag effect in silicon has previously been given in Refs.[21,23], and is not reproduced here.These earlier ab initio works captured the strong drag behavior using a partially decoupled solution of the e and ph BTEs.Nevertheless, they found excellent agreement with experimental measurements.Our calculations corroborate their finding that the partially decoupled solution is indeed sufficient to the capture the strong drag effect in silicon -the green and blue squares denote the partially decoupled solutions, and they coincide nearly perfectly with corresponding full drag curves.
For this case, the RTA and dragless full solutions of the e BTE also give essentially the  same result and are not shown on the plot to reduce clutter.However, this is not guaranteed to hold true for all materials, and, in general, a full drag solution of the e-ph BTEs should be used.Fig. 4 shows the calculated thermopower of an n-type sample with a high density of 2.7 × 10 19 cm −3 , matching that of 140 in Ref. [15].Good agreement with measured data is again obtained.The percentage difference of the calculated total thermopower from the experimental value near 50 (300) K is around 15 (1)%.Note that the calculated thermopower without the phonon component again significantly underestimates the measured thermopower across the full range of temperatures considered.
In Fig. 5, we present the calculated temperature dependent mobilities for a carrier concentration of 2.75 × 10 14 cm −3 .Excellent agreement is found over the entire temperature  range considered here with experimental measurements (shown in red circles) on similar low-doped samples [52].The effect of phonon drag on this quantity is negligible.In fact, the difference between the dragged full, RTA, dragless full, and partially decoupled solutions is very small, and, to reduce clutter, we do not show the latter two results on the plot.This corroborates the finding in Ref. [21].Fig. 6 shows the calculated temperature dependent mobilities for a degenerate carrier concentration of 2 × 10 19 cm −3 and the comparison to experimental data from samples with similar electron concentrations from Ref. [53].The calculated mobilities including e-chimp scattering (solid blue curve) are between a factor of 4 to 5 higher than the measured values (red circles).The RTA, dragless full, and partially decoupled solutions give nearly the same values as the dragged full solution and, to reduce clutter, we do not   show these points on this plot.The large discrepancy between the calculated and the measured mobilities could be due to a combination of multiple reasons.First, the measurements were done on compensated samples and the acceptor and donor compositions were not reported in Ref. [53].Our calculations have assumed that the mobile electron concentration is equal to the ionized donor density, while the acceptor density has been taken to be zero.In the actual samples, acceptors provide additional charge scattering centers thus lowering the mobility.Second, the e-chimp scattering employed in the calculation assumes a static, Thomas-Fermi screened Coulomb interaction treated in the Born approximation.It is known that the Thomas-Fermi model leads to a significant over-screening of the interaction in the degenerate limit [38].The validity of the Born approximation is also stretched in this limit and more sophisticated non-perturbative ap-proaches might be better suited.Lastly, the consideration of the e-plasmon scattering, which is currently not included in our calculation, has been shown to reduce the mobility of silicon significantly at high carrier concentrations [39].A more rigorous treatment of the e-chimp scattering along with e-e scattering is planned for a future version of elphbolt.
We also show in Fig. 6 the calculated mobility with the e-chimp interaction turned off (dashed blue curve).Here we envision that the mobile charge carriers are created in a region in which charged dopants do not exist, which could happen, for example, through modulation doping [55], or in GaN/AlGaN heterostructures [56].In this case, the phonon drag effect leads to a large increase in the mobility -at 300 (50) K, the drag gain of mobility is about a factor of 2.5 (50).Such high gains in the mobility can, potentially, be exploited along with those previously identified in the thermopower [21,23,57] to boost the thermoelectric figure-of-merit.Note that in this case the phonon drag effect is not fully captured by the partially decoupled solution which predicts about half the value given by the dragged full solution at 50 K.At this temperature, the RTA and dragless full solutions undershoot the value of the dragged full solution by about a factor of 5.
Finally, in Fig. 7 we compare the calculated phonon thermal conductivities against measurements on high purity samples of natural silicon.The measured values (red circles) are taken from Ref. [54].Remarkable agreement is found over the full temperature range considered for the low concentration case (solid blue curve).We also plot the results for a high concentration case in dashed green.The ph-e interactions cause a weaker suppression at high temperatures compared to that in the low temperature limit even at a high carrier concentration of 2 × 10 19 cm −3 .This happens because the ph-ph scattering rates dominate over the ph-e ones at high temperatures.This is consistent with the findings in Ref. [58].At low temperatues, the low energy acoustic phonons progressively contribute more to the thermal conductivity, while, at the same time, the ph-e scattering rates begin to dominate over the ph-ph ones for these modes.This results in a stronger suppression of the thermal conductivity at low temperatures for the high doped case.
The electron drag effect on the phonon thermal conductivity is found to be small in this material, as has been shown earlier in Refs.[21,23].This has also been shown to be true for gallium arsenide [25] and silicon carbide [26].The reason for this has been discussed in these references.The dragged full, dragless full, partially decoupled, and the RTA solutions all give nearly the same results and only the first type is shown on the plot.

Discussion
In this work, we discussed the theory and implementation of elphbolt -an efficient code for solving the coupled electron-phonon Boltzmann transport equations.This gives ab initio access to the thermal, charge, and thermoelectric transport properties in materials, and the effect of e-ph drag on them.The code is distributed as Free/Libre software under the GNU General Public License version 3 that gives the user the right to use, modify, and distribute the original and their modified versions of the software.The code combines the object-oriented and the procedural programming styles, has a clean, modular structure, features coarray parallelization, and is well-documented.This makes the code easy to extend.For future releases, we plan to include a more rigorous treatment of impurity scattering, electron-electron interactions, provisions for including four-phonon interactions [59], quadrupolar corrections, and magnetotransport.Table 1: Breakdown of the computational time needed for the various expensive parts of the code for silicon.An n-type carrier concentration of 2.75 × 10 14 cm −3 at 300K using a (50, 150) mesh set is considered here.

Figure 1 :
Figure 1: The logo of elphbolt.The red symbol is a phonon self-energy diagram and the yellow bolt signifies both the electron transport and the efficiency of the code.

Figure 2 :
Figure 2: Workflow of elphbolt showing the various input data required for a coupled e-ph BTEs calculation.

Figure 5 :
Figure 5: Temperature dependence of the mobility of silicon for an n-type carrier concentration of 2.75 × 10 14 cm −3 .The red circles are measurements on various different samples with carrier concentrations ranging from 3.5 × 10 13 to 1.4 × 10 14 cm −3 [52].

Figure 7 :
Figure 7: Temperature dependence of the phonon thermal conductivity of silicon for ntype carrier concentrations of 2.75 × 10 14 and 2 × 10 19 cm −3 .The red circles are measurements on high purity samples with natural isotopic mix [54].
task e-ph vertex ph-ph vertex e-ph BTEs e BTE ph BTE time (cpu-hours)