LinReTraCe: The Linear Response Transport Centre

Matthias Pickem 1, Emanuele Maggio 1, Jan M. Tomczak 1

1 Institute of Solid State Physics, TU Wien, Austria

## Abstract

We describe the ”Linear Response Transport Centre” (LinReTraCe), a package for
the simulation of transport properties of solids. LinReTraCe captures quantum
(in)coherence effects beyond semi-classical Boltzmann techniques, while
incurring similar numerical costs. The enabling algorithmic innovation is a
semi-analytical evaluation of Kubo formulae for resistivities and the
coefficients of Hall, Seebeck and Nernst. We detail the program’s
architecture, its interface and usage with electronic-structure packages such
as WIEN2k, VASP, and Wannier90, as well as versatile tight-binding settings.

###### Contents

  1. 1 Introduction
    1. 1.1 Methodological context
    2. 1.2 Transport coefficients
  2. 2 Transport equations and optical elements
    1. 2.1 Kubo kernel expressions
    2. 2.2 Boltzmann approximation
    3. 2.3 Quasi-particle renormalizations
    4. 2.4 Matrix elements
      1. 2.4.1 Dipole optical elements
      2. 2.4.2 Peierls approximation
  3. 3 Implementation
    1. 3.1 Requirements
    2. 3.2 Installation
    3. 3.3 Flowchart
    4. 3.4 Energy File
      1. 3.4.1 Density functional theory codes
      2. 3.4.2 Wannier90
      3. 3.4.3 Tight binding models
      4. 3.4.4 Generic interface
    5. 3.5 Scattering File
    6. 3.6 Configuring and running LinReTraCe
    7. 3.7 Output File
  4. 4 Technical details
    1. 4.1 Matrix elements on irreducible grids
    2. 4.2 Chemical potential search
      1. 4.2.1 Impurity levels
      2. 4.2.2 Homogeneous doping
    3. 4.3 Polygamma evaluation
    4. 4.4 Code scaling
  5. 5 Test cases
    1. 5.1 Silicon: Benchmark to BoltzTraP
    2. 5.2 One-band metal in two dimensions
    3. 5.3 Two-band insulator in three dimensions
      1. 5.3.1 Temperature scan
      2. 5.3.2 Chemical potential scan
    4. 5.4 Honeycomb lattice (tight-binding calculation)
    5. 5.5 FeAs2 (WIEN2k calculation)
    6. 5.6 Tl-doped PbTe (WIEN2k calculation)
  6. 6 Conclusion
  7. A Cheat sheet
  8. B Advanced LinReTraCe configuration
    1. B.1 General settings
    2. B.2 Temperature mode settings
    3. B.3 Scattering files
  9. C Dimensional analysis
    1. C.1 Internal units
    2. C.2 External units
    3. C.3 Lower dimensionality
      1. C.3.1 two dimensions
      2. C.3.2 one dimension
  10. D Interband - intraband limit

##  1 Introduction

Signatures of external electrical fields, magnetic fields, thermal gradients,
or combinations thereof can be used to monitor induced charge and energy
currents in materials. Analyzing these transport properties, the ability to
conduct charge, heat, or entropy, not only provides fundamental insight, but
also quantifies potential functionalities.

A key to the interpretation of measurements and to predict or optimize such
properties are transport simulations. Here, we describe the Linear Response
Transport Centre LinReTraCe (github.com/linretrace/linretrace), a software
package that facilitates the computation of a variety of transport
observables. The unique feature of LinReTraCe is the treatment of thermal and
lifetime broadening on an equal footing [1, 2, 3], while still incurring
numerical costs as low as semi-classical Boltzmann approaches in the
relaxation-time approximation. We exploit that linearizing the dynamics of
many-body renormalizations (self-energy) allows for a semi-analytical (instead
of numerical) evaluation of leading contributions in Kubo’s linear response
theory [4, 1, 2]. LinReTraCe’s principle input are electronic excitation
energies and associated quasi-particle weights and lifetimes, as well as
optical transition matrix elements. Being agnostic to the input’s origin,
LinReTraCe can be used in a variety of settings, including electronic
structures from tight-binding or Wannier projections[5], density functional
theory (DFT) [6], many-body perturbation theory[7, 8], dynamical mean-field
theory (DMFT) [9, 10], or approaches beyond [11, 12]. Scattering amplitudes
and many-body renormalizations can be phenomenological, extracted from
electronic self-energies (obtained, e.g., from DMFT), or could incorporate
results from electron-phonon codes [13, 14, 15, 16, 17, 18, 19]. In this
release we include interfaces to the DFT codes WIEN2k[20, 21] and VASP[22],
the band interpolation tool of BoltzTraP2 [23, 24], maximally localized
Wannier functions of Wannier90[25], as well as tools for general tight-binding
systems. We further provide templates for the implementation of interfaces to
other codes. With an emphasis on numerical accuracy and scalability,
LinReTraCe will be of value also for high-throughput studies[26, 27, 28, 29,
30, 31] and for tight-binding descriptions of very large unit-cells.

###  1.1 Methodological context

To highlight the merits of LinReTraCe, let us group previous packages for
electronic transport properties of solids into two categories: semi-classical
Boltzmann and Kubo linear response codes. Owing to their numerical efficiency
and ease of handling, Boltzmann codes [23, 24, 14, 13, 32, 33] have become
popular tools. Typically, they are used in conjunction with band-theory, where
excitation energies are sharply defined. The selection of carriers
participating in conduction is then solely determined from the thermal
broadening (activation) via the Fermi function. In the usually employed
relaxation-time approximation (RTA), electron scattering then only results in
amplitude-scaling prefactors for a given momentum k and state $n$, e.g., for
the conductivity
$\sigma\propto\sum_{\mathbf{k}n}\Gamma_{\mathbf{k}n}^{-1}\times\cdots$, cf.
Eq. (38). Often, the scattering rate $\Gamma$ is moreover assumed to be equal
for all states and momenta. Then, the lifetime $\tau=\frac{\hbar}{2\Gamma}$
scales the Boltzmann conductivity globally, while the Seebeck and Hall
coefficient become independent of $\tau$. With this assumption, insight into
electronic transport can be gained without explicit knowledge of scattering
amplitudes and their physical origin. As a consequence, RTA Boltzmann
transport kernels are relatively simple and the algorithmic complexity of,
e.g., BoltzTraP[23, 24] and BoltzWann[32] (and the difference between them),
lies in how they achieve convergence in the sampling of the Brillouin zone.

There are however circumstances, when the approximations inherent to band
theory and the semi-classical treatment of transport fail. Inadequacies of
band-theory for strongly correlated materials are well-documented: Examples
not only include Mott insulators, where correlation effects fully invalidate
the band-picture[34]. Also the electronic structure of correlated metals[35,
36, 37, 38] and correlated narrow-gap semiconductors [39, 40] are severely
altered and have to be accounted for with methodologies that include dynamical
renormalizations (self-energies). However, even if such many-body corrections
are captured on the level of electronic structure theory, plugging them into a
semi-classical transport methodology may still lead to severe pathologies:
Irrespective of the size of the scattering rate, the resistivity, the Seebeck,
and the Hall coefficient of a clean semiconductor diverge in the zero
temperature limit within Boltzmann’s relaxation time approximation. A
diverging activation law for the $T\rightarrow 0$ resistivity is physically
admissible—but it is never observed. A diverging Seebeck coefficient,
$|S(T\rightarrow 0)|\rightarrow\infty$, instead, violates the third law of
thermodynamics (there can be no entropy transport at $T=0$ for non-degenerate
ground-states).

A quantum mechanical description of transport using Kubo’s linear response
theory[41] overcomes these artefacts by correctly treating effects of finite
lifetimes (incoherence) of charge carriers[42, 1, 2]. To allow, beyond the
thermal, also for a lifetime broadening of excitations, Kubo formulas require
an integration over energies. This evaluation may become expensive for large
systems and is hazardous at low temperatures and when the scattering rate is
small. With LinReTraCe we conquer this bottleneck by performing frequency
integrations analytically instead of numerically. This step becomes possible
after linearizing the dynamics of many-body renormalizations (the self-energy
$\Sigma$), yielding the LinReTraCe input: the scattering rate $\Gamma$, the
quasi-particle weight $Z$, and possible static offsets $\Re\Sigma$. This
approximation is warranted if the self-energy varies slowly inside the narrow
energy-window probed by transport (a few $k_{B}T$). This is typically the case
for conventional and moderately correlated systems. In strongly correlated
metals, LinReTraCe investigations, while superior to Boltzmann studies, should
focus on low enough temperatures. For Mott insulators, with pole-like
structures in the self-energy, an evaluation of the full Kubo kernels is
required. Following common practice of such Kubo implementations[43, 44, 45,
46, 47, 48, 49, 38, 36, 37, 50], LinReTraCe also neglects particle-hole
scattering, so-called vertex corrections.111Vertex corrections mediated, e.g.,
by the Coulomb interaction (leading to, e.g., excitons or $\pi$-tons[51]) or
disorder scattering (effects of, e.g., weak localization[52]), that dress the
excited particle-hole pairs, are actively researched [53, 54, 55, 51, 56, 57],
but are often small in practice[56]. Formally, they vanish in the (DMFT) limit
of infinite dimensions[9], at least in the one-orbital case[58]. The ensuing
analytical transport functions are then not only numerically inexpensive and
stable, they also reveal valuable microscopic information: In particular they
show that the scattering is not a mere prefactor (scaling the amplitude of
conduction), but a relevant energy scale that engages in a complex interplay
with other energies of the system, e.g., the charge gap in a semiconductor[2].
In all, LinReTraCe combines the best of both (Boltzmann & Kubo) worlds: an
efficient and stable evaluation of transport observables that treats thermal
and lifetime broadening on an equal footing. For a detailed derivation of the
theory, see Refs. [2, 3].

###  1.2 Transport coefficients

Physical observables of transport processes are described by charge and heat
currents

$\displaystyle j_{e}^{\alpha}$
$\displaystyle=\mathcal{L}_{11}^{\alpha\beta}E^{\beta}-\frac{1}{T}\mathcal{L}_{12}^{\alpha\beta}\partial_{\beta}T$
(1) $\displaystyle j_{q}^{\alpha}$
$\displaystyle=-\mathcal{L}_{21}^{\alpha\beta}E^{\beta}-\frac{1}{T}\mathcal{L}_{22}^{\alpha\beta}\partial_{\beta}T,$
(2)

generated by an external electric field $\mathbf{E}$ and a temperature
gradient $\nabla T$ perturbing the system. This non-equilibrium state can be
described within linear response theory. There, the coupling constants (the
Onsager coefficients) become

$\displaystyle\mathcal{L}_{ab}^{\alpha\beta}$ $\displaystyle=\frac{\pi\hbar
e^{\left(4-a-b\right)}}{V}\sum_{\begin{subarray}{c}n,m\\\
\mathbf{k},\sigma\end{subarray}}\mathcal{K}_{ab}(\mathbf{k},n,m)M^{\alpha\beta}(\mathbf{k},n,m)$
(3) $\displaystyle\mathcal{L}_{ab}^{B,\alpha\beta\gamma}$
$\displaystyle=\frac{4\pi^{2}\hbar
e^{\left(5-a-b\right)}}{3V}\sum_{\begin{subarray}{c}n,m\\\
\mathbf{k},\sigma\end{subarray}}\mathcal{K}_{ab}^{B}(\mathbf{k},n,m)M^{B,\alpha\beta\gamma}(\mathbf{k},n,m)$
(4)

with $a,b\in\\{1,2\\}$ defining the coupling in the Cartesian
$\alpha,\beta,\gamma\in\\{x,y,z\\}$ direction and the summations running over
all possible band combinations $(n,m)$, momenta in the Brillouin zone
$\mathbf{k}\in\mathrm{BZ}$ and spins $\sigma$. The coefficients
$\mathcal{L}_{ab}^{B,\alpha\beta\gamma}$ are the linear correction terms in
the presence of a magnetic flux ($\mathbf{B}$) in the $\gamma$-direction.
Employing the Onsager-Casimir relation

$\mathcal{L}_{ab}^{\alpha\beta}(\mathbf{B})=\mathcal{L}_{ba}^{\beta\alpha}(-\mathbf{B})$
(5)

on the expansion

$\mathcal{L}_{ab}^{\alpha\beta}(\mathbf{B})=\mathcal{L}_{ab}^{\alpha\beta}+\mathcal{L}_{ab}^{B,\alpha\beta\gamma}B^{\gamma}+\mathcal{O}(\mathbf{B}^{2}),$
(6)

which we treat up to linear order in $\mathbf{B}$, one finds the connections

$\displaystyle\mathcal{L}_{ab}^{\alpha\beta}$
$\displaystyle=\mathcal{L}_{ba}^{\beta\alpha}$ (7)
$\displaystyle\mathcal{L}_{ab}^{B,\alpha\beta\gamma}$
$\displaystyle=-\mathcal{L}_{ba}^{B,\beta\alpha\gamma}$ (8)

when the $\mathbf{B}$-field is the only source of time-reversal symmetry
breaking.222In magnetized materials the Onsager relations take the form of
$\mathcal{L}_{ba}^{\beta\alpha}(H,M)=\mathcal{L}_{ab}^{\alpha\beta}(-H,-M)$
instead ($H$: magnetic field strength; $M$: magnetization density), leading to
the anomalous Hall effect and a linear magneto resistance [59, 60, 61]. The
Onsager coefficients combine the transport kernel functions, describing the
(in our case: free) particle-hole propagation of the excited system,

$\displaystyle\mathcal{K}_{ab}(\mathbf{k},n,m)$
$\displaystyle=\int_{-\infty}^{\infty}d\omega\;\omega^{(a+b-2)}\left(-\frac{\partial
f}{\partial\omega}\right)A_{\mathbf{k}n}(\omega)A_{\mathbf{k}m}(\omega)$ (9)
$\displaystyle\mathcal{K}_{ab}^{B}(\mathbf{k},n,m)$
$\displaystyle=\int_{-\infty}^{\infty}d\omega\;\omega^{(a+b-2)}\left(-\frac{\partial
f}{\partial\omega}\right)A^{2}_{\mathbf{k}n}(\omega)A_{\mathbf{k}m}(\omega)$
(10)

with the optical matrix elements $M^{(B),\alpha\beta(\gamma)}$, describing the
coupling of the external perturbation to the system. From them, the electrical
conductivity ($\sigma$), electrical resistivity ($\rho$), Peltier coefficient
($\Pi$), Seebeck coefficient ($S$), thermal conductivity ($\kappa$), Hall
conductivity ($\sigma^{B}$), Hall coefficient ($R_{H}$), Nernst coefficient
($\nu$)333We use the historical convention for the sign of the Nernst
coefficient[62]., Hall mobility ($\mu_{H}$) and its analogue, the thermal
mobility ($\mu_{T}$) [63] can be calculated. Note that, here, we limit the
transport kernels to symmetric contributions. Anti-symmetric terms leading to
anomalous transport from a non-trivial topological state [64, 60, 65, 66] are
neglected. For a detailed dimensionality analysis of all involved quantities,
see Appendix C.

$\displaystyle\sigma_{\alpha\beta}$
$\displaystyle=\mathcal{L}_{11}^{\alpha\beta}$ (11)
$\displaystyle\rho_{\alpha\beta}$
$\displaystyle=\left(\mathcal{L}_{11}^{-1}\right)^{\alpha\beta}$ (12)
$\displaystyle\Pi_{\alpha\beta}$
$\displaystyle=-\left(\mathcal{L}_{11}^{-1}\right)^{\alpha
i}\mathcal{L}_{12}^{i\beta}$ (13) $\displaystyle S_{\alpha\beta}$
$\displaystyle=-\frac{1}{T}\left(\mathcal{L}_{11}^{-1}\right)^{\alpha
i}\mathcal{L}_{12}^{i\beta}$ (14) $\displaystyle\kappa_{\alpha\beta}$
$\displaystyle=\frac{1}{T}\left[\mathcal{L}_{22}^{\alpha\beta}-\mathcal{L}_{12}^{\alpha
i}\left(\mathcal{L}_{11}^{-1}\right)^{ij}\mathcal{L}_{12}^{j\beta}\right]$
(15) $\displaystyle\sigma^{B}_{\alpha\beta\gamma}$
$\displaystyle=\mathcal{L}_{11}^{B,\alpha\beta\gamma}$ (16) $\displaystyle
R_{H,\alpha\beta\gamma}$
$\displaystyle=\left(\mathcal{L}_{11}^{-1}\right)^{\alpha
i}\mathcal{L}_{11}^{B,ij\gamma}\left(\mathcal{L}_{11}^{-1}\right)^{j\beta}$
(17) $\displaystyle\nu_{\alpha\beta\gamma}$
$\displaystyle=-\frac{1}{T}\left(\mathcal{L}_{11}^{-1}\right)^{\alpha
i}\left[\mathcal{L}_{11}^{B,ij\gamma}\mathcal{L}_{12}^{jk}-\mathcal{L}_{12}^{B,ij\gamma}\mathcal{L}_{11}^{jk}\right]\left(\mathcal{L}_{11}^{-1}\right)^{k\beta}$
(18) $\displaystyle\mu_{H,\alpha\beta\gamma}$
$\displaystyle=\left(\mathcal{L}_{11}^{-1}\right)^{\alpha
i}\mathcal{L}_{11}^{B,i\beta\gamma}$ (19)
$\displaystyle\mu_{T,\alpha\beta\gamma}$
$\displaystyle=\left(\mathcal{L}_{12}^{-1}\right)^{\alpha
i}\mathcal{L}_{12}^{B,i\beta\gamma}$ (20)

##  2 Transport equations and optical elements

###  2.1 Kubo kernel expressions

The equations implemented in LinReTraCe result from the analytic integration
of Eqs. (9-10), assuming a Lorentzian spectral function $A(\omega$)

$A_{\mathbf{k}n}(\omega)=\frac{Z_{\mathbf{k}n}}{\pi}\frac{\Gamma_{\mathbf{k}n}}{(\omega-
a_{\mathbf{k}n})^{2}+\Gamma_{\mathbf{k}n}^{2}}.$ (21)

