elphbolt: Ab Initio Electron–Phonon BTE Solver

• elphbolt is an ab initio framework for solving coupled electron–phonon Boltzmann transport equations, capturing mutual drag effects between electrons and phonons.
• It integrates DFT, DFPT, and Wannier interpolation to enforce Kelvin–Onsager reciprocity and accurately compute transport coefficients in metals, semiconductors, and 2D Dirac materials.
• The implementation leverages modern Fortran coarray parallelization to efficiently handle ultra-dense k- and q-point meshes for high-throughput transport calculations.

elphbolt is an ab initio computational framework designed for the efficient and thermodynamically consistent solution of the coupled electron–phonon Boltzmann transport equations (BTEs). Implemented in modern Fortran with explicit coarray parallelization, elphbolt enables predictive calculations of coupled charge, lattice, and thermoelectric transport coefficients—including both phonon drag (phonons driven out of equilibrium influencing electrons) and electron drag (nonequilibrium electrons affecting the phononic system). By leveraging density functional theory (DFT), density functional perturbation theory (DFPT), and maximally localized Wannier-function interpolation, elphbolt systematically captures full mutual electron–phonon drag with strict enforcement of Kelvin–Onsager reciprocity and is capable of resolving subtle transport phenomena in bulk metals, semiconductors, and 2D Dirac materials across a wide range of carrier concentrations and temperatures (Protik et al., 2021, Kazemian et al., 25 Aug 2025).

1. Theoretical Foundation: Coupled Electron–Phonon Boltzmann Equations

elphbolt formulates transport in response to external temperature gradients ($\nabla T$) and electric fields ($\mathbf{E}$) via the linearized, coupled BTEs for electron ($f_{m\mathbf{k}}$) and phonon ($n_{s\mathbf{q}}$) distribution functions:

• Electron: $$f_{m\mathbf{k}} \approx f_{m\mathbf{k}}^0[1+(1-f_{m\mathbf{k}}^0)\Psi_{m\mathbf{k}}]$$
• Phonon: $$n_{s\mathbf{q}} \approx n_{s\mathbf{q}}^0[1+(1+n_{s\mathbf{q}}^0)\Phi_{s\mathbf{q}}]$$ with $$\Psi_{m\mathbf{k}} = -\beta\nabla T\cdot I_{m\mathbf{k}} - \beta\mathbf{E}\cdot J_{m\mathbf{k}}$$, $$\Phi_{s\mathbf{q}} = -\beta\nabla T\cdot F_{s\mathbf{q}} - \beta\mathbf{E}\cdot G_{s\mathbf{q}}$$, and $\beta=1/(k_BT)$.

