www.nature.com/npjcompumats ARTICLE OPEN The elphbolt ab initio solver for the coupled electron-phonon Boltzmann transport equations 1✉ 1✉ Nakib H. Protik , Chunhua Li2, Miguel Pruneda1, David Broido 2 and Pablo Ordejón 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. npj Computational Materials (2022)8:28 ; https://doi.org/10.1038/s41524-022-00710-0 1234567890():,; INTRODUCTION Frederikse13 and Geballe and Hull14,15. Around the same time, Ab initio computation of the transport properties of materials an influential theory was developed by Conyers Herring16 to allows a deeper understanding and probing of the rich underlying explain the experimental findings. Herring’s theory, which physics. It is also crucial for the predictive designing of materials contains several free-parameters, partially decouples the e and for industrial applications. The advances in density functional the ph BTEs, retaining an approximate term to account for the theory (DFT)1,2 and density functional perturbation theory (DFPT)3 drag effect. From a physical point of view, however, Herring’s have enabled accurate electron (e) and phonon (ph) band theory violates a deep, thermodynamic connection between the structure calculations. Furthermore, various freely available codes Seebeck and the Peltier effects known as the Kelvin–Onsager exist that allow fast computation of the ph–ph4 and e–ph5,6 relationship17. Semi-analytical work on 2D systems was carried out interactions. On the transport side, significant progress has been in 1987 by Cantrell, Butcher, and coworkers18,19. More recently, in made that allows computations of the mode (band/branch and 2014, a semi-analytical model with ab initio fitted parameters and wave vector) resolved lifetimes of electrons and phonons. The a partially decoupled framework was employed by Mahan, phonon thermal and the electronic charge conductivity can now Lindsay, and Broido to calculate the drag effect on the thermo- be calculated ab initio in the relaxation time approximation power of silicon20. Soon after, fully ab initio drag calculations with (RTA) or, going one step further, from a full solution of the partially decoupled solutions of the e and ph BTEs were devised corresponding single-species Boltzmann transport equation by Zhou et. al. in 2015 for silicon21. The code used from that work (BTE)4,6–9. In an interacting e–ph gas, however, the transport of was released in 202022. Ab initio solutions to partially decoupled e the two systems is intimately connected, and a unified, coupled and ph BTEs were also developed by Fiorentini and Bonini in 2016 e–ph transport may, instead, be needed for an accurate for silicon23 and Macheda and Bonini24 in 2018 for diamond. description of the underlying physics. The idea of a complete separation of the e and ph BTEs dates Finally, in 2020, a fully coupled e–ph BTEs solution was devised back to the 1930s and is known as Bloch’s Assumption10. A and applied to gallium arsenide using model e–ph interactions by paradoxical consequence of this assumption is that when one Protik and Broido25. In the same year, that method was improved solves the e (ph) BTE, the phonon (electron) system is taken to to include fully ab initio e–ph interactions to calculate the drag remain in equilibrium. This framework was promptly questioned effect in silicon carbide by Protik and Kozinsky26. by Peierls11, who argued that there must exist a momentum- Here we present elphbolt (short for electron–phonon mixing between the electrons and the phonons, which would, in Boltzmann transport), a code that features major improvements turn, cause both systems to move under the influence of a driving over the methods given in refs. 25,26. Moreover, we release the field. A theory of the coupled e–ph transport was developed by code as Free/Libre software under the GNU General Public Gurevich in 194612, which allowed calculation of the effect of the License version 3, bringing the capabilities for calculating the non-equilibrium electrons (phonons) on the phonons (electrons). e–ph drag physics via an ab initio solution of the fully coupled The first is called the electron drag and the latter, phonon drag. In e–ph BTEs, almost a century after Peierls’ conception of the idea, an interacting electron–phonon gas, however, the two effects are to the broader transport physics community. Our code is hosted inseparable and a mutual electron–phonon drag effect occurs. A on github27. decade later the first experimental evidence for the drag effect on Using an ab initio and Kelvin–Onsager relationship conserving the thermopower of germanium and silicon was found by solution of the coupled e–ph BTEs, elphbolt gives access to the: 1 Catalan Institute of Nanoscience and Nanotechnology (ICN2), CSIC and BIST, Campus Bellaterra, 8193 Barcelona, Spain. 2Department of Physics, Boston College, Chestnut Hill, Boston, MA 02467, USA. ✉email: nakib.haider.protik@gmail.com; pablo.ordejon@icn2.cat Published in partnership with the Shanghai Institute of Ceramics of the Chinese Academy of Sciences N.H. Protik et al. 2 ● mode resolved phonon thermal conductivity; where e is the electronic charge, V is the primitive unit cell volume, ● mode resolved phonon Peltier coefficient; τ and τ 0 label the basis atoms with masses mτ and mτ 0 , Z* is the ● mode resolved electronic charge conductivity; Born effective charge tensor, ε0 is the permittivity of free space, ● mode resolved electronic thermal conductivity; and ϵ∞ is the high-frequency dielectric tensor. This additional ● mode resolved electronic Seebeck and Peltier coefficients; term is only required for polar materials. ● and the effect of e–ph drag on all of the above quantities. We obtain the phonon branch energies ℏωsq and the This code is suitable for the study of 3d and 2d insulators, eigenstates jsqi by diagonalizing Dss0 ðqÞ with the unitary matrix semiconductors, semimetals, and metals. uq. Here ωsq is the phonon angular frequency. The phonon group In the sections below, we present the theory behind velocities are calculated using elphbolt, the implementation, and the outlook. 1 vsq ¼ hsqj∇q DðqÞjsqi: (5) 2ωsq RESULTS Equipped with these electronic and phononic quantities, we In this section, we give results for the calculated thermopower, move on to the calculation of the various interactions in the mobility, and thermal conductivity of n-doped silicon. First, the basic electron–phonon system. ingredients of the theory are described in detail. We present the We start with the e–ph interactions. The vertices (matrix elementary interactions considered in this work. This is followed by a elements) in Bloch space are given by30 discussion of the BTEs in the forms in which they are implemented, sffiffiffiffiffiffiffiffiffiffiffiffiffiffiffiffiffi _ 1 X X 0 0 0 y Fr;smn and the various types of solutions that elphbolt offers. Lastly, we kq ¼ gsmn expðik  Re þ iq  Rph Þ gsRemRphn Unn0 k0 Um 0 mk us0 sq þ gkq : 2mτ _ωsq Ne Nph R R s0 m 0 n0 present the transport coefficients that are obtained from the solution e ph of the BTEs along with a brief discussion of the Kelvin–Onsager (6) reciprocal relationship connecting the thermoelectric coefficients. where gsmnRe Rph are the e–ph matrix elements in Wannier space. The additional final term is the long-range correction due to the so- Electrons, phonons, and interactions called Fröhlich interaction33,34 required for polar materials. For these 1234567890():,; First, we need the electrons and phonons, and their interactions materials, the gRe Rph is constructed to be the short-range part of the calculated on arbitrarily fine wave vector meshes. These are full interaction. The Fröhlich term takes the following form: achieved by the Wannier interpolation techniques described in 2 P qffiffiffiffiffiffiffiffiffiffiffiffiffih i P ðq þ GÞ  Zτ  ξ τ;sq Fr;smn _ y refs. 5,28–30. We refer the readers to these seminal works for the gkq ¼ i εe0 V 2mτ _ωsq Uk Uk nm 0 ðq þ GÞ  ϵ1  ðq þ GÞ τ G≠q (7) details of the calculation of the Wannier functions and obtaining the Wannier representations of the various physical quantities ´ exp½iðq þ GÞ  rτ : starting from their Bloch representations. In particular, the where ξτ,sq = uτ,sq is the eigendisplacement of the basis atom τ due expressions for the Hamiltonian, dynamical matrix, and the e–ph to the phonon mode sq, rτ is the position of the atom τ in the matrix elements are given in refs. 5,30. Below we simply provide primitive unit cell, and the sum over the reciprocal lattice vectors G the expressions for the Wannier to Bloch transformations that are is performed using the Ewald sum technique by multiplying each computed within elphbolt. term of the summation by the factor exp½ðq þ GÞ  ϵ1  The Hamiltonian in the Bloch representation can be obtained by ðq þ GÞ=ð4αÞ, where the parameter α is set to 1 Bohr−2. the Fourier transformation: Recently, it has been shown that quadrupolar corrections may 1 X be important for an accurate description of the e-acoustic Hmn ðkÞ ¼ expðik  Re ÞHmn ðRe Þ; (1) Ne R phonon interactions35,36. Particularly, the electron mobility can e be strongly affected for those piezoelectric materials in which where m and n are band indices, k is an arbitrary wave vector, the transport-active band extremum is at the BZ center. H(Re) is the Hamiltonian in the Wannier representation, living on Currently, we do not have the capabilities for including the the real space spanned by {Re}, and Ne is the number of real space quadrupolar corrections. As such, care must be taken when unit cells. dealing with strongly polar materials such as cubic GaAs and Diagonalizing H(k), we obtain the electronic band energies ϵmk. wurtzite GaN that feature zone-centered band extrema. We plan The unitary matrix Uk digonalizing the Hamiltonian contains the to include support for quadrupolar corrections in a future eigenstates jmki. The band velocities are obtained from the release of the code. Hellmann–Feynman theorem31: In terms of the e–ph matrix elements described above, the 1 temperature-dependent transition rates of the electrons due to vmk ¼ hmkj∇k HðkÞjmki; (2) the phonon absorption (+) and emission (−) processes are _ given by37 where ℏ is the reduced Planck constant. 8 9 ( 0 ) Similarly, the dynamical matrix in the Bloch representation is