This form is motivated by a linearization of the self-energy (for a
discussion, see Ref. [2])444 Eq. (22) shows the expansion (to linear order) of
the self-energy around the Fermi level $(\omega=0$). Typically, an expansion
around the finite quasi-particle energy $a_{\mathbf{k}n}$ is more accurate and
can be provided by the user.

$\Sigma_{\mathbf{k}n}(\omega)\approx\Re\Sigma_{\mathbf{k}n}(0)+(1-Z_{\mathbf{k}n}^{-1})\omega-i\Gamma^{0}_{\mathbf{k}n}.$
(22)

Eq. (21) describes a quasi-particle peak of weight $Z_{\mathbf{k}n}$ with
renormalized scattering rates
$\Gamma_{\mathbf{k}n}=Z^{\vphantom{0}}_{\mathbf{k}n}\Gamma^{0}_{\mathbf{k}n}$
and renormalized energies
$a_{\mathbf{k}n}=Z_{\mathbf{k}n}(\varepsilon_{\mathbf{k}n}^{0}-\mu+\Re\Sigma_{\mathbf{k}n}(0))$.
Then, as detailed in Ref. [2], the kernels Eqs. (9,-10) can be calculated
analytically. For the intra-band kernels ($n\equiv m$) one finds, with
polygamma functions $\psi_{n}(z)$ [67] evaluated at
$z=\frac{1}{2}+\frac{\beta}{2\pi}(\Gamma+ia)$:

$\displaystyle\mathcal{K}_{11}(\mathbf{k},n)=\frac{Z^{2}\beta}{4\pi^{3}\Gamma}\bigg{[}\Re\psi_{1}(z)-\frac{\beta\Gamma}{2\pi}\Re\psi_{2}(z)\bigg{]}$
(23)
$\displaystyle\mathcal{K}_{12}(\mathbf{k},n)=\frac{Z^{2}\beta}{4\pi^{3}\Gamma}\bigg{[}a\Re\psi_{1}(z)-\frac{a\beta\Gamma}{2\pi}\Re\psi_{2}(z)-\frac{\beta\Gamma^{2}}{2\pi}\Im\psi_{2}(z)\bigg{]}$
(24)
$\displaystyle\mathcal{K}_{22}(\mathbf{k},n)=\frac{Z^{2}\beta}{4\pi^{3}\Gamma}\bigg{[}(a^{2}+\Gamma^{2})\Re\psi_{1}(z)+\frac{\beta\Gamma}{2\pi}\left(\Gamma^{2}-a^{2}\right)\Re\psi_{2}(z)-\frac{a\beta\Gamma^{2}}{\pi}\Im\psi_{2}(z)\bigg{]}$
(25)
$\displaystyle\mathcal{K}^{B}_{11}(\mathbf{k},n)=\frac{Z^{3}\beta}{16\pi^{4}\Gamma^{2}}\bigg{[}3\Re\psi_{1}(z)-\frac{3\beta\Gamma}{2\pi}\Re\psi_{2}(z)+\frac{\beta^{2}\Gamma^{2}}{4\pi^{2}}\Re\psi_{3}(z)\bigg{]}$
(26)
$\displaystyle\mathcal{K}^{B}_{12}(\mathbf{k},n)=\frac{Z^{3}\beta}{16\pi^{4}\Gamma^{2}}\bigg{[}3a\Re\psi_{1}(z)-\frac{3a\beta\Gamma}{2\pi}\Re\psi_{2}(z)-\frac{\beta\Gamma^{2}}{2\pi}\Im\psi_{2}(z)+\frac{a\beta^{2}\Gamma^{2}}{4\pi^{2}}\Re\psi_{3}(z)+\frac{\beta^{2}\Gamma^{3}}{4\pi^{2}}\Im\psi_{3}(z)\bigg{]}$
(27) $\displaystyle\begin{aligned}
\mathcal{K}^{B}_{22}(\mathbf{k},n)=\frac{Z^{3}\beta}{16\pi^{4}\Gamma^{2}}\bigg{[}&(3a^{2}+\Gamma^{2})\Re\psi_{1}(z)-\frac{\beta\Gamma}{2\pi}(3a^{2}+\Gamma^{2})\Re\psi_{2}(z)-\frac{a\beta\Gamma^{2}}{\pi}\Im\psi_{2}(z)\\\
&-\frac{\beta^{2}\Gamma^{2}}{4\pi^{2}}(\Gamma^{2}-a^{2})\Re\psi_{3}(z)+\frac{a\beta^{2}\Gamma^{3}}{2\pi^{2}}\Im\psi_{3}(z)\bigg{]}\end{aligned}$
(28)

The inter-band ($n\neq m$) kernels become

$\displaystyle\begin{aligned}
\mathcal{K}_{11}(\mathbf{k},n,m)=&\frac{Z_{1}Z_{2}\beta}{2\pi^{3}\left[(a_{1}-a_{2})^{2}+(\Gamma_{1}-\Gamma_{2})^{2}\right]\left[(a_{1}-a_{2})^{2}+(\Gamma_{1}+\Gamma_{2})^{2}\right]}\\\
\times\Bigg{[}&\Re\Big{\\{}\left[(a_{1}-a_{2})^{2}+\Gamma_{2}^{2}-\Gamma_{1}^{2}-2i\Gamma_{1}(a_{2}-a_{1})\right]\Gamma_{2}\psi_{1}\left(z_{1}\right)\Big{\\}}\\\
+&\Re\Big{\\{}\left[(a_{2}-a_{1})^{2}+\Gamma_{1}^{2}-\Gamma_{2}^{2}-2i\Gamma_{2}(a_{1}-a_{2})\right]\Gamma_{1}\psi_{1}\left(z_{2}\right)\Big{\\}}\Bigg{]}\end{aligned}$
(29) $\displaystyle\begin{aligned}
\mathcal{K}_{12}(\mathbf{k},n,m)=&\frac{Z_{1}Z_{2}\beta}{2\pi^{3}\left[(a_{1}-a_{2})^{2}+(\Gamma_{1}-\Gamma_{2})^{2}\right]\left[(a_{1}-a_{2})^{2}+(\Gamma_{1}+\Gamma_{2})^{2}\right]}\\\
\times\Bigg{[}&\Re\Big{\\{}\left(a_{1}-i\Gamma_{1}\right)\left[\left(a_{1}-a_{2}\right)^{2}+\Gamma_{2}^{2}-\Gamma_{1}^{2}-2i\left(a_{2}-a_{1}\right)\Gamma_{1}\right]\Gamma_{2}\psi_{1}(z_{1})\Big{\\}}\\\
+&\Re\Big{\\{}\left(a_{2}-i\Gamma_{2}\right)\left[\left(a_{2}-a_{1}\right)^{2}+\Gamma_{1}^{2}-\Gamma_{2}^{2}-2i\left(a_{1}-a_{2}\right)\Gamma_{2}\right]\Gamma_{1}\psi_{1}(z_{2})\Big{\\}}\Bigg{]}\end{aligned}$
(30) $\displaystyle\begin{aligned}
\mathcal{K}_{22}(\mathbf{k},n,m)=&\frac{Z_{1}Z_{2}\beta}{2\pi^{3}\left[(a_{1}-a_{2})^{2}+(\Gamma_{1}-\Gamma_{2})^{2}\right]\left[(a_{1}-a_{2})^{2}+(\Gamma_{1}+\Gamma_{2})^{2}\right]}\\\
\times\Bigg{[}&\Re\Big{\\{}\left(a_{1}-i\Gamma_{1}\right)^{2}\left[\left(a_{1}-a_{2}\right)^{2}+\Gamma_{2}^{2}-\Gamma_{1}^{2}-2i\left(a_{2}-a_{1}\right)\Gamma_{1}\right]\Gamma_{2}\psi_{1}\left(z_{1}\right)\Big{\\}}\\\
+&\Re\Big{\\{}\left(a_{2}-i\Gamma_{2}\right)^{2}\left[\left(a_{2}-a_{1}\right)^{2}+\Gamma_{1}^{2}-\Gamma_{2}^{2}-2i\left(a_{1}-a_{2}\right)\Gamma_{2}\right]\Gamma_{1}\psi_{1}\left(z_{2}\right)\Big{\\}}\Bigg{]},\end{aligned}$
(31)

where $\psi_{1}(z_{1/2})$ is evaluated at
$z_{1/2}=\frac{1}{2}+\frac{\beta}{2\pi}(\Gamma_{1/2}+ia_{1/2})$. These kernels
represent a generalization of Eqs. (23-25) and per Eq. (9) obey band-swapping
symmetry
$\mathcal{K}_{ab}(\mathbf{k},n,m)\equiv\mathcal{K}_{ab}(\mathbf{k},m,n)$. See
Appendix D for the generic intra-band limit ($f_{m}\rightarrow f_{n}$,
$f=Z,\Gamma,a$).

The magnetic inter-band ($n\neq m$) kernels evaluate to

$\displaystyle\mathcal{K}_{11}^{B}(\mathbf{k},n,m)=$
$\displaystyle\frac{Z_{1}Z_{2}^{2}\Gamma_{1}\Gamma_{2}^{2}\beta}{2\pi^{4}}$
(32) $\displaystyle\times$
$\displaystyle\Bigg{[}\Re\Big{\\{}\frac{1}{\Gamma_{1}}\psi_{1}(z_{1})\frac{1}{\left[a_{1}-a_{2}-i(\Gamma_{1}+\Gamma_{2})\right]^{2}\left[a_{1}-a_{2}-i(\Gamma_{1}-\Gamma_{2})\right]^{2}}\Big{\\}}$
$\displaystyle-$
$\displaystyle\Re\Big{\\{}\frac{\beta}{4\pi\Gamma_{2}^{2}}\psi_{2}(z_{2})\frac{1}{\left[a_{2}-a_{1}-i(\Gamma_{1}+\Gamma_{2})\right]\left[a_{2}-a_{1}+i(\Gamma_{1}-\Gamma_{2})\right]}\Big{\\}}$
$\displaystyle+$
$\displaystyle\Im\Big{\\{}\frac{1}{\Gamma_{2}^{2}}\psi_{1}(z_{2})\frac{(a_{2}-a_{1}-i\Gamma_{2})}{\left[a_{2}-a_{1}-i(\Gamma_{1}+\Gamma_{2})\right]^{2}\left[a_{2}-a_{1}+i(\Gamma_{1}-\Gamma_{2})\right]^{2}}\Big{\\}}$
$\displaystyle+$
$\displaystyle\Re\Big{\\{}\frac{1}{2\Gamma_{2}^{3}}\psi_{1}(z_{2})\frac{1}{\left[a_{2}-a_{1}-i(\Gamma_{1}+\Gamma_{2})\right]\left[a_{2}-a_{1}+i(\Gamma_{1}-\Gamma_{2})\right]}\Big{\\}}\Bigg{]}$

$\displaystyle\mathcal{K}_{12}^{B}(\mathbf{k},n,m)=$
$\displaystyle\frac{Z_{1}Z_{2}^{2}\Gamma_{1}\Gamma_{2}^{2}\beta}{2\pi^{4}}$
(33) $\displaystyle\times$
$\displaystyle\Bigg{[}\Re\Big{\\{}\frac{1}{\Gamma_{1}}\psi_{1}(z_{1})\frac{(a_{1}-i\Gamma_{1})}{\left[a_{1}-a_{2}-i(\Gamma_{1}+\Gamma_{2})\right]^{2}\left[a_{1}-a_{2}-i(\Gamma_{1}-\Gamma_{2})\right]^{2}}\Big{\\}}$
$\displaystyle-$
$\displaystyle\Im\Big{\\{}\frac{1}{2\Gamma_{2}^{2}}\psi_{1}(z_{2})\frac{1}{\left[a_{2}-a_{1}-i(\Gamma_{1}+\Gamma_{2})\right]\left[a_{2}-a_{1}+i(\Gamma_{1}-\Gamma_{2})\right]}\Big{\\}}$
$\displaystyle-$
$\displaystyle\Re\Big{\\{}\frac{\beta}{4\pi\Gamma_{2}^{2}}\psi_{2}(z_{2})\frac{(a_{2}-i\Gamma_{2})}{\left[a_{2}-a_{1}-i(\Gamma_{1}+\Gamma_{2}]\right)\left(a_{2}-a_{1}+i(\Gamma_{1}-\Gamma_{2})\right)}\Big{\\}}$
$\displaystyle+$
$\displaystyle\Im\Big{\\{}\frac{1}{\Gamma_{2}^{2}}\psi_{1}(z_{2})\frac{(a_{2}-i\Gamma_{2})(a_{2}-a_{1}-i\Gamma_{2})}{\left[a_{2}-a_{1}-i(\Gamma_{1}+\Gamma_{2})\right]^{2}\left[a_{2}-a_{1}+i(\Gamma_{1}-\Gamma_{2})\right]^{2}}\Big{\\}}$
$\displaystyle+$
$\displaystyle\Re\Big{\\{}\frac{1}{2\Gamma_{2}^{3}}\psi_{1}(z_{2})\frac{(a_{2}-i\Gamma_{2})}{\left[a_{2}-a_{1}-i(\Gamma_{1}+\Gamma_{2})\right]\left[a_{2}-a_{1}+i(\Gamma_{1}-\Gamma_{2})\right]}\Big{\\}}\Bigg{]}$

$\displaystyle\mathcal{K}_{22}^{B}(\mathbf{k},n,m)=$
$\displaystyle\frac{Z_{1}Z_{2}^{2}\Gamma_{1}\Gamma_{2}^{2}\beta}{2\pi^{4}}$
(34) $\displaystyle\times$
$\displaystyle\Bigg{[}\Re\Big{\\{}\frac{1}{\Gamma_{1}}\psi_{1}(z_{1})\frac{(a_{1}-i\Gamma_{1})^{2}}{\left[a_{1}-a_{2}-i(\Gamma_{1}+\Gamma_{2})\right]^{2}\left[a_{1}-a_{2}-i(\Gamma_{1}-\Gamma_{2})\right]^{2}}\Big{\\}}$
$\displaystyle-$
$\displaystyle\Im\Big{\\{}\frac{1}{\Gamma_{2}^{2}}\psi_{1}(z_{2})\frac{(a_{2}-i\Gamma_{2})}{\left[a_{2}-a_{1}-i(\Gamma_{1}+\Gamma_{2})\right]\left[a_{2}-a_{1}+i(\Gamma_{1}-\Gamma_{2})\right]}\Big{\\}}$
$\displaystyle-$
$\displaystyle\Re\Big{\\{}\frac{\beta}{4\pi\Gamma_{2}^{2}}\psi_{2}(z_{2})\frac{(a_{2}-i\Gamma_{2})^{2}}{\left[a_{2}-a_{1}-i(\Gamma_{1}+\Gamma_{2})\right]\left(a_{2}-a_{1}+i(\Gamma_{1}-\Gamma_{2})\right)}\Big{\\}}$
$\displaystyle+$
$\displaystyle\Im\Big{\\{}\frac{1}{\Gamma_{2}^{2}}\psi_{1}(z_{2})\frac{(a_{2}-i\Gamma_{2})^{2}(a_{2}-a_{1}-i\Gamma_{2})}{\left[a_{2}-a_{1}-i(\Gamma_{1}+\Gamma_{2})\right]^{2}\left[a_{2}-a_{1}+i(\Gamma_{1}-\Gamma_{2})\right]^{2}}\Big{\\}}$
$\displaystyle+$
$\displaystyle\Re\Big{\\{}\frac{1}{2\Gamma_{2}^{3}}\psi_{1}(z_{2})\frac{(a_{2}-i\Gamma_{2})^{2}}{\left[a_{2}-a_{1}-i(\Gamma_{1}+\Gamma_{2})\right]\left[a_{2}-a_{1}+i(\Gamma_{1}-\Gamma_{2})\right]}\Big{\\}}\Bigg{]}$

where $\psi_{n}(z_{1/2})$ is, again, evaluated at
$z_{1/2}=\frac{1}{2}+\frac{\beta}{2\pi}(\Gamma_{1/2}+ia_{1/2})$.

###  2.2 Boltzmann approximation

The semiclassical Boltzmann transport expressions can be obtained by expanding
the polygamma functions of our approach around $\Gamma=0$

$\psi_{n}\left(\frac{1}{2}+\frac{\beta}{2\pi}(\Gamma+ia)\right)=\psi_{n}\left(\frac{1}{2}+\frac{i\beta
a}{2\pi}\right)+\frac{\beta\Gamma}{2\pi}\psi_{n+1}\left(\frac{1}{2}+\frac{i\beta
a}{2\pi}\right)+\mathcal{O}\left(\Gamma^{2}\right),$ (35)

and recognizing that

$\displaystyle f_{\mathrm{FD}}(a)$
$\displaystyle=\frac{1}{2}-\frac{1}{\pi}\Im\psi\left(\frac{1}{2}+\frac{i\beta
a}{2\pi}\right)$ (36) $\displaystyle-\partial_{a}f_{\mathrm{FD}}(a)$
$\displaystyle=\frac{\beta}{2\pi^{2}}\Re\psi_{1}\left(\frac{1}{2}+\frac{i\beta
a}{2\pi}\right).$ (37)

The Boltzmann intra-band expressions are then recovered from the expansions as
the leading terms in $\Gamma$. Our quantum mechanical formalism thus contains
the semi-classical description as the coherent (infinite lifetime) limit. To
ensure that the Boltzmann inter-band kernels reduce to the Boltzmann intra-
band expressions in the limit of degenerate states (with the same lifetime),
they have to include terms beyond the leading order. Instead, if the Boltzmann
approximation only takes into account the leading terms, the limit
$a_{1}\rightarrow a_{2};\Gamma_{1}\rightarrow\Gamma_{2}$ will yield
inconsistent results. To our knowledge this ensemble of inter-band Boltzmann
expressions has not been derived previously:

$\displaystyle\mathcal{K}^{\mathrm{Boltzmann}}_{11}(\mathbf{k},n)=-\frac{Z^{2}}{2\pi\Gamma}\partial_{a}f_{\mathrm{FD}}(a)$
(38)
$\displaystyle\mathcal{K}^{\mathrm{Boltzmann}}_{12}(\mathbf{k},n)=-\frac{aZ^{2}}{2\pi\Gamma}\partial_{a}f_{\mathrm{FD}}(a)$
(39)
$\displaystyle\mathcal{K}^{\mathrm{Boltzmann}}_{22}(\mathbf{k},n)=-\frac{\left(a^{2}+\Gamma^{2}\right)Z^{2}}{2\pi\Gamma}\partial_{a}f_{\mathrm{FD}}(a)$
(40)
$\displaystyle\mathcal{K}^{\mathrm{Boltzmann},B}_{11}(\mathbf{k},n)=-\frac{3Z^{3}}{8\pi^{2}\Gamma^{2}}\partial_{a}f_{\mathrm{FD}}(a)$
(41)
$\displaystyle\mathcal{K}^{\mathrm{Boltzmann},B}_{12}(\mathbf{k},n)=-\frac{3aZ^{3}}{8\pi^{2}\Gamma^{2}}\partial_{a}f_{\mathrm{FD}}(a)$
(42)
$\displaystyle\mathcal{K}^{\mathrm{Boltzmann},B}_{22}(\mathbf{k},n)=-\frac{\left(3a^{2}+\Gamma^{2}\right)Z^{3}}{8\pi^{2}\Gamma^{2}}\partial_{a}f_{\mathrm{FD}}(a)$
(43)

$\displaystyle\begin{aligned}
&\mathcal{K}^{\mathrm{Boltzmann}}_{11}(\mathbf{k},n,m)=-\frac{Z_{1}Z_{2}}{\pi\left[(a_{1}-a_{2})^{2}+(\Gamma_{1}-\Gamma_{2})^{2}\right]\left[(a_{1}-a_{2})^{2}+(\Gamma_{1}+\Gamma_{2})^{2}\right]}\\\
&\times\Bigg{[}\left[(a_{1}-a_{2})^{2}+\Gamma_{2}^{2}-\Gamma_{1}^{2}\right]\Gamma_{2}\partial_{a_{1}}f_{\mathrm{FD}}(a_{1})+\left[(a_{2}-a_{1})^{2}+\Gamma_{1}^{2}-\Gamma_{2}^{2}\right]\Gamma_{1}\partial_{a_{2}}f_{\mathrm{FD}}(a_{2})\Bigg{]}\end{aligned}$
(44) $\displaystyle\begin{aligned}
&\mathcal{K}^{\mathrm{Boltzmann}}_{12}(\mathbf{k},n,m)=-\frac{Z_{1}Z_{2}}{\pi\left[(a_{1}-a_{2})^{2}+(\Gamma_{1}-\Gamma_{2})^{2}\right]\left[(a_{1}-a_{2})^{2}+(\Gamma_{1}+\Gamma_{2})^{2}\right]}\\\
&\times\Bigg{[}\left\\{\left[(a_{1}-a_{2})^{2}+\Gamma_{2}^{2}-\Gamma_{1}^{2}\right]a_{1}-2\Gamma_{1}^{2}(a_{1}-a_{2})\right\\}\Gamma_{2}\partial_{a_{1}}f_{\mathrm{FD}}(a_{1})\\\
&\;\;\;+\left\\{\left[(a_{2}-a_{1})^{2}+\Gamma_{1}^{2}-\Gamma_{2}^{2}\right]a_{2}-2\Gamma_{2}^{2}(a_{2}-a_{1})\right\\}\Gamma_{1}\partial_{a_{2}}f_{\mathrm{FD}}(a_{2})\Bigg{]}\end{aligned}$
(45) $\displaystyle\begin{aligned}
&\mathcal{K}^{\mathrm{Boltzmann}}_{22}(\mathbf{k},n,m)=-\frac{Z_{1}Z_{2}}{\pi\left[(a_{1}-a_{2})^{2}+(\Gamma_{1}-\Gamma_{2})^{2}\right]\left[(a_{1}-a_{2})^{2}+(\Gamma_{1}+\Gamma_{2})^{2}\right]}\\\
&\times\Bigg{[}\left\\{\left[(a_{1}-a_{2})^{2}+\Gamma_{2}^{2}-\Gamma_{1}^{2}\right](a_{1}^{2}-\Gamma_{1}^{2})-4\Gamma_{1}^{2}a_{1}(a_{1}-a_{2})\right\\}\Gamma_{2}\partial_{a_{1}}f_{\mathrm{FD}}(a_{1})\\\
&\;\;\;+\left\\{\left[(a_{2}-a_{1})^{2}+\Gamma_{1}^{2}-\Gamma_{2}^{2}\right](a_{2}^{2}-\Gamma_{2}^{2})-4\Gamma_{2}^{2}a_{2}(a_{2}-a_{1})\right\\}\Gamma_{1}\partial_{a_{2}}f_{\mathrm{FD}}(a_{2})\Bigg{]}\end{aligned}$
(46)

$\displaystyle\mathcal{K}_{11}^{\mathrm{Boltzmann},B}$
$\displaystyle(\mathbf{k},n,m)=-\frac{Z_{1}Z_{2}^{2}}{2\pi^{2}\Gamma_{2}}$
(47) $\displaystyle\times$
$\displaystyle\Bigg{[}\Re\Big{\\{}\frac{2\Gamma_{2}^{3}}{\left[a_{1}-a_{2}-i(\Gamma_{1}+\Gamma_{2})\right]^{2}\left[a_{1}-a_{2}-i(\Gamma_{1}-\Gamma_{2})\right]^{2}}\Big{\\}}\partial_{a_{1}}f_{\mathrm{FD}}(a_{1})$
$\displaystyle+$
$\displaystyle\Im\Big{\\{}\frac{2\Gamma_{1}\Gamma_{2}(a_{2}-a_{1}-i\Gamma_{2})}{\left[a_{2}-a_{1}-i(\Gamma_{1}+\Gamma_{2})\right]^{2}\left[a_{2}-a_{1}+i(\Gamma_{1}-\Gamma_{2})\right]^{2}}\Big{\\}}\partial_{a_{2}}f_{\mathrm{FD}}(a_{2})$
$\displaystyle+$
$\displaystyle\Re\Big{\\{}\frac{\Gamma_{1}}{\left[a_{2}-a_{1}-i(\Gamma_{1}+\Gamma_{2})\right]\left[a_{2}-a_{1}+i(\Gamma_{1}-\Gamma_{2})\right]}\Big{\\}}\partial_{a_{2}}f_{\mathrm{FD}}(a_{2})\Bigg{]}$

$\displaystyle\mathcal{K}_{12}^{\mathrm{Boltzmann},B}$
$\displaystyle(\mathbf{k},n,m)=-\frac{Z_{1}Z_{2}^{2}\beta}{2\pi^{2}\Gamma_{2}}$
(48) $\displaystyle\times$
$\displaystyle\Bigg{[}\Re\Big{\\{}\frac{2\Gamma_{2}^{3}(a_{1}-i\Gamma_{1})}{\left[a_{1}-a_{2}-i(\Gamma_{1}+\Gamma_{2})\right]^{2}\left[a_{1}-a_{2}-i(\Gamma_{1}-\Gamma_{2})\right]^{2}}\Big{\\}}\partial_{a_{1}}f_{\mathrm{FD}}(a_{1})$
$\displaystyle-$
$\displaystyle\Im\Big{\\{}\frac{\Gamma_{1}\Gamma_{2}}{\left[a_{2}-a_{1}-i(\Gamma_{1}+\Gamma_{2})\right]\left[a_{2}-a_{1}+i(\Gamma_{1}-\Gamma_{2})\right]}\Big{\\}}\partial_{a_{2}}f_{\mathrm{FD}}(a_{2})$
$\displaystyle+$
$\displaystyle\Im\Big{\\{}\frac{2\Gamma_{1}\Gamma_{2}(a_{2}-i\Gamma_{2})(a_{2}-a_{1}-i\Gamma_{2})}{\left[a_{2}-a_{1}-i(\Gamma_{1}+\Gamma_{2})\right]^{2}\left[a_{2}-a_{1}+i(\Gamma_{1}-\Gamma_{2})\right]^{2}}\Big{\\}}\partial_{a_{2}}f_{\mathrm{FD}}(a_{2})$
$\displaystyle+$
$\displaystyle\Re\Big{\\{}\frac{\Gamma_{1}(a_{2}-i\Gamma_{2})}{\left[a_{2}-a_{1}-i(\Gamma_{1}+\Gamma_{2})\right]\left[a_{2}-a_{1}+i(\Gamma_{1}-\Gamma_{2})\right]}\Big{\\}}\partial_{a_{2}}f_{\mathrm{FD}}(a_{2})\Bigg{]}$

$\displaystyle\mathcal{K}_{22}^{\mathrm{Boltzmann},B}$
$\displaystyle(\mathbf{k},n,m)=-\frac{Z_{1}Z_{2}^{2}}{2\pi^{2}\Gamma_{2}}$
(49) $\displaystyle\times$
$\displaystyle\Bigg{[}\Re\Big{\\{}\frac{2\Gamma_{2}^{3}(a_{1}-i\Gamma_{1})^{2}}{\left[a_{1}-a_{2}-i(\Gamma_{1}+\Gamma_{2})\right]^{2}\left[a_{1}-a_{2}-i(\Gamma_{1}-\Gamma_{2})\right]^{2}}\Big{\\}}\partial_{a_{1}}f_{\mathrm{FD}}(a_{1})$
$\displaystyle-$
$\displaystyle\Im\Big{\\{}\frac{2\Gamma_{1}\Gamma_{2}(a_{2}-i\Gamma_{2})}{\left[a_{2}-a_{1}-i(\Gamma_{1}+\Gamma_{2})\right]\left[a_{2}-a_{1}+i(\Gamma_{1}-\Gamma_{2})\right]}\Big{\\}}\partial_{a_{2}}f_{\mathrm{FD}}(a_{2})$
$\displaystyle+$
$\displaystyle\Im\Big{\\{}\frac{2\Gamma_{1}\Gamma_{2}(a_{2}-i\Gamma_{2})^{2}(a_{2}-a_{1}-i\Gamma_{2})}{\left[a_{2}-a_{1}-i(\Gamma_{1}+\Gamma_{2})\right]^{2}\left[a_{2}-a_{1}+i(\Gamma_{1}-\Gamma_{2})\right]^{2}}\Big{\\}}\partial_{a_{2}}f_{\mathrm{FD}}(a_{2})$
$\displaystyle+$
$\displaystyle\Re\Big{\\{}\frac{\Gamma_{1}(a_{2}-i\Gamma_{2})^{2}}{\left[a_{2}-a_{1}-i(\Gamma_{1}+\Gamma_{2})\right]\left[a_{2}-a_{1}+i(\Gamma_{1}-\Gamma_{2})\right]}\Big{\\}}\partial_{a_{2}}f_{\mathrm{FD}}(a_{2})\Bigg{]}$

In all, these Boltzmann approximations stem from a Taylor expansion of the
terms in square brackets in Eqs. (23-28). Notably, there, terms linear in
$\Gamma$ cancel exactly in the expansion of the _intra_ band kernels, thus
higher order terms ‘only‘ become important if $\beta\Gamma$ is significant in
size. Consequently Boltzmann results are accurate (comparable to Kubo) at
elevated temperatures despite failing to fulfill fundamental theorems of
thermodynamics in the zero temperature limit.

###  2.3 Quasi-particle renormalizations

The kernel expressions in Sec. 2.1, based on the self-energy linearization Eq.
(22), clearly exhibit a non-trivial dependency on the quasi-particle
renormalization: $Z$-factors emerge, both, as part of the overall pre-factor
(à la Boltzmann) as well as in the argument of the polygamma functions via a
renormalization of the energy ($a$) and the scattering ($\Gamma$):

Consider a _metallic_ system with bare dispersion
$\varepsilon_{0}(\mathbf{k})$ and a bare scattering rate $\Gamma_{0}$ at
temperatures where the Boltzmann approximation is accurate. If the same quasi-
particle renormalization is applied to all states, a partial cancellation of
the pre-factors in Eqs. (38-43) occurs due to renormalization of the
scattering rate $\Gamma=Z\Gamma_{0}$. The Boltzmann response kernels then
verify

$\mathcal{K}^{(B)}_{ij}\propto Za^{i+j-2}\partial_{a}f_{\mathrm{FD}}(a).$ (50)

The conductivity kernel ($i=j=1$) that is even in the energy $a=Z\epsilon_{0}$
is to good approximation unaffected by $Z$ due to a compensation of two
effects: $Z$ simultaneously decreases the weights of the selected states, and
pushes more states into the (thermal) selection window through band-
narrowing.555This is exact for $0<Z\leq 1$, in the limit of infinite bandwidth
and a flat density of states. In realistic scenarios, however, there can be a
notable $Z$-dependence for strong renormalizations and narrow band-widths,
elevated temperatures, or a strongly energy-dependent density of states. Odd
kernels ($\mathcal{L}_{12}$, $\mathcal{L}^{B}_{12}$) on the other hand
distinguish between electron and hole contributions via the sign of $a$.
Through this differentiation the overall summation will be tilted to either
direction depending on the asymmetry of the system. Energy renormalization in
this context then can be thought of as an amplification of this asymmetry,
increasing the non-interacting signal666By assuming a linearized density of
states centered around the thermal selection window,
$D(\varepsilon)=D_{0}+\alpha\varepsilon$, it is apparent that only the linear
(constant) term is responsible for finite values of kernels that are odd
(even) in the energy $a$. Quasi-particle renormalizations drop out for the
constant term and amplify the effect of the linear term, leading to the
increase of $\frac{1}{Z}$ of Eq. (51).

$\mathcal{K}_{12}^{(B)}(Z)\propto\frac{1}{Z}\mathcal{K}_{12}^{(B)}(Z=1)$ (51)

providing a correlation mechanism to boost the Seebeck[68, 69] and Nernst
coefficient, realized, e.g., in heavy-fermion systems[35]. The above arguments
in general do not hold for insulating systems where we find a more nuanced
interplay of band gap, energies, chemical potential and quasi-particle
renormalization. At the very least $Z<1$ will result in a band gap reduction
$\Delta=Z\Delta_{0}$, which affects thermal activation, and hence conduction,
exponentially.

###  2.4 Matrix elements

####  2.4.1 Dipole optical elements

Given the Fermi velocities (matrix elements of the momentum operator)

$v_{\mathbf{k}nn^{\prime}}^{\alpha}=\frac{1}{m}\left\langle\mathbf{k}n^{\prime}|\mathcal{P}_{\alpha}|\mathbf{k}n\right\rangle$
(52)

with $\alpha$ indicating a Cartesian direction, $m$ the electron mass, and
$\langle\hbox{{r}}|\hbox{{k}}n\rangle=\chi_{\mathbf{k}n}(\hbox{{r}})$ a band-
momentum basis, the amplitude of optical dipole ($\hbox{{q}}=0$) transitions
is given by

$M^{\alpha\beta}(\hbox{{k}},n,m)=v_{\mathbf{k}nm}^{\alpha*}v_{\mathbf{k}mn}^{\beta}\\\
$ (53)

and can be calculated within band-theory; for WIEN2k’s implementation see Ref.
[70].

####  2.4.2 Peierls approximation

In tight-binding or model settings, in which there is no access to
wavefunctions, the above matrix elements cannot be calculated. Instead, one
couples the electromagnetic vector potential directly to the lattice fermions
using the Peierls substitution approach[71]. Following this (approximate)
procedure, Fermi velocities are momentum-derivatives of the (one-particle)
Hamiltonian. Performing the derivative in the band-basis, there are only
intra-band velocities,
$v_{\mathbf{k}nm}^{\alpha}\propto\delta_{nm}v_{\mathbf{k}n}^{\alpha}$, for
which

$v^{\alpha}_{\mathbf{k}n}=\frac{1}{\hbar}\partial_{\mathbf{k}_{\alpha}}\varepsilon_{n}(\mathbf{k}),$
(54)

and

$\displaystyle M^{\alpha\beta}(\hbox{{k}},n,n)$
$\displaystyle=v_{\mathbf{k}n}^{\alpha}v_{\mathbf{k}n}^{\beta}.$ (55)

Using the band-curvatures
$c^{\alpha\beta}_{k}=1/\hbar\,\partial_{k_{\alpha}}\partial_{k_{\beta}}\varepsilon(\hbox{{k}})$,
also the matrix elements for magnetic quantities can be derived. One finds[72]

$\displaystyle M^{B,\alpha\beta\gamma}(\hbox{{k}},n,n)$
$\displaystyle=\varepsilon_{\gamma ij}v^{\alpha}_{\mathbf{k}n}c^{\beta
i}_{\mathbf{k}n}v^{j}_{\mathbf{k}n}$ (56)

with $\varepsilon_{ijk}$ the Levy-Civita symbol in three dimensions. These
expressions for $M^{(B)}$ are common to LinReTraCe, BoltzTraP and other codes.

Clearly, however, taking momentum-derivatives does not commute with a general
basis transformation $U$:
$\partial_{k}(U^{\dagger}(\hbox{{k}})H(\hbox{{k}})U(\hbox{{k}}))\neq
U^{\dagger}(\hbox{{k}})(\partial_{k}H(\hbox{{k}}))U(\hbox{{k}})$. Therefore,
transport properties using the Peierls approach will be basis-dependent[73].
One can show[43] that the Peierls approximation is best (closest to the true
dipole element) the more localized the basis is. In the tight-binding and
Wannier mode, LinReTraCe therefore performs the momentum-derivative in the
local/Wannier basis
$\chi_{\mathbf{R}l}(\hbox{{r}})=\langle\hbox{{r}}|\hbox{{R}}l\rangle$:

$v_{\mathbf{k}ll^{\prime}}^{\alpha}=\frac{1}{\hbar}\partial_{k_{\alpha}}H^{ll^{\prime}}(\hbox{{k}})-i(\rho_{l^{\prime}}^{\alpha}-\rho_{l}^{\alpha})H^{ll^{\prime}}(\hbox{{k}})$
(57)