The coupled BTEs are decomposed into out-scattering (RTA, “0”), self-scattering ("S”), and drag (“D”) terms:

• $\nabla T$:
  • Electron: $\mathbf{E}$0
  • Phonon: $\mathbf{E}$1
• $\mathbf{E}$2:
  • Electron: $\mathbf{E}$3
  • Phonon: $\mathbf{E}$4

Collision integrals for electron–phonon, phonon–phonon, and phonon–electron scattering are explicitly constructed from first-principles electron–phonon matrix elements $\mathbf{E}$5, dynamical matrices $\mathbf{E}$6, and higher-order force constants (Protik et al., 2021, Kazemian et al., 25 Aug 2025). The drag terms, formulated as $\mathbf{E}$7 and $\mathbf{E}$8, link the nonequilibrium response of one subsystem to the other, allowing calculation of full mutual drag effects.

2. Numerical Strategy and Algorithmic Implementation

elphbolt implements an iterative "single-iterator" scheme indexed on the phonon BTE: for each outer phonon BTE iteration, the electron BTE is solved to convergence. This ensures reciprocal, self-consistent solutions where both electron and phonon distributions reach mutual equilibrium under nonequilibrium drive. The approach enforces Kelvin–Onsager reciprocity (e.g., $\mathbf{E}$9) at every stage.

Four solution regimes are supported:

• RTA (only “0” terms, no scattering-in)
• Self-scattering only (“S”)
• Fully coupled, dragless (“S” terms with off-diagonal coupling suppressed)
• Fully coupled, dragged (including “D” terms for mutual drag)

Delta-function energy conservation in scattering integrals is implemented via analytical tetrahedron or triangular integration (triangular required for 2D), eliminating artificial smearing. Convergence is determined by reduction of BTE residuals below strict numerical thresholds (Protik et al., 2021).

Fine sampling of the Brillouin zone is achieved by exploiting crystal point-group symmetries and restricting the active bands to transport-relevant windows around the Fermi energy (for electrons) or low-energy acoustic regions (for phonons), thus allowing practical treatment of ultra-dense $f_{m\mathbf{k}}$0-, $f_{m\mathbf{k}}$1-space grids requisite for accurate drag-resolved transport.

3. Workflow, Inputs, and Data Preparation

The workflow in elphbolt integrates the following ab initio quantities:

• Second- and third-order interatomic force constants (IFC2, IFC3) from DFPT and finite-difference methods (Quantum Espresso, ShengBTE).
• Wannier-interpolated tight-binding Hamiltonians, dynamical matrices, and electron–phonon matrix elements from EPW.
• Input control (input.nml): $f_{m\mathbf{k}}$2- and $f_{m\mathbf{k}}$3-mesh sizes, carrier doping, temperature, solution regime, convergence criteria.

The detailed computational steps are:

1. Generation of DFT, DFPT, and Wannier-interpolation data for the crystal.
2. Configuration of solver parameters via input files.
3. Initialization, symmetry reduction, and construction of transport-active mesh windows.
4. Calculation of temperature-dependent scattering rates using integration methods tailored to system dimensionality.
5. Iterative solution of coupled BTEs for electron and phonon deviations.
6. Output generation: mode- and total-resolved transport quantities, including electronic conductivity, Seebeck, Peltier, open/closed-circuit thermal conductivities, and spectral data for further analysis (Protik et al., 2021, Kazemian et al., 25 Aug 2025).

4. Parallelization, Computational Performance, and Usage

elphbolt exploits Fortran 2018's native coarray features to parallelize $f_{m\mathbf{k}}$4- and $f_{m\mathbf{k}}$5-point loops across distributed-memory nodes, distributing the computational workload by equal partitioning of irreducible Brillouin zone points. This balances memory and CPU utilization, supporting high-throughput calculations over massive sampling meshes.

For example, an $f_{m\mathbf{k}}$6-type Si system ($f_{m\mathbf{k}}$7 cm$f_{m\mathbf{k}}$8, $f_{m\mathbf{k}}$9 K, $n_{s\mathbf{q}}$0 $n_{s\mathbf{q}}$1-mesh, $n_{s\mathbf{q}}$2 $n_{s\mathbf{q}}$3-mesh) required $n_{s\mathbf{q}}$43000 CPU-hours for full electron–phonon BTE convergence; key components of the runtime stem from e-ph and ph–ph vertex construction and iterative BTE solving. Six outer phonon-BTE iterations typically yield full drag convergence (Protik et al., 2021).

Best practices include $n_{s\mathbf{q}}$5- and $n_{s\mathbf{q}}$6-mesh convergence to at least $n_{s\mathbf{q}}$7 and $n_{s\mathbf{q}}$8 for Si, with command-line driven coarray runs and monitoring for BTE residual thresholds below $n_{s\mathbf{q}}$9.

5. Transport Properties and Physical Insights

elphbolt computes, from the converged deviation functions:

• Electron conductivity $$f_{m\mathbf{k}} \approx f_{m\mathbf{k}}^0[1+(1-f_{m\mathbf{k}}^0)\Psi_{m\mathbf{k}}]$$0 and Seebeck $$f_{m\mathbf{k}} \approx f_{m\mathbf{k}}^0[1+(1-f_{m\mathbf{k}}^0)\Psi_{m\mathbf{k}}]$$1:

$$f_{m\mathbf{k}} \approx f_{m\mathbf{k}}^0[1+(1-f_{m\mathbf{k}}^0)\Psi_{m\mathbf{k}}]$$2

• Electronic Peltier coefficient $$f_{m\mathbf{k}} \approx f_{m\mathbf{k}}^0[1+(1-f_{m\mathbf{k}}^0)\Psi_{m\mathbf{k}}]$$3
• Lattice (phonon) thermal conductivity $$f_{m\mathbf{k}} \approx f_{m\mathbf{k}}^0[1+(1-f_{m\mathbf{k}}^0)\Psi_{m\mathbf{k}}]$$4:

$$f_{m\mathbf{k}} \approx f_{m\mathbf{k}}^0[1+(1-f_{m\mathbf{k}}^0)\Psi_{m\mathbf{k}}]$$5

• Phonon Peltier and total Seebeck/Peltier via Onsager consistency

The framework rigorously captures phenomena such as:

• Phonon drag peaks at low $$f_{m\mathbf{k}} \approx f_{m\mathbf{k}}^0[1+(1-f_{m\mathbf{k}}^0)\Psi_{m\mathbf{k}}]$$6 in the Seebeck coefficient
• Mobility enhancement in modulation-doped configurations
• Crossover from phonon- to electron-dominated heat transport as a function of doping and temperature
• Strong, nontrivial modifications to $$f_{m\mathbf{k}} \approx f_{m\mathbf{k}}^0[1+(1-f_{m\mathbf{k}}^0)\Psi_{m\mathbf{k}}]$$7, $$f_{m\mathbf{k}} \approx f_{m\mathbf{k}}^0[1+(1-f_{m\mathbf{k}}^0)\Psi_{m\mathbf{k}}]$$8, and $$f_{m\mathbf{k}} \approx f_{m\mathbf{k}}^0[1+(1-f_{m\mathbf{k}}^0)\Psi_{m\mathbf{k}}]$$9 due to mutual drag (Protik et al., 2021, Kazemian et al., 25 Aug 2025).

6. Applications and Physical Case Studies

Applications include:

• Bulk metals: Quantitative reproduction of ab initio results showing phonons carry up to 40% of thermal conductivity when electron–phonon scattering is included.
• Semiconductors: Prediction and interpretation of thermal conductivity suppression and subtle phonon drag features across broad doping and temperature ranges; validation versus time-domain thermoreflectance (TDTR) and isotope engineering data.
• 2D Dirac crystals (e.g., graphene): Resolution of mode-specific heat transport, demonstrating underestimation of $$n_{s\mathbf{q}} \approx n_{s\mathbf{q}}^0[1+(1+n_{s\mathbf{q}}^0)\Phi_{s\mathbf{q}}]$$0 in RTA, and highlighting the dominance of the flexural ZA branch at 300 K. The treatment of higher-order four-particle electron–phonon and phonon–phonon processes, as required for accurate transport at low Fermi energies, is outlined as an advanced extension (Kazemian et al., 25 Aug 2025).

7. Limitations, Extensions, and Open Challenges

elphbolt currently omits electron–electron and electron–plasmon scattering, four-phonon interactions, dynamical quadrupolar corrections, cluster-defect scattering, and magnetotransport. Planned extensions target rigorous treatment of these effects, expansion to correlated systems via DFT+U, hybrid functionals or GW/DMFT, and the inclusion of mode-to-mode, energy- and field-dependent Peierls–Boltzmann equations with higher-order (four-particle) processes and dynamically screened coupling for polar and Dirac systems.

Additional open challenges include interfacing with ultrafast, time-dependent transport solvers for non-equilibrium pump–probe scenarios, and embedding of ab initio mode-resolved lifetimes into multiscale, device-level simulations for next-generation electronic, thermoelectric, and phononic functionalities (Kazemian et al., 25 Aug 2025).