where $H(\hbox{{k}})$ is the Fourier transform of $H(\hbox{{R}})$. In this
generalized Peierls approach[43], the second term arises for unit-cells with
more than one atom, with intra-cell coordinates $\rho_{l}$ of the atom hosting
orbital $l$. This extra-term777Sometimes it is absorbed into the Hamiltonian
${H}(\hbox{{k}})$[74]. in particular ensures that an arbitrary extension of
the unit-cell (conventional cell or equivalent supercells) gives the same
result as calculations for the primitive unit-cell. Velocities evaluated in
the local basis (orbitals indexed with $l$) are then rotated into the band-
basis (band-index $n$). Because of the mentioned non-commutation of momentum-
derivative and basis-transformation, the generalized Peierls approach may
yield inter-band transitions à la Eq. (53) that are absent in Eq. (55). Eq. 57
was derived using the minimal coupling prescription for the vector potential
in the continuum instead of the Peierls substitution for the lattice[43]. A
derivation of the magnetic response beyond the Peierls approach has not yet
been performed. Magnetic transport functions in LinReTraCe therefore rely, via
Eq. (56), on group velocities and band-curvatures.

##  3 Implementation

###  3.1 Requirements

_Hardware_
— processor and compiler combination supporting quadruple precision (e.g., x86
or x86-64 architectures)

_Programming languages_
— Fortran 95 and Python 3

_Required dependencies_
— HDF5 ($\geq$ 1.12.1)
— h5py, numpy, scipy, ase ($\geq 3.18.0$), spglib ($\geq$ 1.9.5)

_Optional dependencies_
— MPI (Fortran 95), matplotlib, boltztrap2 ($\geq$ 20.7.1), cmake, pip, git

###  3.2 Installation

The program package can be obtained directly from the command line

$ git clone https://github.com/linretrace/linretrace

or by downloading the files from the github page with a browser. The pre- and
post-processing is handled via Python3 where HDF5 data files are interfaced
with `h5py` and computationally costly work is performed with `scipy` and
`numpy`. `Ase` [75] and `Spglib` [76] are used to detect crystal symmetries
and create irreducible momentum meshes, if they are not provided by the DFT
code. In order to achieve maximal operability, `BoltzTraP2` [24], as well as
`matplotlib` are recommended. `BoltzTraP2` is interfaced to interpolate bands,
necessary to generate band derivatives and curvatures used in the magnetic
optical elements, and `matplotlib` is used for graphical plotting. A quick way
to obtain all (necessary and optional) dependencies is via the pip package
manager

$ pip install matplotlib h5py numpy scipy ase spglib $ pip install boltztrap2

After that, the LinReTraCe Python scripts can be executed directly from the
source folder. Alternatively, in the LinReTraCe folder, execute

$ python3 setup.py install

to make them globally available.888If root access is not available, use
python3 setup.py install --user instead. This command has to be rerun after
every version update. The main program is written in MPI parallelized
Fortran95 and requires the HDF5 library.999 The LinReTraCe installation
includes an HDF5 wrapper written by one of the authors, see
https://github.com/linretrace/hdf5_wrapper for more details. There, you can
also find an installation guide for the required HDF5 library. Ensure that the
HDF5 library and LinReTraCe use the same compiler. The LinReTraCe installation
is based on Makefiles, where all necessary details, following below, are also
documented in

$ make help

The installation requires a configuration file `make_include` in the
`linretrace` folder. Here, the Fortran compiler (FC, FCDG), Fortran flags
(FFLAGS), Fortran precompiler flags (FPPFLAGS) and library paths (HDF5) are to
be provided. For a single core installation FC and FCDG are identical. Multi
core installations require the MPI compiler in FC to be compatible with FCDG.
The MPI installation is ‘activated’ via the FPPFLAGS variable `-DMPI`. An
exemplary make configuration `make_include` for a multi-core intel setup looks
like

FC = mpiifort FCDG = ifort FFLAGS = -O3 FPPFLAGS = -DMPI HDF5 =
-I/opt/hdf5-1.12.1_icc/include HDF5 += -L/opt/hdf5-1.12.1_icc/lib
-lhdf5_fortran -lhdf5hl_fortran

while a single-core gfortran setup could be

FC = gfortran FCDG = gfortran FFLAGS = -O3 HDF5 =
-I/opt/hdf5-1.12.1_gcc/include HDF5 += -L/opt/hdf5-1.12.1_gcc/lib
-lhdf5_fortran -lhdf5hl_fortran

The provided configuration can be tested via

$ make validate

which performs a small list of library and precision tests. The LinReTraCe
executable is then compiled with

$ make linretrace

which creates `bin/linretrace`.

$ make install

can be used to copy the linretrace binary to the `bin` folder in the user’s
home directory. For some compilers (notably `gfortran`) an explicit link to
the dynamic HDF5 library is necessary before executing the binary, made
available with

$ export LD_LIBRARY_PATH=/opt/hdf5-1.12.1_gcc/lib:$LD_LIBRARY_PATH

e.g., in the `.bashrc` file. A more detailed step-by-step installation guide
is provided in the repository’s `documentation/tutorial.pdf`, also available
here. Finally, a test suite can be run to verify the correct Python3/Fortran
setup and installation of the code via

$ make test

###  3.3 Flowchart

LinReTraCe calculations follow the flow chart of Fig. 1 and require two
mandatory files: an energy file and a config file. The former is a direct
result of one of the various interfaces and contains all the necessary
energies $\varepsilon(\mathbf{k},n)$, optical elements
$M^{(B),\alpha\beta(\gamma)}$, and other auxiliary data. Depending on the data
source only some of these optical elements can be generated, as, e.g., the
magnetic field optical elements require band velocities and curvatures, not
available in stand-alone WIEN2k calculations. The transport calculation itself
is configured via a config text file for which `lconfig` provides a minimal
starting point. More elaborate options, e.g., impurity states, need to be
added by hand, see Sec. 3.6. In the config file one has access to simplistic
scattering rate (and quasi-particle weight) dependencies: polynomials in
temperature. More control over these dependencies can be gained via a
scattering file where one has the option to specify all available data points
individually, see Sec. 3.5. The results of the transport calculation are then
saved in the HDF5 output file containing all the Onsager coefficients
$\mathcal{L}_{ab}^{\alpha\beta}$, $\mathcal{L}_{ab}^{B,\alpha\beta\gamma}$
among other auxiliary information, including the configuration, structure
information, etc. `lprint` provides easy access to these data containers as
well as their combinations that form the physical transport quantities. In the
next sections we go into more detail about each step, for a quick set of
instructions see the cheat sheet in Appendix A.

Figure 1: Flow chart of the LinReTraCe package: The energy file can be
generated by interfacing various electronic structure codes or modelling your
own dependencies. Combined with a config file and, _optionally_ , full
scattering dependencies, this constitutes the input of the core program
bin/linretrace. The HDF5 output file can either be accessed effortlessly by
lprint or via any external HDF5 library. Colors of boxes are as follows.
Green: possible sources of (one-particle) electronic structure input; blue:
LinReTraCe input/output files; yellow: possible ways to access data
(hdf5/screen); red: LinReTraCe executables and commands.

###  3.4 Energy File

The center piece of the LinReTraCe input is the LRTC energy file. It contains
the band energies $\varepsilon(\mathbf{k},n)$, associated optical elements
$M^{\alpha\beta}(\mathbf{k},n,m)$, $M^{B,\alpha\beta\gamma}(\mathbf{k},n,m)$
and all relevant information on the unit cell and the momentum mesh in HDF5
format. While the band-diagonal elements (energies; intra-band
$M^{\alpha\beta}(\mathbf{k},n,n)$; intra-band
$M^{B,\alpha\beta\gamma}(\mathbf{k},n,n)$) are saved as fully continuous
datasets, we explicitly separate the off-diagonal elements into a tree
structure. Each momentum point then is represented by an HDF5 group whose
datasets contain only its personal directional and band dependencies. This is
done so that the (momentum parallelized) Fortran program is able to load in
the elements in a more efficient manner. The overall tree structure of these
files then looks as follows (‘…’ signals more datasets in these groups).

/.bands Group /.bands/energyBandMax Dataset ... /.kmesh Group /.kmesh/nkp
Dataset ... /.unitcell Group /.unitcell/volume Dataset ... /energies Dataset
/kPoint Group /kPoint/0000000001 Group /kPoint/0000000001/moments Dataset
/kPoint/0000000001/momentsBfield Dataset ... /momentsDiagonal Dataset
/momentsDiagonalBfield Dataset

All the information that can be accessed in these pre-processed files can be
listed with our multi-purpose interface `lprint`. Point the script at a
specific input file and execute

$ lprint <LRTCinput file> list

We support the output of a structure information summary, the (spin-resolved)
density of states and high-symmetry band paths via

$ lprint <LRTCinput file> info $ lprint [-p] <LRTCinput file> dos $ lprint
[-p] <LRTCinput file> path

Where the `-p` (`--plot`) flag results in a graphical on-screen output. The
path option relies on the `ase` library [75], which finds the high-symmetry
points of the internal unit cell and displays them on-screen. Then, a custom
$\mathbf{k}$-path can be selected and plotted with

$ lprint [-p] <LRTCinput file> path GKMG # Gamma-K-M-Gamma

####  3.4.1 Density functional theory codes

Natively, we provide full support for the density functional theory codes
WIEN2k[21] and VASP [22], and include connections to the BoltzTraP2 [24]
library. In its basic form, to generate the LRTC energy file, the interface
`ldft` simply has to be pointed at the folder of the (converged) DFT
calculation

$ ldft <DFT folder> [--output lrtc-dft.hdf5]

This automatically detects if a valid calculation is present, initializes the
according interface and detects the type of calculation that was performed
(spin polarized, unpolarized, with or without spin-orbit-coupling).The
`--output` flag sets the output file name, if none is provided we default to
”`lrtc-dft.hdf5`”.

##### WIEN2k (version $\geq 18.x$).

We rely on the following files to be present101010 LinReTraCe assumes a serial
WIEN2k run to be present. In case a parallelized calculation was performed,
WIEN2K’s join_vectorfiles tool must be used to combine the energy files.

  * •

`case.scf` (number of electrons) [fallback: `case.in2(c)`]111111If the
converged WIEN2k calculation was saved with save_lapw the case.scf file will
not be present. In this case we fall back to case.in2 / case.in2c to retrieve
the number of electrons. However, when using the BolzTraP2 interface, case.scf
must be present, so you will have to copy the saved scf-file.

  * •

`case.struct` (atoms, symmetry operations)

  * •

`case.klist` (k-points, multiplicities)

  * •

`case.energy` (band structure)

The interface to the optical matrix elements (generated with `x optic` [70])
is done via the `case.symmat` files. The energy and optic file endings vary
depending on the type of the calculation performed, which we detect
automatically. Please note that `ldft` expects _exactly_ 3, 6 or 9 columns in
the symmat files. These correspond to the following `case.inop` configurations

  * •

3: `Re xx`, `Re yy`, `Re zz`

  * •

6: `Re xx`, `Re yy`, `Re zz`, `Re xy`, `Re xz`, `Re yz`,

  * •

9: `Re xx`, `Re yy`, `Re zz`, `Re xy`, `Re xz`, `Re yz`, `Im xy`, `Im xz`, `Im
yz`

The read-in of the optical elements $M$ is activated via the `--optic` flag

$ ldft <WIEN2k folder> \--optic

Alternatively, Peierls velocities can be used for $M^{\alpha\beta}$ via
`BoltzTraP2`. The $M^{B,\alpha\beta\gamma}$ for magnetic quantities are only
provided via the Peierls approach, see below.

##### VASP (version $\geq 5.x$).

Identically, simply point `ldft` to the VASP calculation folder where the file
`vasprun.xml` must be present:

$ ldft <VASP folder>

If, as in VASP, no optical elements $M^{\alpha\beta}$ are provided by the
electronic structure package, the interface to `BoltzTraP2` is required to
generate the necessary input for transport calculations. For the magnetic
elements, $M^{B,\alpha\beta\gamma}$, `BoltzTraP2` is always needed in
conjunction with `ldft`.

##### BoltzTraP2 interpolation.

We employ the Python3 library functionality of BoltzTrap2 to interpolate the
band structure and construct derivatives and curvatures. These are then used
to build the optical matrix elements according to Eq. (55-56). The
interpolation is activated via

$ ldft <DFT folder> \--interp [m]

The `--interp` optionally takes an additional interpolation parameter
corresponding to the `BoltzTraP2` `-m` option in the `btp2 interpolate`
command. The default value is $3$, meaning for every $\mathbf{k}$-point of the
input, `BoltzTraP2` tries to sample 3 additional irreducible
$\mathbf{k}$-points to generate the coefficients of the trigonometric
polynomial representing the band structure. If the interpolation was done on
an irreducible grid, we properly symmetrize the constructed optical elements
$M^{\alpha\beta}$ and $M^{B,\alpha\beta\gamma}$, see Sec. 4.1. Since VASP does
not provide the symmetry matrices at default settings, we extract the space
group via the `ase` library [75]. For complicated structures and symmetry-
broken calculations, we interactively ask the user to confirm the detected
space group number or to assign the correct one.

##### Band truncation.

For all input combinations we provide the possibility to truncate the energy
levels around the Fermi level in order to reduce computational costs.

$ ldft <DFT folder> \--trunc -5 +5

for example restricts the energy interval from $-5$ to $+5$eV around the DFT
Fermi level and truncates every band that is fully outside of that interval.
Please note that, contrary to, e.g., WIEN2k, we set the chemical potential of
a fully gapped system to the center-point of the fundamental gap, i.e.,
$\mu_{\mathrm{DFT}}=[\min(E_{\mathrm{valence}})+\max(E_{\mathrm{conduction}})]/2$.

####  3.4.2 Wannier90

Wannier90 [25], a program to calculate maximally localized Wannier
functions[5], is based on minimizing the total spread of the Wannier function
in real space. The output we are interested in are the real-space hopping
parameters, $H_{ll^{\prime}}(\hbox{{R}})$, which we extract from the files

  * •

`case.nnkp` (lattice vectors, k-points, projections)

  * •

`case_hr.dat` (hopping $H_{ll^{\prime}}(R)$)

  * •

`case.wout` (units)

The multi-orbital Hamiltonian in reciprocal space and its derivatives
(velocities $v$ and curvatures $c$) are then constructed from Fourier
transforms

$\displaystyle H_{ll^{\prime}}(\mathbf{k})$
$\displaystyle=\sum_{\mathbf{R}}e^{i\mathbf{k}\cdot\mathbf{R}}H_{ll^{\prime}}(\mathbf{R})$
(58) $\displaystyle v_{\mathbf{k}ll^{\prime}}^{\alpha}$
$\displaystyle=i\sum_{\mathbf{R}}R^{\alpha}e^{i\mathbf{k}\cdot\mathbf{R}}H_{ll^{\prime}}(\mathbf{R})$
(59) $\displaystyle c_{\mathbf{k}ll^{\prime}}^{\alpha\beta}$
$\displaystyle=-\sum_{\mathbf{R}}R^{\alpha}R^{\beta}e^{i\mathbf{k}\cdot\mathbf{R}}H_{ll^{\prime}}(\mathbf{R})$
(60)

where $R^{\alpha}$ is the Cartesian component $\alpha\in\\{x,y,z\\}$ of the
unit-cell vector R. For unit-cells with more than one atom, the velocities
include an extra term from the generalized Peierls approach[43], see Eq. (57).
We diagonalize the Hamiltonian matrices to go into the Kohn-Sham basis via

$\varepsilon(\mathbf{k})=U^{-1}(\mathbf{k})H(\mathbf{k})U(\mathbf{k}).$ (61)

Applying the same unitary matrices onto the velocity and curvature matrices
allows us to construct the transition matrices in the band basis. This
procedure is more justified[43] than performing the momentum derivatives in
the Kohn-Sham basis (see above). To interface a Wannier90 calculation simply
point `lwann` at the corresponding interface

$ lwann <Wannier90 folder> [--charge <N>] [--output lrtc-wann.hdf5]

where we try to extract the required charge `N` in the wannierized bands
ourselves by rounding the calculated charge (at Wannier90’s $\mu=0$) to its
nearest integer value. This can be adjusted by providing the charge yourself
with the optional `--charge` flag. The default output name `lrtc-wann.hdf5`
can be adjusted with the `--output` flag. Since wannierizations are usually
calculated on a restricted momentum grid, we provide options to refine the
momentum grid directly:

$ lwann <Wannier90 folder> \--kmesh <nkx> <nky> <nkz>

can be used to increase the reducible grid, where the new mesh has to conform
to the old mesh’s symmetry. Since reducible calculation can be costly we also
provide an additional WIEN2k sub-interface that allows for setting up large
irreducible meshes:

$ lwann <Wannier90 folder> \--wien2k

Here the additional files

  * •

`case.klist`

  * •

`case.struct`

must be present in the same folder with the same file prefix (`case`) as used
by Wannier90, as is standard when using wien2wannier[77]. `Case.klist`
(generated with WIEN2k via `x kgen`) provides a new (irreducible) momentum
mesh with associated multiplicities while the symmetry operations in
`case.struct` are used to symmetrize the calculated velocities and curvatures,
see Sec. 4.1.

####  3.4.3 Tight binding models

In the same vein we provide an interface to generate input data from arbitrary
tight-binding parameter sets. To this end we internally interface `spglib`
[76] to generate the unit cell. Akin to Wannier90 we require the primitive
lattice vectors and the hopping parameters. In addition, a list of
(in)equivalent atoms inside the unit cell is required to determine the unit-
cell symmetries. With this information the tight binding input can be created
via

$ ltb <tb file> <nkx> <nky> <nkz> <filling>

where we additionally provide the desired number of $\mathbf{k}$-points in
each direction and an initial band filling.121212 For insulators and
semiconductors, the nominal filling should be given. The generated energy file
can then be used also for doping studies, as the number of carriers (and also
the band gap) can be modified via the config.lrtc file, see Appendix B.  The
underlying equations are identical to Sec. 3.4.2. An exemplary tight binding
file looks as follows

begin hopping # a1 a2 a3 orb1 orb2 hopping.real [hopping.imag] 0 0 0 1 1 0.3 #
on site energy +1 0 0 1 1 1.0 # nearest neighbor hopping 0 +1 0 1 1 1.0 -1 0 0
1 1 1.0 0 -1 0 1 1 1.0 end hopping begin atoms # sort rx ry rz 1 0 0 0 #
fractional coordinates end atoms begin real_lattice # x y z 5 0 0 # a1 lattice
vector in units of Angstroem 0 5 0 # a2 0 0 1 # a3 end real_lattice

A number of different structures are saved as templates in the LinReTraCe
repository. The hopping parameters are specified within the `begin hopping`
and `end hopping` markers where each row represents the directional hopping
amplitude along the lattice vectors (first three columns), the associated
orbital indices (next two columns, numbering starts at $1$) and finally the
hopping amplitude (in units eV) in the last (two) column(s). The optional, 7th
column allows for imaginary contributions to inter-orbital
(`orb1`$\neq$`orb2`) hoppings. Important to mention here is that, contrary to
Wannier90, we use the sign convention commonly used in the strongly correlated
electron community131313Hoppings are positive, with the gain of kinetic energy
accounted for by a global minus sign, e.g.,
$\epsilon(\hbox{{k}})=-2t\sum_{\alpha=1,2,3}\cos(\hbox{{k}}_{\alpha})$ with
$t>0$.

$\displaystyle H(\mathbf{k})$
$\displaystyle=-\sum_{\mathbf{R}}e^{i\mathbf{k}\cdot\mathbf{R}}(1-2\delta_{\mathbf{R},\mathbf{0}}\delta_{l,l^{\prime}})H_{ll^{\prime}}(\mathbf{R}).$
(62)

The fractional atomic positions in the unit cell are listed within the `begin
atoms` and `end atoms` markers where each row represents one atom. The first
column counts the atomic ‘species’ starting at one, whereas the next three
columns describe the fractional position within the unit cell with respect to
the provided lattice vectors, i.e., all directional values are in the range
$0\leq r_{i}<1$.141414see https://spglib.github.io/spglib/python-spglib.html
for details (last accessed 12.06.2022). The atomic positions determine the
symmetries of the unit cell. Note that if the symmetry of the real-space
hopping Hamiltonian is lower than the crystal-symmetry inferred from the
atoms, the construction of the irreducible $\mathbf{k}$-mesh will fail. In
that case, ”virtual” atoms can be placed to manually break excess unit-cell
symmetries. Finally, the lattice vectors in units of Ångstrom are saved as
rows in between the `begin real_lattice` and `end real_lattice` markers.

Currently, the ltb interface only supports spin-unpolarized tight-binding
files. A generalization is straight forward. The filling argument is needed to
pre-compute the chemical potential. Note that the counting is such that a
fully filled orbital hosts $2$ electrons. Again, for doping studies, the
number of electrons can later be adjusted via the config file.

####  3.4.4 Generic interface

`linterface` contains the class `StructureFromArrays`, derived from the
abstract base class `ElectronicStructure`. An object is instantiated with the
number of k-points in the three reciprocal lattice directions, list of real
space lattice vectors, and the total charge in the system. After loading in
the multiplicity, the energies, and optical elements in form of lists or
`numpy` arrays with the corresponding `loadData` method, the output method
`outputData` takes care of calculating the chemical potential, the setting of
required flags and the file output itself. Since this interface is agnostic to
the origin of the input, it can be used to interface any data source,
including other density functional codes, dynamical mean field theory codes
and so on. If developers write an interface to their electronic structure
code, we will be happy to include it in the package.

###  3.5 Scattering File

For typical runs, information on the scattering rate and quasi-particle
weights is set-up in the config file, see below. However, to include the full
state- and momentum-dependence of the scattering amplitude $\Gamma$, the
quasi-particle weight $Z$ and band-shifts $\Re\Sigma$, an otherwise optional
LRTC scattering file has to be created. Therewith, $\Gamma$, $Z$, and
$\Re\Sigma$ extracted from, e.g., a many-body calculation can be directly used
to compute transport properties. To keep this route as generic as possible the
user is required to interact with Python3 code, where at its core a scattering
object is instantiated. First, the LRTC energy file from the previous section
is read in to initialize states and momenta. Then, the only steps necessary
are the definition of the calculation axis (chemical potential or temperature
scan), the load in of user defined `numpy` arrays that describe the scattering
dependencies, and the final output. The work flow is as follows:

  1. 1.

Copy `lscat_template` from the installation folder into your working
direction.

  2. 2.

Insert the linretrace folder into the system path.151515If the Python package
was installed via python3 setup.py install this step is not necessary.

  3. 3.

Reference to correct energy file.

  4. 4.

Define calculation axis ($\mu$ or $T$-scan).

  5. 5.

Define scattering rates as a numpy array.
Optionally: define quasi particle weights and/or band shifts.

  6. 6.

Execute script to generate LRTC scattering file:

$ python3 lscat_template

A minimalistic script can look as simple as

import sys import numpy as np sys.path.insert(0,’/home/user/linretrace’) #
Step 1: installation folder from scattering.fullscattering import
FullScattering scatobj = FullScattering(’energies.hdf5’) # Step 2: # reference
LRTC energy file # main dependencies: spins, kpoints, bands spins, nkp, nbands
= scatobj.getDependencies() kgrid = scatobj.getMomentumGrid() # [ nkp, 3 ]
energies = scatobj.getEnergies() # [ spins, nkp, nbands ] mudft =
scatobj.mudft # chemical potential # Step 3: define temperature range in
Kelvin nt = 100 scatobj.defineTemperatures(tmin = 10, \ tmax = 300, \ nt = nt,
\ tlog = True) # Step 4: Create array that defines scattering rates in eV
gamma = np.zeros((nt, spins, nkp, nbands), dtype=np.float64) # Gamma
gamma[0,...] = (energies-mudft)**2. / 10000. gamma[1:,...] = gamma[0,...] #
Optional Step 5: Create arrays that define Z and ReSigma(0) qpweight =
np.ones_like(gamma, dtype=np.float64) # Z bandshift = np.zeros_like(gamma,
dtype=np.float64) # ReSigma scatobj.defineScatteringRates(gamma, qpweight,
bandshift) scatobj.createOutput(’scattering_file.hdf5’)

How to handle the other dependencies is illustrated in code snippets in the
provided `lscat` and `lscat_template`. In case scattering-rates and quasi-
particle weights do not depend on band and momentum, and follow a simple
temperature dependence, they can be more conveniently specified in the
LinReTraCe configuration file, see the next section.

###  3.6 Configuring and running LinReTraCe

LinReTraCe is configured via a free format text configuration file. A
minimalist starting point for this config file can be generated with
`lconfig`. Here, through interactive questioning a basic config file is
generated and saved as `config.lrtc`. From there, more elaborate options can
be added manually. For more advanced settings, see Appendix B, for a full
documentation of all possibilities, see the configuration specification in
`documentation/configspec`.

In its basic form the configuration file sets the run mode and defines which
quantities should be calculated and at which precision. LinReTraCe can either
scan through a range of temperatures for a given number of electrons (`RunMode
= temp`) or scan through various chemical potentials (`RunMode = mu`) at fixed
temperature. Further, one is also able to modify some information from the
energy file on-the-fly, like the bandgap or the number of electrons in the
system, e.g., for doping studies.

The run modes are then configured in their own respective section of the
config file. Without a scattering file, the temperature mode is configured by
specifying a temperature range and provides the option to include impurity
states as well as homogeneous doping, see Sec. 4.2. Scattering rates can be
specified for a simple polynomial temperature dependence
($\Gamma(T)=\Gamma_{0}+\Gamma_{1}T+\Gamma_{2}T^{2}+\cdots$) when acting on all
bands and momenta equally. For the inclusion of arbitrary scattering rates, a
scattering file needs to be created, see Section 3.5 and Appendix B.3. The
chemical potential mode instead works at a fixed temperature and only requires
information on the range over which the scan is performed and which scattering
rate and quasi-particle renormalization to use.

For a temperature scan, the output of `lconfig` looks as follows

[General] # input/output configuration RunMode = temp # scan through
temperatures EnergyFile = lrtc-dft.hdf5 # any energy file from Sec. 3.3
OutputFile = output.hdf5 # output file name BFieldMode = T # calculate L11B
L12B L22B Interband = F # T/F enable/disable inter-band Intraband = T # enable
intra-band Boltzmann = T # Kubo AND Boltzmann responses QuadResponse = T #
kernel evaluation: quad precision FermiOccupation = F # digamma function as
occupation [TempMode] # temperature mode sub configuration [[Scattering]]
TMinimum = 100.0 # temperature range [K] TMaximum = 700.0 TPoints = 100
TLogarithmic = T # logarithmic temperature steps ScatteringCoefficients = 1e-5
0 1e-8 # G0 G1 G2 ... # G(T) = G0 + G1 * T + G2 * T**2 ...
QuasiParticleCoefficients = 1 # renormalization Z0 ... # Z(T) = Z0 ...

The configuration of the chemical potential mode looks like

[General] RunMode = mu ... [MuMode] # mu mode sub configuration [[Scattering]]
Temperature = 300.0 # fixed temperature [K] MuMinimum = -5.0 # chemical
potential range [eV] MuMaximum = +5.0 # with respect to mu_DFT MuPoints = 100
ScatteringRate = 1e-5 # fixed Gamma [eV] QuasiParticleWeight = 1.0 # fixed Z
[0 < Z <= 1]

Once the input and configuration has been prepared, LinReTraCe is run by
executing the binary with the generated config file as argument. The MPI
installation is invoked via

$ mpirun -np <cores> bin/linretrace config.lrtc

and the single core installation via

$ bin/linretrace config.lrtc

The internal flow of this program is listed in Table 1.

| function name | file
---|---|---
initialize MPI interface | mpi_initialize | mpi_org.F90
read config file | read_config | config.f90
initialize HDF5 interface | hdf5_init | hdf5_wrapper.F90
preprocess energy file | read_preproc_energy | input.f90
_optional_ : preprocess scattering file | read_preproc_scattering_hdf5 | input.f90
or: | read_preproc_scattering_text | input.f90
_optional_ : read $\mu(T)$ | read_muT_hdf5 | input.f90
or: | read_muT_text | input.f90
distribute MPI $\mathbf{k}$-load | mpi_genkstep | mpi_org.F90
read $\varepsilon(\mathbf{k},n)$, $M(\mathbf{k},n,n)$ | read_energy | input.f90
_optional_ : read $M^{B}(\mathbf{k},n,n)$ | read_energy | input.f90
allocate data structures | allocate_response | response.F90
_optional_ : read $M(\mathbf{k},n,m)$ | read_full_optical_elements | input.f90
_optional_ : read $M^{B}(\mathbf{k},n,m)$ | read_full_magnetic_elements | input.f90
initialize HDF5 output | output_auxiliary | output.F90
LOOP: or temperature $T$ |  |
LOOP: or chemical potential $\mu$ |  |
|  _optional_ : read $\Gamma(k,n,\mu|T)$ | read_scattering_hdf5 | input.f90
|  _if necessary_ : calculate $\mu(T)$ | find_mu | root.F90
| calculate $\psi_{n}(\mathbf{k},n)$ | calc_polygamma | response.F90
| initialize data structures | initialize_response | response.F90
| LOOP: momentum $\mathbf{k}$ |  |
|  |  | response_intra_km | response.F90
|  | calculate $\mathcal{K}_{ab}M^{\alpha\beta}$, $\mathcal{K}^{B}_{ab}M^{B,\alpha\beta\gamma}$ | response_inter_km | response.F90
|  |  | response_intra_km_Q | response.F90
|  |  | $\vdots$ | response.F90
|  | save / output $\mathcal{L}_{ab}^{\alpha\beta}$, $\mathcal{L}^{B,\alpha\beta\gamma}_{ab}$ | response_h5_output | response.F90
|  |  | response_h5_output_Q | response.F90
output $\mu,T,n_{e},n_{h},\cdots$ |  |
close MPI and HDF5 interface |  |
Table 1: Internal program flow of LinReTraCe.

During execution, first an options summary, then continuous run information is
printed to the standard output. After a successful exit, the generated HDF5
output file contains the calculated Onsager coefficients as well as all
relevant config information and some auxiliary datasets, see the next section.

###  3.7 Output File

In the HDF5 output file, all the calculated Onsager coefficients are saved as
a combination of their identifier `L11`, `L12`, `L22`, `L11B`, `L12B`, `L22B`
and the type of response that was used: `intra`, `inter`, `intraBoltzmann`,
`interBoltzmann`. Non-magnetic datasets contain the momentum- and band-summed
quantities with array shape `[steps, spins, 3, 3]` and magnetic datasets
contain the momentum- and band-summed quantities with array shape `[steps,
spins, 3, 3, 3]` where the last 2 (3) dimensions refer to the Cartesian
directions $[\alpha,\beta]$ ($[\alpha,\beta,\gamma]$). Please note that we
also include the possibility to output the Onsager coefficients with either
their full dependency (k, $n$), only momentum-k-summed and only
band-$n$-summed, see Appendix B. The standard tree structure of the output
file looks as follows

/.config Group /.quantities Group /.quantities/mu Dataset
/.quantities/tempAxis Dataset ... /.scattering Group ... /.structure Group
/.structure/charge Dataset ... /.unitcell Group /.unitcell/vol Dataset ...
/L11 Group /L11/intra Group /L11/intra/sum Dataset /L12 Group /L12/intra Group
/L12/intra/sum Dataset /L22 Group /L22/intra Group /L22/intra/sum Dataset

While extracting information directly from these files is possible (necessary
for the full dependence output), it can be quite cumbersome. To this end we
provide the post-processing interface `lprint`. Simply point `lprint` to the
specific output file (or use ”latest” to search for the latest valid output
file) and execute

$ lprint <LRTCoutput file> list

to list all available physical datasets to print (`olist` to list all the
Onsager coefficients). The configuration is listed via

$ lprint <LRTCoutput file> config

and datasets can be output (text/plot) by referring to the keys listed in the
”list” option. To name a few selected datasets: the chemical potential can be
plotted via

$ lprint [-p] <LRTCoutput file> mu

Transport tensors are evaluated from the Onsager coefficients via Eqs. 11-20,
and can be extracted/plotted using the keywords from Tab. 2. `lprint` can be
supplied with Cartesian directions chained after each other, e.g., the `xx`
and `yy` entries of the conductivity tensor can be obtained via

$ lprint [-p] <LRTCoutput file> c-intra xx yy

If no direction is supplied, all combinations xx, xy, xz, yx, … are plotted.
Magnetic tensors, e.g., the Nernst coefficient require three directions, where
the third defines the Cartesian direction of the applied magnetic field, cf.
Eq. (6). E.g., for the Nernst coefficient:

$ lprint [-p] <LRTCoutput file> n-intra xyz

Spin resolved output is obtained by adding ’u’ or ’d’ to the directions (`uxx
dxx` for the up and down component of the `xx` entry, respectively).

Identical datasets of different files can be compared with the `-c`
(`--compare`) flag:

$ lprint [-p] <out1> L12-inter dzz -c <out2> <out3>

This flag in particular helps assessing the convergence of transport
observables with the number of $\mathbf{k}$-points, cf. Section 5.2.
Switching to the alternative print/plot axis, inverse temperature instead of
temperature: $T[\mathrm{K}]\rightarrow\beta[\mathrm{eV}^{-1}]$, or, carrier
concentration instead of chemical potential: $\mu[\mathrm{eV}]\rightarrow
n[\mathrm{cm}^{-3}]$, is done via `-x` (`--axis`).

$ lprint [-p] -x <LRTCoutput file> s-total xx

The full functionality is described in the help option

$ lprint --help

quantity | key | component
---|---|---
Onsager coefficients $\mathcal{L}_{ij}$ | Lij- | intra, inter
Onsager coefficients $\mathcal{L}^{B}_{ij}$ | LijB- | intra, inter
conductivity $\sigma$ | c- | intra, inter, total
resistivity $\rho$ | r- | intra, inter, total
Peltier $\Pi$ | p- | intra, inter, total
Seebeck $S$ | s- | intra, inter, total
power factor $S^{2}\sigma$ | pf- | intra, inter, total
thermal conductivity $\kappa$ | tc- | intra, inter, total
thermal resistivity $R_{\lambda}$ | tr- | intra, inter, total
Hall conductivity $\sigma^{B}$ | cb- | intra, inter, total
Hall coefficient $R_{H}$ | rh- | intra, inter, total
Nernst coefficient $\nu$ | n- | intra, inter, total
Hall mobility $\mu_{H}$ | muh- | intra, inter, total
thermal mobility $\mu_{T}$ | mut- | intra, inter, total
Table 2: Syntax of lprint: The Onsager and transport tensors are accessed by
combining the key and the component, e.g., s-intra for the intra-band Seebeck
coefficient. The ”total” contributions only exists if both intra- and inter-
band transitions are present. For the Onsager coefficients Lij(B) possible
indices are $ij\in\\{11,12,22\\}$.

##  4 Technical details

###  4.1 Matrix elements on irreducible grids

Periodic unit cells can be assigned a point group that is commonly represented
by a set of $n$ square matrices $P_{i}$ ($\mathrm{det}P_{i}=\pm 1$) that
describe all applicable symmetry operations acting in real space. In
reciprocal space this can be exploited as one can reduce the number of momenta
necessary to represent the full Brillouin zone (for a symmetry reduction
algorithm scaling linearly with the number of points see [78]). Each so-called
irreducible $\mathbf{k}$-point $\mathbf{k}_{\mathrm{irr}}$ then represents a
set of momenta $\\{\mathbf{k}_{i}\\}$ that are all connected to each other via
the transposed matrices $P_{i}^{T}$

$\mathbf{k}_{i}=P_{i}^{T}\mathbf{k}_{\mathrm{irr}}.$ (63)

The number of unique momenta generated from $\mathbf{k}_{\mathrm{irr}}$ is
called its multiplicity $m$, where each point in the set will be generated
exactly $\frac{n}{m}$ times. Naturally, the energies within this set remain
unchanged

$\forall
i,\mathbf{k}_{irr}:\,\,\varepsilon(\mathbf{k}_{irr})=\varepsilon(P_{i}^{T}\mathbf{k}_{irr})$
(64)

which then also applies to the kernel functions Eqs. (9-10). Crucially, owing
to their directionality, this does not hold for the associated optical
elements, band velocities and band curvatures, hence Eqs. (3-4) must not be
evaluated directly on the irreducible grid. Instead, matrix elements have to
be symmetrized: By averaging over all connected optical elements

$M^{\mathrm{symmetrized}}_{opt}(\mathbf{k}_{\mathrm{irr}})=\frac{1}{n}\sum^{\mathrm{n}}_{i=1}M_{opt}(P_{i}^{T}\mathbf{k}_{\mathrm{irr}})$
(65)

one is able to absorb all the required symmetry information. Please note that
these schemes require the momentum mesh to respect the same point group
symmetries as the unit cell itself, e.g., a cubic crystal structure requires
$n_{kx}=n_{ky}=n_{kz}$. If this were not the case,
$P_{i}^{T}\mathbf{k}_{\mathrm{irr}}$ generates points outside the initial
grid.161616An incommensurate reducible grid could still be reduced, by
deselecting invalid momenta. Also an exact symmetry mapping of every single
irreducible $\mathbf{k}$-point to all $m$ connected momenta would solve this
issue. Our current implementation, however, requires the user to make a
sensible choice for the grid.

While density functional theory codes like WIEN2k provide dipole matrix
elements on an irreducible grid, optical elements as listed in Sec. 2.4.2 need
to be symmetrized explicitly. Since Eq. 65 relies on information from the full
Brillouin, the symmetrization is implemented via real-space rotations. Band
velocities then transform as

$v(P_{i}^{T}\mathbf{k}_{\mathrm{irr}})=\left[K^{-1}P_{i}^{-1}K\right]v(\mathbf{k}_{\mathrm{irr}})$
(66)

whereas band curvatures $c$ transform as

$c(P_{i}^{T}\mathbf{k}_{\mathrm{irr}})=\left[K^{-1}P_{i}^{-1}K\right]c(\mathbf{k}_{\mathrm{irr}})\left[K^{-1}P_{i}^{-1}K\right]^{T}$
(67)

as they correspond to a single and twofold momentum derivative, respectively.
In the same vein the optical elements $M$ themselves transform as the
curvatures in Eq. (67). Here $K$ is the matrix formed by the reciprocal
lattice vectors (the rows of $K$). This transforms the point group matrix into
Cartesian directions and is explicitly necessary for non-orthogonal unit
cells. The generated velocities and curvatures are then combined to optical
elements via Eqs. (55-56), over which the symmetrization is performed. To
showcase the correctness of these equations, see Fig. 2 where we compare a
graphene-inspired honeycomb lattice on a reducible and corresponding
irreducible momentum grid.

Figure 2: Reducible calculations compared to calculations on an irreducible
grid with and without optical element symmetrization. Left: Density of states
of the honeycomb lattice with hopping of $t_{AB}=t_{BA}=1$eV and lattice
length $|a|=1$Å. Middle and right: Onsager coefficients
$\mathcal{L}_{11}^{xx}$ and $\mathcal{L}_{12}^{xyz}$ performed for $180\times
180$ reducible ($8191$ irreducible) $\mathbf{k}$-points at a constant chemical
potential of $\mu=0.20$eV (vertical dashed line) and scattering rate
$\Gamma=10^{-5}$eV. Only a properly symmetrized irreducible grid leads to
consistent data. Note that for this 2D system, Onsager coefficients needed to
be multiplied with the (fictitious) $c$-lattice constant to yield proper
units, as, e.g., $\mathcal{L}_{11}$ is linked to a conductance
$[\sigma^{2D}]=1/\Omega$ instead of a conductivity $[\sigma]=1/(\Omega m)$,
see Appendix C.3.

###  4.2 Chemical potential search

Determining the chemical potential is a common root finding problem: The
numerical search of $\mu$ can be represented by

$N-\sum_{\mathbf{k},n}f(\varepsilon_{\mathrm{k},n}-\mu)\overset{!}{=}0.$ (68)

For the case of no band renormalizations ($Z\equiv 1$) the occupation
$f(\varepsilon_{\mathrm{k},n}-\mu)$ is either determined from the Fermi
function

$f(\varepsilon_{\mathrm{k},n}-\mu)=f_{\mathrm{FD}}(\varepsilon_{\mathrm{k},n}-\mu)=\frac{1}{1+e^{\beta(\varepsilon_{\mathrm{k},n}-\mu)}}$
(69)

or from the lifetime-broadened spectrum Eq. (21), entailing (see Ref. [4, 2])

$f(\varepsilon_{\mathrm{k},n}-\mu)=\frac{1}{\pi}-\Im\psi\left(\frac{1}{2}+\frac{\beta}{2\pi}(\Gamma_{\mathrm{k},n}+i(\varepsilon_{\mathrm{k},n}-\mu))\right).$
(70)

In some cases, employing root-finding algorithms on Eq. (68) may lead to
problems. While metallic systems suffer mostly from too coarse momentum grids,
in gapped systems the very small number of thermally activated carriers can be
difficult to track.171717Which is why some relaxation-time approximation
Boltzmann codes typically use a fixed chemical potential instead of searching
for it. Due to the additional $\Gamma$-smearing in our formalism, the number
of active carriers increases and the problem is absent for reasonably large
scattering rates ($\Gamma\geq 10^{-6}$eV) and reasonable band gaps
($\Delta<10$eV) at all temperatures when using Eq. (70). Using the Fermi-Dirac
distribution for insulators, on the other hand, the root-finding is strongly
restricted in the temperatures that can be safely captured, irrespective of
the numerical accuracy, see double and quadruple precision calculations (blue
and green lines) in Fig. 3, respectively.

Figure 3: Chemical potentials in a gapped system determined via Eq. (68) with
double (blue) and quadruple precision (green) compared to the chemical
potential via the reformulated Eq. (71) (red). In all three cases the Fermi
distribution is used. For the given band gap $\Delta=0.1$eV the ‘standard’
root finding breaks down below $T=40$K and $T=10$K as the floating point
accuracy is exhausted at double and quadruple precision, respectively. The
reformulated problem on the other hand is stable down to $T=0.05$K (below
plotted range), see text.

In order to circumvent this problem, the root-finding problem can be
reformulated to

$\underbrace{\sum_{\mathbf{k},n\geq\mathrm{CB}}f(\varepsilon_{\mathrm{k},n}-\mu)}_{\mathrm{activated\;electrons}}-\underbrace{\sum_{\mathbf{k},n\leq\mathrm{VB}}f(-(\varepsilon_{\mathrm{k},n}-\mu))}_{\mathrm{activated\;holes}}\overset{!}{=}0$
(71)

in fully gapped systems with $\mu$ inside the gap: The chemical potential is
determined by balancing the electrons in the conduction bands with the holes
in valence bands. As a consequence one is not limited by machine precision
anymore and can exploit the full floating point range. Nonetheless, due to
finite bit length, temperatures are still bounded. The lowest possible
achievable temperature corresponds to resolving density contributions down to
the smallest positive number representable in quadruple precision:
$2^{-16494}$. If the occupation is determined via the Fermi function and the
chemical potential is in the middle of the band gap $\Delta$, it follows from

$\frac{1}{e^{\beta\varepsilon}+1}\approx
e^{-\beta\varepsilon}=e^{-\beta\frac{\Delta}{2}}$

that the lowest temperature bound is

$T^{\mu}_{\mathrm{bound}}[\mathrm{K}]=\frac{\Delta[\mathrm{eV}]}{2\ln(2)16494k_{B}}\approx
0.5\Delta[\mathrm{eV}].$ (72)

The chemical potential determined via this refined root-finding problem is
also illustrated in Fig. 3.

####  4.2.1 Impurity levels

LinReTraCe allows the inclusion of passive impurity states: Neglecting
explicit contributions to the transport functions, the charge of these extra
states affects the transport data merely through the position of the chemical
potential. The above refinement algorithm is especially important when
employing impurity states inside the gap.181818Without them, one could
approximate the chemical potential obtained via the Fermi function, by
extrapolating $\mu$ towards the band-gap mid-point at $T=0$[4]. Including
impurity states is straight forward:

  1. 1.

calculate (intrinsic) electron occupation $f(\mu)$ according to Eq. (69) or
Eq. (70)

  2. 2.

calculate (extrinsic) impurity contribution $n_{\mathrm{imp}}$ according to
Eq. (73) or Eq. (74)

  3. 3.

calculate total occupation:
$n(\mu)=f(\mu)-n^{D}_{\mathrm{imp}}+n^{A}_{\mathrm{imp}}$ (given a donor
(acceptor) level, the chemical potential has to increase (decrease) to
compensate)

  4. 4.

determine $\mu$ according to Eq. (68) or Eq. (71)191919The reformulated root-
finding problem of Eq. (71) becomes:
$\sum_{\mathbf{k},n\in\mathrm{CB}}f(\varepsilon_{\mathrm{k},n}-\mu)-\sum_{\mathbf{k},n\in\mathrm{VB}}f(-(\varepsilon_{\mathrm{k},n}-\mu))-n_{\mathrm{imp}}^{D}+n_{\mathrm{imp}}^{A}\overset{!}{=}0.$

Here the impurity contribution differs between donor (D) and acceptor (A)
levels

$\displaystyle n^{D}_{\mathrm{imp}}$
$\displaystyle=\frac{\rho_{D}}{1+ge^{\beta(\mu-E^{D}_{\mathrm{imp}})}}$ (73)
$\displaystyle n^{A}_{\mathrm{imp}}$
$\displaystyle=\frac{\rho_{A}}{1+ge^{\beta(E^{A}_{\mathrm{imp}}-\mu)}}$ (74)

where $\rho^{D/A}$ is the impurity density, $g$ is the impurity degeneracy and
$E_{\mathrm{imp}}$ is the impurity position. For a donor level in the vicinity
of the conduction band edge the chemical potential has to increase to
compensate for the additional impurity occupation. Fig. 4 illustrates the
effect for varying impurity densities. Using the Fermi distribution for
intrinsic states, the chemical potential approaches the center point between
the impurity level $E_{\mathrm{imp}}^{D}$ and the closest conduction state
(instead of the band-gap middle point, realized for $\rho_{D}=0$). From a
transport perspective the effective band gap transforms from
$\Delta=E_{c}-E_{v}$ at high temperatures to $\Delta=E_{c}-E_{\mathrm{imp}}$
at low temperatures. The transition temperature is controlled by the impurity
density and partially by the degeneracy of the impurity level. In addition to
single impurity levels LinReTraCe also offers finite-size impurity bands with
various shapes including box (constant), half-circle, and squared sine, see
Appendix B.2.

Figure 4: Comparison of the chemical potential for different impurity
densities $\rho_{D}$ for a donor level located at $E_{\mathrm{imp}}^{D}$
inside a semiconducting gap of 100meV. Intrinsic states were chosen to be
described by the Fermi distribution ($\Gamma=0$).

####  4.2.2 Homogeneous doping

Besides explicit impurity states, LinReTraCe also supports generic doping.
Contrary to impurity levels where the chemical potential converges for
$T\rightarrow 0$ to a point inside the gap, any global doping forces the
chemical potential eventually to move outside the gap. Leading up to this, the
root finding works identically as in Sec. 4.2.1. Here the ‘impurity
contribution‘ is simply the doping which, now, is not affected by temperature
and the position of the chemical potential. Please note that, technically,
this kind of doping is more nuanced than a simple change of the total electron
occupation. For the refinement algorithm to work inside the gap, an underlying
integer filling is mandatory. Thus instead of changing the filling, we
introduced an explicit `Doping` keyword in the config, see Appendix B.

In Fig. 5 we showcase the same underlying band-structure as in Sec. 4.2.1 with
various electron doping levels. Note the differences to Fig. 4: For the
largest shown doping level ($\delta_{e}=10^{-6}$) the crossing into the
conduction band (shaded gray) already happens at around $T=100$K.

Figure 5: Comparison of the chemical potential for different levels of
electron doping $\delta_{e}$ in a semiconductor. This example uses the Fermi
distribution ($\Gamma=0$) to find the chemical potential.

###  4.3 Polygamma evaluation

While all Fermi-function related equations can be implemented efficiently with
native Fortran functions up to quadruple precision, the evaluation of the
digamma and polygamma functions is more delicate. The digamma function is the
derivative of the natural logarithm of the Gamma function whose series
representation is closely related to the harmonic numbers

$\displaystyle\psi(z)$
$\displaystyle=\frac{d}{dz}\ln\Gamma(z)=\frac{\Gamma^{\prime}(z)}{\Gamma(z)}$
(75) $\displaystyle\psi(z)$
$\displaystyle=-\gamma+\sum_{n=1}^{\infty}\left(\frac{1}{n}-\frac{1}{n+z-1}\right)$
(76)

where $\gamma$ is the Euler–Mascheroni constant. The polygamma function
$\psi_{m}$ ($m>0$) is the $m^{\mathrm{th}}$ derivative of the digamma function
whose series expansion follows directly from the differentiation

$\displaystyle\psi_{m}(z)$
$\displaystyle=\frac{d^{m}}{dz^{m}}\psi(z)\quad\quad m\in\mathbb{N}^{+}$ (77)
$\displaystyle\psi_{m}(z)$
$\displaystyle=(-1)^{m+1}m!\sum_{n=1}^{\infty}\frac{1}{(n+z-1)^{m+1}}.$ (78)

Eqs. (76,78) are used extensively in our analytic derivation of the transport
kernels [2]. Numerical implementations instead use different expressions
depending on the location of the argument $z$ in the complex plane. In
particular, for large $\Re(z)$, an asymptotic Bernoulli expansion makes the
evaluation very efficient[67]:

$\psi(z)=\ln(z)-(2z)^{-1}-\sum_{k=1}^{n-1}\frac{B_{2k}}{2k}z^{-2k}+\mathcal{O}(z^{-2n})\qquad|\arg(z)|\leq\pi-\epsilon,\,\epsilon>0$
(79)

with the Bernoulli numbers $B_{2k}$. In LinReTraCe, we adapt the cernlib [79]
Fortran routine `wpsipg` (v1.2) for polygamma functions with complex argument,
originally described by Kölbig [80]. We upgraded the routine to quadruple
precision and commensurately increased the Bernoulli expansion order (to up to
$k=16$).202020 At $T=5.8$K and $\Gamma=0$ this setting reproduces the
analytical Fermi function $f_{FD}(\omega)$ with an error of less than
$10^{-27}$ for any $\omega$.

###  4.4 Code scaling

As the code is meant to be a hybrid—designed to solve models and large
realistic materials—the scaling behavior of the parallelization is important:
LinReTraCe is parallelized over the number of Brillouin zone momenta $n_{k}$
(see Sec. 5.2 for a momentum convergence test). The resulting scaling is
illustrated in Fig. 6. As expected, a purely linear behavior (runtime
$t\propto n_{k}$) emerges.212121The calculations contain $100$ temperature
steps and were performed on a single node of the Vienna Scientific Cluster 4.
Each node consists of two sockets with one Intel Xeon Platinum 8174 Processor
@ 3.1 GHz (formerly Skylake) on each socket, leading to $48$ physical cores.
We used the Intel Fortran compiler ifort and mpiifort (2020 release) for the
single core and multi core installation, respectively, with -O3 optimization.
While a single core installation experiences almost no overhead, the MPI
implementation requires roughly $200$ data points per core to become
efficient. The quadruple precision evaluation of kernels roughly doubles the
runtime.

In intra-band calculations the runtime will necessarily further scale linearly
with the number of bands, while the inter-band portion of the code will scale
with the square of the number of bands, as each band permutation must be
evaluated. In normal circumstances, i.e., for small primitive unit cells, the
number of bands is usually limited to between $\mathcal{O}(10)$ and
$\mathcal{O}(100)$ where necessarily the momentum mesh must remain dense. On
the contrary, in super cells one reallocates the computational load from the
momentum mesh to the number of bands. The increased unit cell size requires
fewer momenta to achieve convergence, which is counterbalanced by an increased
number of atoms in the cell and thus an increased number of bands. To combat
this intrinsic scaling problem one is encouraged to truncate the number of
bands to a specific energy range around the Fermi level, see Section 3.4.1.
Super cells especially benefit from this drastic decrease of the number of
bands while maintaining numerical accuracy.

Figure 6: Runtime $t$ of a one-band model over a wide range of momentum grids
with $n_{k}$ points. Scaling is linear in $n_{k}$ for single-core and MPI
runs. Communication and input overheads delay linear scaling for MPI. ”quad”
and ”double” indicate the used internal precision.

##  5 Test cases

In this section we are going to apply the previously described transport
methodology. Our aim is to highlight some of the features implemented in
LinReTraCe on systems that cover a wide range of phenomena, without, however,
exhausting the code’s full functionality. We start off with a documented
benchmark of silicon in Sec. 5.1, after which we consider various prototypical
behavior of simple models in Sec. 5.2, Sec. 5.3, and 5.4 before transitioning
to further realistic crystal structures in Sec. 5.5 and Sec. 5.6.

Input and configuration files of the silicon benchmark are part of the test
suite (see Sec. 3.2), while the data to reproduce all other examples can be
found in the repository github.com/linretrace/linretrace_examples.222222For
storage reasons, the hosted HDF5 input files may contain data on a coarser
momentum mesh than used for the results shown here. For the model systems,
however, the user can directly reproduce our results by generating their own
HDF5 input files using the number of $\mathbf{k}$-points specified in the
text. The same is possible for the material test cases, but requires
performing ones own DFT calculations.

###  5.1 Silicon: Benchmark to BoltzTraP

We start off with a simple benchmark test of silicon. As input we employ a
Wien2K LDA calculation with $36\times 36\times 36$ k-points resulting, as
expected, in an underestimated band gap $\Delta_{\mathrm{LDA}}=0.48$eV,
compared to the experimental value of $\Delta_{\mathrm{exp}}=1.12$eV at room
temperature. The benchmark test is performed with the following commands, see
Sec. 3 for a detailed documentation:

$ ldft testsuite/tests/Si --interp 3 --output Si_peierls.hdf5 $ lconfig
INTERACTIVE GENERATION OF LRTC CONFIG FILES energy file: Si_peierls.hdf5
output file: Si_muscan.hdf5 calculate intra-band quantities [y,n]: y calculate
inter-band quantities [y,n]: n calculate boltzmann quantities [y,n]: y
calculate magnetic field quantities [y,n]: n calculate quadruple precision
[y,n]: y run mode (temp [1], mu [2]): 2 use Scattering file [y,n]: n MU-MODE
configuration: temperature [K]: 300 chemical potential starting point: -1.5
chemical potential ending point: +1.5 chemical potential points: 300
scattering rate: 1e-3 quasiparticle weight: 1 config.lrtc successfully
created. $ mpirun -np 1 bin/linretrace config.lrtc $ lprint -p Si_muscan.hdf5
c-intra xx yy zz # conductivities $ lprint -p Si_muscan.hdf5 s-intra xx yy zz
# Seebeck coefficents

Following the execution of LinReTraCe, a verbose summary of the configuration
and calculation progress is printed to the standard output of the terminal.
The graphical output of the two lprint commands are shown in Fig. 7.

(a) conductivity $\sigma$

(b) Seebeck coefficient $S$

Figure 7: Default lprint output via matplotlib of a $\mu$-scan calculation for
silicon. Please note that the chemical potential $\mu_{\mathrm{DFT}}$ is
centered inside the gap. As the system is three dimensional, the conductivity
is in units of $\left[\frac{1}{\Omega\mathrm{m}}\right]$ and the Seebeck
coefficient in $\left[\frac{\mathrm{V}}{\mathrm{K}}\right]$, see Appendix C
for a full dimensional analysis.

The transport data can now be quantitatively compared to BoltzTraP2 for which
the following commands were used:

$ btp2 interpolate testsuite/tests/Si -m 3 # identical to --interp 3 $ btp2
integrate interpolation.bt2 300 # T = 300K

The relevant output data of BoltzTraP2 ($\sigma/\tau_{0}$, $S$) is saved in
text format. While Seebeck coefficients can be compared directly,
$\sigma/\tau_{0}$ must be multiplied with the relaxation time to be comparable
to the physical conductivity of LinReTraCe

$\sigma=\frac{\sigma}{\tau_{0}}\frac{\hbar}{2\Gamma}.$

We find reasonably good agreement between the two codes, see Fig. 8. This is
expected as the chosen temperature places the system well into the regime in
which the semi-classical approximation holds. The transport data of BoltzTraP2
is smoother than LinReTraCe’s as the former suffers less from the coarse
momentum grid. The reason is that its transport tensors are computed via a
numerically broadened transport density of states [24]. LinReTrace on the
other hand evaluates the sum over all momenta without additional numerical
(not physical) broadening, see Eqs. (3-4). This may lead to more spiky
transport functions, but will provide a less-biased description, in particular
near the band-edges.

Figure 8: (a) Conductivity and (b) Seebeck coefficient of the silicon
benchmark calculation, comparing LinReTraCe and BoltzTraP2. The temperature
was set to $T=300$K and the scattering rate to $\Gamma=10^{-3}$eV. The energy
range of the DFT bandgap is highlighted in gray.

###  5.2 One-band metal in two dimensions

We consider the simple electronic structure

$\epsilon^{0}_{\mathbf{k}n}=-2t\sum_{\alpha=x,y}\cos\left(k_{\alpha}\right)$
(80)

with hopping amplitude $t=0.25$eV, lattice spacing $a_{x}=a_{y}=1$Å and
optical elements determined with the Peierls approximation. In order to
introduce particle-hole asymmetry (to generate a finite thermopower) we set
the system’s charge to $N=1.2$ and determine the chemical potential with the
digamma occupation Eq. (70). Figure 9 shows the conductivity and Seebeck
coefficient for a temperature-independent scattering rate $\Gamma=10^{-4}$eV,
no quasi-particle renormalization $Z=1$, and various momentum grids. Once
$\mathbf{k}$-convergence is reached, the conductivity is essentially
temperature independent. Then, as expected from the Mott formula [81], a
linear-in-$T$ behaviour appears in the Seebeck coefficient. As the system is
above half-filling, the carriers are of hole-type, and the Seebeck coefficient
is positive. Unsurprisingly, a dense momentum grid is required to reach
convergence at low temperatures. For the largest momentum grid employed here
($n_{k}=n_{k_{x}}\times n_{k_{y}}=5000\times 5000$) the results are converged
down to approximately $T_{c}=6$K. Please note that due to the chemical
potential search, discrepancies for coarse $\mathbf{k}$-meshes are the results
of a mixture of chemical potential and kernel sampling errors. In all, even if
the $\mathbf{k}$-convergence in LinReTraCe is more stable than in Boltzmann
codes, a thorough check of the Brillouin-zone discretization is mandatory—at
least when the response is metallic.

Figure 9: (a) Resistivity and (b) Seebeck coefficient for the two-dimensional
electronic structure of Eq. (80) with $\Gamma=10^{-4}$eV, system charge
$N=1.2$ and various numbers of $\mathbf{k}$-points $n_{k}$. The convergence
temperature $T_{c}$ (inset in (a)) scales as $T_{c}\propto
n_{k}^{-0.5}=n_{k_{x}}^{-1}$ (fitted, dashed line). A dense momentum grid is
required to reach convergence (fat gray line) for reasonably low temperatures.
There, we observe an almost temperature independent resistivity and a linear-
in-$T$ behavior of the Seebeck coefficient. The identical behavior is observed
in the three dimensional equivalent where instead the convergence temperature
scales as $T_{c}\propto n_{k}^{-0.33}=n_{k_{x}}^{-1}$ (not shown).

###  5.3 Two-band insulator in three dimensions

Next, we consider the two-band electronic structure used in Ref. [2]

$\epsilon^{0}_{\mathbf{k}n}=-\sum_{\alpha=x,y,z}2t_{n}\cos\left(k_{\alpha}\right)+(-1)^{n}(6t_{n}+\Delta_{0}/2)$
(81)

with valence band hopping amplitude $t_{1}=0.25$eV and lattice spacing
$a_{x}=a_{y}=a_{z}=1$Å. Here, we use the Peierls approach in the band-basis,
see Eq. (55), to compute matrix elements, limiting the response to intra-band
transitions. For the temperature scan we use a conduction band hopping
$t_{2}=-0.30$eV and a band gap of $\Delta_{0}=0.1$eV while for the chemical
potential scan we use $t_{2}=-0.25$eV and $\Delta_{0}=1$eV. On the calculated
temperature range, full momentum-grid convergence is achieved for $60\times
60\times 60$ $\mathbf{k}$-points. Contrary to the metallic system, the
transport coefficients of this insulator _must_ be evaluated with quadruple
precision for $\Gamma>0$ in order to avoid numeric instabilities in the
polygamma functions at low temperatures. Instead, evaluating the Boltzmann
kernels, a much denser momentum grid, that scales similarly to the metallic
system of Sec. 5.2 is required for convergence (not shown). Therefore, the
lifetime broadening in our more sophisticated treatment of scattering
amplitudes actually facilitates the numerical evaluation compared to semi-
classical approaches.

(a) Prototypical $\mu_{\psi/\mathrm{FD}}(T)$

(b) $\mu_{FD}(T)$ guidance via impurity states

Figure 10: Prototypical temperature dependencies in a semiconductor of the
resistivity $\rho$ (top), the coefficients of Seebeck $S$, Hall $R_{H}$, and
Nernst $\nu$ (second to fourth panel), as well as the chemical potential
(bottom; the band gap is delimited by thick black lines): (a) Effects stemming
purely from the occupation distribution (Fermi-Dirac ”FD” or digamma
”$\psi$”). (b) Effects achievable by guiding the chemical potential with
impurity states. Red: donor level $E^{D}=0.04$eV,
$\rho_{D}=10^{-8}\frac{1}{\mathrm{unit}\;\mathrm{cell}}$; Blue: acceptor level
$E^{A}=-0.02$eV, $\rho_{A}=10^{-14}\frac{1}{\mathrm{unit}\;\mathrm{cell}}$. In
all cases $\Delta_{0}=100$meV, $\Gamma=10^{-4}$eV; $Z=1$. The dashed line in
panel (a) shows the result of Boltzmann transport calculations.

####  5.3.1 Temperature scan

The introduced band asymmetry leads to an imbalance of electronic carriers as
the chemical potential lies above $\mu=0$, see bottom panel of Fig. 10a. The
thermal activation across the charge gap is reflected in all considered
physical observables where one naturally finds $\rho,R_{H}\propto
e^{-\frac{\Delta_{0}}{2kBT}}$ and $S,\nu\propto\frac{1}{T}$ at high enough
temperatures. At low enough temperatures, however, a saturation regime is
observed in $\rho$ and $R_{H}$, stemming from the Kubo transport kernels, see
Ref. [2] for details. Further, entropy transport is thermodynamically
consistent, as $S\propto T$ and $\nu\propto T$ towards absolute zero. In
Boltzmann theory, instead, the thermoelectric responses $S$ and $\nu$
unphysically diverge towards zero Kelvin (see dashed lines). Also $\rho$ and
$R_{H}$ diverge—in disagreement with experimental findings for narrow-gap
semiconductors[40].

Taking lifetime broadening into account also for the occupation, the chemical
potential $\mu_{\psi}(T)$ is necessarily (since $|t_{2}|>|t_{1}|$) forced
towards the conduction band from which a second activated regime is generated
below $100$K, see Fig. 10a (top panel). This _intrinsic_ chemical potential
behavior can also be achieved, in principle, through a guidance of
$\mu_{\mathrm{FD}}(T)$ via _explicit_ impurity states, see Fig. 10b. Here,
sharper transitions between regimes and even changes in the dominant type of
carriers can be generated. Please note the different plotting scales and the
sign changes of the Seebeck and Hall coefficient. Also note that the impurity-
engineered chemical potential can account, at best, for the intermediate and
high temperature response. The saturation of the resistivity and the Hall
coefficient, as well as the vanishing of the thermoelectrical responses
crucially relies on LinReTraCe’s Kubo kernels (for an in-depth discussion see
Ref. [2]). Hence, we propose LinReTraCe as the method of choice for future
transport studies of narrow-gap semiconductors.

(a) $\sigma(\mu)$ and $S(\mu$) for various temperatures with
$\Gamma=10^{-4}$eV.

(b) $\sigma(\mu)$ and $S(\mu$) for various temperatures with
$\Gamma=10^{-2}$eV.

Figure 11: Chemical potential scans for the symmetric insulator
$t_{1}=-t_{2}=0.25$ eV for (a) small $\Gamma=10^{-4}$eV and (b) large
$\Gamma=10^{-2}$eV scattering rates. Blue-to-green shades represent
conductivities; red-to-orange shades Seebeck coefficients. In the Boltzmann
temperature regime ($\Gamma=10^{-4}$eV; $T>250$K) we find perfect agreement
with the Goldsmid rule $2eS_{\mathrm{max}}T=\Delta_{0}$. Lower temperatures
and increased scattering leads to a departure of the Boltzmann regime and the
transport responses qualitatively change shape. The white (grey shaded)
background indicates the gap (valence/conduction band region). The gap is
$\Delta_{0}=1$eV and we used $Z=1$.

####  5.3.2 Chemical potential scan

The chemical potential scans in Fig. 11 illustrate the in-gap behavior of the
(intra-band) conductivity and the (intra-band) Seebeck coefficient. For small
scattering rates and elevated temperatures (Fig. 11a), we recover the usual
‘S’-shaped curve for the Seebeck coefficient: $S$ crosses zero at $\mu=0$,
where the system is particle-hole symmetric. Slightly above (below) a minimum
(maximum) develops, corresponding to the dominant type of carriers ($\mu>0$:
electrons; $\mu<0$: holes). In the ‘Boltzmann’ regime we find perfect
agreement with the Goldsmid rule $2eS_{\mathrm{max}}T=\Delta_{0}$ [82],
relating the maximal Seebeck coefficient $S_{\mathrm{max}}$ to the system’s
gap. At lower temperatures, however, deviations from this behavior can be
observed as a plateau around $\mu=0$ develops that expands as we lower the
temperature further. This effect stems from the Kubo kernels Eqs. (23-28) and
signals the transition from the activated to a saturation regime. For a
detailed discussion, see Ref. [2]. Larger scattering rates (Fig. 11b) already
lead to a deviation from the ‘Boltzmann’ behavior at the highest temperatures,
where we additionally find a strong suppression towards the band edges at
$\mu=\pm 0.5$eV.

###  5.4 Honeycomb lattice (tight-binding calculation)

To showcase the post-processing capabilities of LinReTraCe we revisit the
honeycomb structure of Sec. 4.1: We define the lattice vectors of length $1$Å
via

$a_{1}=\left(\vbox{\halign{\hfil$#$\hfil\cr 1\cr
0\cr}}\right);\;a_{2}=\frac{1}{2}\left(\vbox{\halign{\hfil$#$\hfil\cr
1\cr\sqrt{3}\cr}}\right)$ (82)

resulting in the reciprocal lattice vectors

$b_{1}=\frac{2\pi}{\sqrt{3}}\left(\vbox{\halign{\hfil$#$\hfil\cr\sqrt{3}\cr-1\cr}}\right);\;b_{2}=\frac{4\pi}{\sqrt{3}}\left(\vbox{\halign{\hfil$#$\hfil\cr
0\cr 1\cr}}\right),$ (83)

fulfilling $a_{i}\cdot b_{j}=2\pi\delta_{i,j}$. The special points K and K’
(see Fig. 13) are then located at

$\displaystyle K$ $\displaystyle=\frac{2}{3}b_{1}+\frac{1}{3}b_{2}$ (84)
$\displaystyle K^{\prime}$ $\displaystyle=\frac{1}{3}b_{1}+\frac{2}{3}b_{2}$
(85)

hence momentum grids that are divisible by three in each direction are
required to include them. We choose $n_{k_{x}}\times n_{k_{y}}=300\times 300$
with nearest neighbor hopping $t_{AB}=t_{BA}=1$eV. The corresponding tight
binding file232323This honeycomb structure, among other structures, including
kagome, 3D bcc, 3D fcc, etc., is available in the code repository’s templates
subfolder. reads:

begin hopping 0 0 0 1 2 1.0 0 0 0 2 1 1.0 +1 0 0 2 1 1.0 0 +1 0 2 1 1.0 -1 0 0
1 2 1.0 0 -1 0 1 2 1.0 end hopping begin atoms 1 0.25 0.25 0 1 0.75 0.75 0 end
atoms begin real_lattice 1 0 0 0.5 0.8660254038 0 0 0 1 end real_lattice

The resulting band-structure and density of states are illustrated in Fig. 12.
The particle-hole symmetry around $\varepsilon=0$ results in symmetric
transport properties $\sigma_{\alpha\beta}(+\mu)=\sigma_{\alpha\beta}(-\mu)$.
Here, we consider $\mu=0.8$eV, as indicated by the horizontal dashed line.
Fig. 13 (left) shows the momentum-resolved optical elements
$M^{\alpha\beta}\propto v_{\alpha}v_{\beta}$ in and around the Brillouin zone.
Combined with the kernel function $\mathcal{K}_{11}$ (middle), we can capture
how each point in the Brillouin zone contributes to the conductivity
$\sigma_{\alpha\beta}=\sum_{\mathbf{k}\in\mathrm{BZ}}\sigma_{\alpha\beta}(\mathbf{k})$
(right panel). Via an interplay of symmetry, the sign, and the values of the
optical elements we find the expected result: $\sigma_{xx}=\sigma_{yy}$ and
$\sigma_{xy}\equiv 0$. This type of analysis can be easily extended to other
transport coefficients and models and can provide valuable microscopic
insight.

Figure 12: Band-structure of the honeycomb lattice (left) and the density of
states (right). The analysis of Fig. 13 is done for $\mu=0.8$eV (dashed line).
Figure 13: Momentum analysis of the honeycomb lattice. Left: Optical elements
in units of $[\mathrm{eV}^{2}\mathrm{\AA}^{2}]$ for the $xx$ (top), $yy$
(middle) and $xy$ (bottom) polarizations (when neglecting the second term in
Eq. (57)). Middle: Kernel function $\mathcal{K}_{11}$, reflecting the Fermi
surface at $\mu=0.8$eV. Right: Momentum-resolved conductivity
$\sigma(\mathbf{k}$) in units of $[\frac{1}{\Omega}]$ which is the left and
middle column multiplied with each other.
$\sum_{\mathbf{k}}\sigma(\mathbf{k})$ results in $\sigma_{xx}=\sigma_{yy}$ and
$\sigma_{xy}\equiv 0$. Figure 14: Hybridization-gap semiconductor FeAs2.
Result of a WIEN2k calculation using the PBE exchange-correlation potential
and $20\times 20\times 40$ k-points. (top) total and partial density of states
(DOS); the DOS around the gap ($\Delta=0.275$eV) is dominated by Fe-character.
(bottom) LinReTraCe intra- and inter-band conductivities in the $x$-direction
and as a function of the chemical potential $\mu$. The Peierls approximation
to the optical elements leads to intra-band results (blue) virtually
indistinguishable to using the full dipole matrix elements (yellow): For the
given temperature ($T=300$K) and scattering rate ($\Gamma=10^{-5}$eV) the
intra-band contributions largely dominate over inter-band effects (green).

###  5.5 FeAs2 (WIEN2k calculation)

As a realistic crystal structure we first consider FeAs2, a hybridization-gap
semiconductor. Band-theory yields a gap of around $\Delta_{0}=0.275$eV when
using the PBE exchange-correlation potential[83, 4]. The top panel of Fig. 14
shows the ensuing density of states (DOS) for a wide energy range, with the
gap centered at $\mu=0$. The states near the gap edges are dominantly of iron
character. The bottom panel shows the resulting conductivities $\sigma$
calculated for $T=300$K, using a scattering rate $\Gamma_{0}=10^{-5}$eV.
Comparing the intra-band conductivity using the WIEN2k dipole matrix elements
(yellow) with that employing group velocities constructed from the BoltzTraP2
band interpolation (blue) illustrates the accuracy of the Peierls
approximation for simple unit cells.242424We cross-checked our result with
BoltzTraP2: Differences are small. We note that BoltzTraP2 calculates
transport distribution functions that are based on an artificially broadened
density of states, leading to a numerical smearing of the gap edges. At its
core, LinReTraCe is designed to stay numerically exact and avoids unphysical
broadening of any transport quantity. If desired, however, the user can apply
a Gaussian broadening of a chemical potential scan in the post-processing. For
explicit comparisons, the BoltzTraP2 conductivities $\sigma/\tau_{0}$
$[\frac{1}{\Omega\mathrm{m}\mathrm{s}}]$ have to be multiplied with
$\tau_{0}=\frac{\hbar}{2\Gamma_{0}}$, $\hbar=6.58211956\cdot 10^{-16}$eVs.
With the full dipole matrix element, we have access also to inter-band
contributions to the conductivity (green). As expected, inter-band
contributions only play a subsidiary role at these elevated temperatures: They
are two orders of magnitude smaller than the intra-band contributions.

Figure 15: FeAs2: WIEN2k intra-band conductivity ($\sigma$), Seebeck
coefficient ($S$) and power factor ($S^{2}\sigma$) for the same temperature
($T=300$K) and scattering rate ($\Gamma=10^{-5}$eV) as in Fig. 14. Here, we
show the dependency on the electron doping
$n=\frac{N(\mu)-N_{\mathrm{neutral}}}{V}$, instead of the chemical potential
$\mu$. The three Cartesian directions ($x$, $y$, $z$) refer to the
orthorhombic unit-cell with volume $V=a_{x}\times a_{y}\times
a_{z}=5.30\mathrm{\AA}\times 5.98\mathrm{\AA}\times 2.88\mathrm{\AA}$. The
gray, shaded area indicates doping levels where the chemical potential has
moved inside the conduction band.

Translating the $\mu$-axis into carrier concentrations $n={(N(\mu)-N)}/{V}$
results in the behavior shown in Fig. 15. For the given temperature, the
shaded gray area marks the region, where the chemical potential has moved
inside the conduction band. Upon entering said region, the conductivity
increases rapidly, while the Seebeck coefficient has its peak amplitude for
$\mu$ inside the gap, cf. Figure 11. Indeed, the behavior of $\sigma$ and $S$
is typically[84], but not always[85], antagonistic. The response is
polarization dependent, with the powerfactor $S^{2}\sigma$ being most
sensitive on the crystal orientation.252525Please note that these result use
the full dipole matrix element, for which WIEN2k has a maximal number of
computable $\mathbf{k}$-points. Results for the metallic regime could
therefore not be checked for momentum-convergence. This limitation does not
apply for Peierls velocities.,262626LinReTraCe only computes thermoelectric
properties from pure electron diffusion. Phonon-drag enhancements—relevant,
e.g., for the related narrow-gap semiconductors FeSb2[86, 87, 88] and
CrSb2[89, 90]—are not included.

Finally, we add an in-gap impurity state and showcase the ensuing temperature
behavior of transport observables in comparison to experiment[63, 91] in Fig.
16. Due to the additional electrons provided by the donor level, located in
the vicinity of the conduction band, the system experiences a transition from
activated behavior with the intrinsic gap $\Delta_{0}$ at high temperatures to
a second regime controlled by a reduced gap $\Delta_{1}$, determined by the
position of the impurity, cf. Fig. 4 and Ref. [2]. At still lower
temperatures, and akin to Sec. 5.3, the prototypical saturation regimes set
in: The resistivity and the Hall coefficient saturate, while the Seebeck and
Nernst coefficients tend towards zero in a linear fashion. Due to the three
distinct crystal directions in FeAs2, the Hall and Nernst coefficients behave
differently depending on the direction of the applied magnetic field.

Figure 16: FeAs2: Intra-band resistivity $\rho$, Seebeck coefficient $S$, Hall
coefficient $R_{H}$ and Nernst coefficient $\nu$ (left) compared to
experimental data (right) [63]. We show results for all three Cartesian
directions. The chemical potential was determined by the Fermi Dirac function
and altered by a donor level (at $E_{D}=0.015$eV below the conduction band
with $\rho=2\cdot 10^{-6}\frac{1}{\mathrm{unit}\;\mathrm{cell}}$), chosen so
that the resistivity approximates the experimental behavior. The scattering
rate includes only polynomial terms: $\Gamma(T)=7\cdot 10^{-5}+5\cdot
10^{-9}T^{2}$ in eV. Note that we underestimate $S$ and overestimate $\nu$.
All shown data use the intra-band Peierls approximation.

###  5.6 Tl-doped PbTe (WIEN2k calculation)

As a final test material we consider Tl-doped PbTe, a prime example for the
enhancement of thermoelectric transport by resonance states [92]. We model the
doping by explicitly constructing a $4\times 4\times 4$ supercell of PbTe and
substituting a single Pb atom with a Tl one. Internal positions of the
resulting Tl0.004Pb0.996Te are then fully relaxed before extracting the
LinReTraCe input. In order to gain access to both the Seebeck and the Nernst
coefficient, we employ the BoltzTraP2 interpolation scheme to obtain Peierls
velocities. As the prepared super cell does not match the experimental
stoichiometry, we are mainly aiming for qualitative aspects.

Through the Tl-(hole)doping an increase of the density of states in the
vicinity of the valence-band edge appears, see Figure 17 and the emerging
resonance pins the chemical potential.

Figure 17: Tl-doped PbTe: Density of states using the LDA exchange-correlation
potential and $17\times 17\times 17$ $\mathbf{k}$-points. Introducing a single
Tl$\Longleftrightarrow$Pb substitutional impurity results in a distorted
electronic structure and additional states appear near the valence-band edge,
pinning the chemical potential. The zoomed inset illustrates the non-trivial
nature of the resonance to which all three elements contribute in some part.

Consequently, the Tl-doping cause a semiconductor-to-metal transition.
Supplementing the electronic structure with a phenomenological, band-
independent $\Gamma(T)=20$meV $+6\cdot 10^{-5}$meV/K${}^{2}T^{2}$ results in
the resistivity, Seebeck, and Nernst coefficient shown in Fig. 18. Comparing
to experimental results[92], plotted on the same scale, we find good agreement
for the resistivity at elevated temperature. However, we do not capture the
resonance’s signature below 200K. This finding suggests that the resistivity
in that temperature range is not controlled by band-structure effects.
Instead, the scattering rate might be markedly different for the resonance
than for other states and display a more complex temperature dependence. We
hence expect that—with scattering rates obtained from more sophisticated
(beyond band-theory) electronic structure methodologies—LinReTraCe could
capture the resistivity also below 200K. Alternatively, a reverse engineering
approach could be employed within LinReTraCe to extract a phenomenological,
possibly band-dependent, scattering rate that reproduces the experimental
resistivity.

The Seebeck and Nernst coefficient, instead, are well reproduced within our
setting, without any additional input, suggesting a simple electronic picture
of thermoelectric transport to hold. Indeed, in the case of metals, the
Seebeck coefficient is—to a first approximation—insensitive to the scattering
rate. The congruence to experiment for $S$ then supports the above claim that
the temperature profile of the resistivity below 200K is controlled by an
intricate scattering rate.

The metallic nature of the transport poses the previously discussed challenges
for the Brillouin-zone discretization. While the resistivity and the Seebeck
coefficient are, for all practical purposes, convergent with the largest
$\mathbf{k}$-mesh used, the Nernst coefficient is notably more sensitive: The
shown data provides a good approximation for the Nernst coefficient above
200K, while below the result is clearly not yet converged. Indeed, for general
reasons, the Nernst coefficient must vanish for $T\rightarrow 0$[2]. We stress
that these limitations regarding the number of usable $\mathbf{k}$-points are
on the side of the electronic structure methodology, while LinReTraCe could
handle larger meshes at acceptable costs.

Figure 18: Tl-doped PbTe: Resistivity (top) and coefficients of Seebeck
(middle) and Nernst (bottom) within LinReTraCe (left) compared to
experiment[92] (right). Even though the resonance feature of the resistivity
below 200K cannot be reproduced with our phenomenological scattering rate, the
Seebeck and Nernst coefficients are matched nearly one-to-one. Indeed, while
the resistivity is directly controlled by our ansatz,
$\Gamma(T)=20$meV$+6\cdot 10^{-5}$meV/K${}^{2}T^{2}$, for the scattering, the
Seebeck coefficient in metals is (to first approximation) independent of
$\Gamma$. There are opposing conventions for the sign of the Nernst
coefficient[62]. Since our Hall signal (not shown) has the same sign as the
experiment[92], but not the Nernst coefficient, we multiply the experimental
$\nu$ with $(-1)$. Please also note that $\nu\rightarrow 0$ for $T\rightarrow
0$ is a must. Indeed, due to the limited number of $\mathbf{k}$-points, the
Nernst simulation is not fully converged below $T=200$K.

##  6 Conclusion

We presented the algorithm, implementation, and usage of LinReTraCe, a package
for the computation of transport properties of solids. The code’s strengths
are (1) based on Kubo linear-response theory it captures effects of finite
lifetimes far beyond semi-classical theories, making a qualitative difference
in semi-conductors; (2) using a semi-analytical evaluation of transport
kernels, LinReTraCe is nonetheless as fast as Boltzmann approaches in the
relaxation-time approximation; (3) LinReTraCe is agnostic to the origin of the
input data, allowing to simulate tight-binding models as well as materials.
Interfaces to WIEN2k, VASP, and Wannier90 are included in the release, while a
template allows connecting LinReTraCe to any other electronic structure code.
The numerical efficiency makes LinReTraCe also attractive for high-throughput
material surveys of, e.g., thermoelectrics.

## Acknowledgements

The authors acknowledge discussions with A. Georges, J. Mitscherling, and J.
Mravlje. Calculations were performed on the Vienna Scientific Cluster (VSC).

##### Author contributions

M.P. and E.M. developed the methodology and the program package. M.P. and
J.M.T. created the manuscript. J.M.T. designed and supervised the research.

##### Funding information

LinReTraCe was funded by the Austrian Science Fund (FWF) through project P
30213.

##  Appendix A Cheat sheet

Code repository: https://github.com/linretrace/linretrace.
Example repository and tutorial: github.com/linretrace/linretrace_examples
https://github.com/linretrace/linretrace/raw/release/documentation/tutorial.pdf.
LinReTraCe workflow:
Create an energy file from WIEN2k, VASP, WANNIER90, a tight binding file or
generic data

$ ldft <wien2k folder> \--optic # wien2k with full dipole elements $ ldft
<vasp folder> \--interp # vasp with intra-band interpolation $ lwann
<wannier90 folder> # wannier90 $ ltb <tb file> <nkx> <nky> <nkz> <charge> #
tight-binding input $ linterface # your customized interface

The full capability of the Python interfaces can be shown via:

$ ldft/lwann/ltb --help

The created energy files can be inspected with `lprint`:

$ lprint <lrtc input> info $ lprint -p <lrtc input> dos $ lprint -p <lrtc
input> path

_Optional:_ Prepare a scattering file by copying `lscat_temp` into your
working directory and adjust it to your needs.

Next, prepare a basic config file with

$ lconfig

Advanced option can be added with the help of the full configuration
specifications found in `documentation/configspec`.

Then, perform LinReTraCe calculation via:

$ mpirun -np <cores> bin/linretrace config.lrtc

When finished, the datasets contained in the output can be inspected via:

$ lprint <lrtc output> list # physical datasets $ lprint <lrtc output> olist #
Onsager coefficients

Printing or plotting (-p option) the data:

$ lprint -p latest mu --gap # chemical potential + gap $ lprint -p latest
c-intra uxx dxx # conductivity $ lprint -p latest r-intra xx # resistivity $
lprint -p latest p-intra xy # Peltier $ lprint -p latest s-intra xx yy zz #
Seebeck $ lprint -p latest rh-intra xyz # Hall $ lprint -p latest n-intra xyz
# Nernst $ lprint -p latest L12B-intra uxyz # L12B Onsager

##  Appendix B Advanced LinReTraCe configuration

Section 3.6 provides a minimal configuration for temperature and chemical
potential runs. Here, we detail more advanced capabilities of LinReTraCe. A
complete documentation of the config specifications can be found in
`documentation/configspec`.

###  B.1 General settings

In the `[General]` section of the config file

[General] Bandgap = <value> ElectronOccupation = <value> FermiOccupation = T/F
# default = T FullOutput = <string> # <string> = full -> momentum resolved,
band resolved # <string> = ksum -> momentum summed, band resolved # <string> =
bsum -> momentum resolved, band summed

one can manipulate the band gap, change the total occupation, and activate the
output to be fully or partially resolved into momentum and band contributions.
The band gap can only be changed if a gap has been detected in the pre-
processing. Hence a properly set occupation in these setups is necessary. The
change in the band gap is achieved via an upwards shift of all the conduction
bands, while the valence bands remain fixed. For spin-polarized systems, two
values need to be provided, one for each spin.

The change in total occupation via the config file triggers a re-search of the
‘DFT’ chemical potential (having a direct impact on chemical-potential run,
that are defined with respect to the updated $\mu_{\mathrm{DFT}}$).

The used occupation function can be switched between the digamma function
$f(\varepsilon)=\frac{1}{\pi}-\Im\psi(\frac{1}{2}+\frac{\beta}{2\pi}(\Gamma+ia))$
and the Fermi function $f(\varepsilon)=[1+e^{\beta a}]^{-1}$ via the
`FermiOccupation` key word. Please note that this affects the chemical
potential as well as all quantities depending on it, including the carrier
concentration, total energy, etc.

The `FullOutput` option can be configured to either output the fully resolved
(`full`), or partially resolved (momentum sum: `ksum`, band sum: `bsum`)
Onsager coefficients. The tree structure of the full dependence output is (in
contrast to the standard output) separated for each step, otherwise only the
dataset shapes and their identifier differ

/step/000001 Group /step/000001/L11 Group /step/000001/L11/intra Group
/step/000001/L11/intra/full Dataset # FullOutput = full
/step/000001/L11/intra/ksum Dataset # or FullOutput = ksum
/step/000001/L11/intra/bsum Dataset # or FullOutput = bsum ...

where the datasets without magnetic field have the array shape
`[nkp,spins,bands,3,3]` and the datasets with magnetic field are of shape
`[nkp,spins,bands,3,3,3]`. Partially momentum (band) summed datasets have `nkp
= 1` (`bands = 1`). The steps are always ordered the same way as the summed
quantities, that is from lowest to highest temperature or from lowest to
highest chemical potential, independent on how the program itself loops over
them. These datasets can only be accessed via external HDF5 libraries, an
exemplary load in via `h5py` looks like

import h5py with h5py.File(’lrtc-output.hdf5’,’r’) as h5: L11full =
h5[’step/000001/L11/intra/full’][()] L11Bfull =
h5[’step/000001/L11B/intra/full’][()] print(L11full.shape) # kpoints, spins,
bands, 3, 3 print(L11Bfull.shape) # kpoints, spins, bands, 3, 3, 3

###  B.2 Temperature mode settings

The temperature sub configuration includes the following options:

[TempMode] ConstantMu = <value> OldMuHdf5 = <LRTCoutput file> OldMuText =
<lprint mu output> NImp = <value> Doping = <value> NominalDoping = T/F #
default = F

`ConstantMu` sets the chemical potential to a fixed absolute value (in eV)
with respect to $\mu_{\mathrm{DFT}}$ for the full temperature run. `OldMuHdf5`
uses the chemical potential values of an old output file. `OldMuText` can be
used to enter manual values via a text file whose format has to be consistent
with the output of `lprint <LRTCoutput file> mu` (column 1 temperature, column
2: chemical potential). Naturally these three option cannot be used in
combination.

Impurity levels and/or homogeneous doping, as illustrated in Sec. 4.2.1 and
Sec. 4.2.2, respectively, can be activated via `Doping` and `NImp`. The value
in `Doping` is in units of $[\mathrm{cm}^{-3}]$ which can be changed to
nominal doping $[1/\mathrm{unit}\;\mathrm{cell}]$ by setting `NominalDoping =
T`. This behavior then also applies to the introduced $n$ impurity states,
activated via `NImp = n`, for which further configuration is required:

[[Impurities]] [[[1]]] # one of Absolute = <values> Valence = <values>
Conduction = <values> Percentage = <values> Bandtype = <string> # <string> =
box, triangle, halfcircle, sine, sine2, sine3, sine4 Bandwidth = <value>
[[[2]]] ...

Separate groups `[[[1]]] ...` have to be introduced to describe each of the
$n$ impurities. The main descriptor is _one_ of the four keywords `Absolute,
Valence, Conduction` or `Percentage`, as they describe how the provided energy
is to be interpreted. The provided string then contains 4 values, that
describe the type of impurity (donor: $+1$/acceptor $-1$), its density,
energy, and degeneracy. The energy is translated into an internal position,
depending on the keyword: `Absolute` uses the energy is as, `Valence`
interprets the energy as positive offset with respect to the top of the
valence band, `Conduction` as a negative offset with respect to the bottom of
the conduction band. `Percentage` puts the impurity state at the fractional
position inside the gap ($0<E_{\mathrm{Percentage}}<1$). Naturally the latter
three options can only be used in gapped systems. Energies refer to bare
values, prior to potential renormalizations with $Z$ and $\Re\Sigma$. Example
inputs for impurity states look as follows

# Identifier = type density energy degeneracy Absolute = -1 1e15 8.35 2 #
acceptor at E = 8.35 eV Conduction = +1 1e18 0.02 1 # donor 20meV below CB
Valence = -1 1e19 0.01 1 # acceptor 10meV above VB Percentage = +1 1e15 0.5 2
# donor at center of gap

These states can be transformed into bands by providing the `Bandtype`
identifier (cf. Section 4.2.1) and an associated `Bandwidth`. The introduced
bandwidth is centered around the provided energy and has the form of the
provided `Bandtype` string. All these impurity states only affect the chemical
potential; they do not contribute to the transport.

###  B.3 Scattering files

Runs with a provided scattering file are configured with

[TempMode] # or [MuMode] [[Scattering]] ScatteringFile = <scattering file>
ScatteringOffset = <value>

as the scattering file already contains the defining calculation axis. Here we
allow for an additional offset value that is added on top of the provided
scattering dependencies, i.e. $\Gamma(k,n,\mu|T)+\Gamma_{\mathrm{offset}}$,
e.g., to mimic additional impurity scattering. This subsection can be used in
the temperature, as well as the chemical potential sub configuration. However
we cross-check the internally set calculation mode of the scattering file with
the calculation mode of the config file.