Interband tunneling effects on materials transport properties using the ﬁ rst principles Wigner distribution

Electronic transport in narrow gap semiconductors is characterized by spontaneous vertical transitions between carriers in the valence and conduction bands, a phenomenon also known as Zener tunneling. However, this effect is not captured by existing models based on the Boltzmann transport equation. In this work, we propose a new fully ﬁ rst principles model for electronic transport using the Wigner distribution function and implement it to solve the equations of motion for electrons. The formalism generalizes the Boltzmann equation to materials with strong interband coupling and include transport contributions from off-diagonal components of the charge current operator. We illustrate the method with a study of Bi 2 Se 3 , showing that interband tunneling dominates the electron transport dynamics at experimentally relevant small doping concentrations, a behavior that is likely shared with other semi- conductors, including topological insulators. Surprisingly, Zener tunneling occurs also between band subvalleys separated by energy much larger than the band gap. © the CC BY (http://creativecommons.org/licenses/by/4.0/). We apply ﬁ study topological We show that


Introduction
Small band gap semiconductors are expected to possess transport characteristics that differ significantly from conventional wide-gap materials due to Zener (or Klein) tunneling, a phenomenon in which carriers undergo vertical transitions across the band gap, which can substantially increase electrical conductivity. This phenomenon is relevant for a number of novel systems, such as small band gap nanotubes or multi-layer graphene-based systems [1,2]. Topological insulators are another important example, thanks to their variety of interesting physical properties and promising applications, such as low-power electronics and quantum computing [3e7]. Topological insulators are characterized by conductive surfaces, while their bulk phases have a small band gap, and they may therefore display Zener tunneling and enhanced bulk conduction, which warrants investigation.
The de-facto tool of choice for first-principles simulation studies of electronic transport properties in crystals is the ab initio Boltzmann Transport Equation (aiBTE), which provides estimates of transport properties in good agreement with experimental measurements (e.g. [8e13]). However, this semiclassical approach is not always sufficient to model electronic transport properties. Zener tunneling, for example, is not captured by semiclassical models, as they omit contributions from off-diagonal components of the flux operators. Sophisticated models based on the nonequilibrium Green's function methods [14] are capable of overcoming the limitations of semiclassical models. However, lack of efficient numerical and theoretical formulations prevents the combination of these approaches with first-principles computations that can treat realistic materials, and neglect the critical influence of scattering on the charge flux. Here we explore a new approach based on the Wigner function: this formalism is particularly appealing as it is exact in principle and, as discussed below, reduces to the aiBTE when off-diagonal flux components are neglected. The use of the Wigner function to study materials' electronic properties has been limited so far to studies of model Hamiltonians (see e.g. Ref. [15] for a recent review). In contrast, first-principles modeling can provide not only a more accurate description of real materials' Hamiltonian, but also an accurate description of electronic scattering that is responsible for the relaxation of the out-of-equilibrium system. This work formulates the generalized transport approach building on recent progress in first-principles computations of electron scattering rates and fills this methodological gap, opening possibility for the study of firstprinciples electronic properties of realistic materials.
In this work, we derive and implement from first-principles an equation of motion for electrons, termed the ab initio Wigner transport equation (aiWTE) based on the single-particle Wigner distribution function. The formalism is capable of describing the space-time evolution of electrons, including effects due to electronic scattering, and captures the off-diagonal contributions of flux operators to transport properties. This equation is explicitly solved for the set of electronic transport coefficients, in particular electrical conductivity and Seebeck coefficient, showing how they can be significantly different from their semiclassical counterpart due to the presence of off-diagonal components of the Wigner distribution function. We apply the formalism to a first-principles study of the bulk transport properties of the topological insulator Bi 2 Se 3 . We show that at small doping concentrations, the aiWTE estimates of bulk electronic transport properties deviate significantly from semiclassical estimates, due to the presence of Zener tunneling, i.e. vertical transitions that couple carriers of valence and conduction bands.

Theory
We start by considering the ground state Hamiltonian H 0 of a crystal, which we assume to be an independent-particle Hamiltonian with eigenvalues ε bs ðkÞ and Bloch states j kbs ðxÞ, where x is the position, k the wavevector, b the band index and s ¼ ± 1 the spin index. The ground-state is perturbed by a constant electric-field E which couples with the carriers' charge e and by the electronphonon interaction H elÀph , so that the total Hamiltonian is H ¼ H 0 þ ex,E þ H elÀph . To derive an equation of motion for such a system, we use the single-particle Wigner function f of the system [16], defined as the Wigner transform of the density matrix r as where t is the time. We build the Wigner function through a transformation of the density matrix r bb 0 ss 0 ðk; k 0 Þ in the reciprocal space representation. Such Wigner transform consists in a rotation of variables k; k 0 / kþk 0 2 ; k 0 À k combined with a Fourier transform on one variable. The Wigner function operates in a phase-space representation, which is especially useful to draw connections between quantum and classical mechanics.
The evolution of the Wigner function [17,18] is found through a Wigner transform of the equation of motion of the density matrix, and has been shown to be vf bb 0 ss 0 ðx; k; tÞ vt ¼ À fff ðx; k; tÞ; Hðx; kÞgg bb 0 ss 0 : where fff ; Hgg is the Moyal bracket (the quantum mechanical extension of the Poisson bracket) and the Moyal product + is defined as where Hðx; kÞ is the Wigner transform of the Hamiltonian, and the arrows indicate that the derivative operator is acting to the left/ right operators. We now simplify the Hamiltonian supposing that the electronphonon interaction is weak and evaluate the Moyal bracket for the single-particle part of the Hamiltonian; the electron-phonon interaction is added later as a perturbation. Since we are interested in macroscopic properties, we can make the assumption that only slow spatial variations of the Wigner function are relevant. Therefore, we expand the Moyal product in Taylor series to the lowest orders of Z and find an equation of motion which we term the ab initio Wigner Transport Equation (aiWTE), that is where f; g is an anticommutator, E bb 0 ss 0 ðkÞ ¼ d bb 0 d ss 0 ε bs ðkÞ is a tensor of electronic energies, d bb 0 ss 0 ðkÞ ¼ ð1 À d bb 0 ÞCkbsjexjkb 0 s 0 D is a tensor of electric dipoles between two Bloch states (typically used to describe optical excitations), and v bb 0 ss 0 ðkÞ is the velocity operator, defined from the commutator of the Hamiltonian and the position operator as v ¼ i Z ½H; r. The electron-phonon scattering operator vf bb 0 ss 0 ðx;k;tÞ vt j coll is added as a perturbation to the aiWTE and is built using scattering rates from the Fermi Golden rule [19e23]. We refer to the Supplementary Information for a more detailed derivation of the aiWTE.
The aiWTE needs to be solved to obtain an estimate of the single-particle Wigner distribution function. As a first comment, we note that the aiBTE is recovered as a limiting case of the aiWTE, when the off-diagonal terms bsb 0 and sss 0 are set to zero. This corresponds physically to a situation when different bands do not couple via Zener tunneling. This can happen, for example, when neither thermal excitation nor dipole interaction provide sufficient energy to allow for the vertical transition of one particle from one band to another. Therefore, the most interesting terms to discuss in the aiWTE are the off-diagonal terms, which introduce the possibility of additional electronic transitions, or couplings, between different electronic states at a given wavevector k. We further note that some of the off-diagonal terms of the aiWTE shown here are absent in other works based on either the density matrix or the Wigner function [24e29]. Additionally, we note that the electronic aiWTE is conceptually similar to a formalism developed for phonon transport in Ref. [19], although here we use a simplified derivation and include the effect of external forces (the electric field).
The aiWTE can be solved to estimate transport coefficients with a technique similar to the one used for the aiBTE. For steady-state transport in a bulk system, the Wigner distribution is stationary in time and independent from the particular position inside the bulk; therefore, we look for a solution in the form f bb 0 ss 0 ðkÞ. Next, we look for the linear response to an applied external electric field E or A. Cepellotti and B. Kozinsky Materials Today Physics 19 (2021) 100412 a temperature gradient VT and write f bb 0 ss 0 ðkÞ ¼ f bs ðkÞ þ f Ei bb 0 ss 0 ðkÞE i or f bb 0 ss 0 ðkÞ ¼ f bs ðkÞ þ f Ti bb 0 ss 0 ðkÞV i T, where i is the direction of the applied perturbations and f bs ðkÞ is the FermieDirac distribution function. We further adopt the relaxation time approximation, which simplifies the treatment of the scattering operator. After inserting these two proposed solutions in the aiWTE, and retaining only terms up to linear order in E or VT, we find two equations for the diagonal part of the aiWTE: ev i bbss ðkÞ vf bbss ðkÞ vε and v i bbss ðkÞ , vf bbss ðkÞ vT where G bs ðkÞ is the carrier's linewidth. These equations are equivalent to the standard aiBTE for describing electronic transport within the relaxation time approximation. However, the aiWTE also gives rise to two equations for the off-diagonal terms of the Wigner distribution, which are and iðε bs ðkÞ À ε b 0 s 0 ðkÞÞf Ti bb 0 ss 0 ðkÞ þ describing the response of the system to electrical and thermal perturbations, respectively. Now the aiWTE is in a form that can be readily solved for f Ei and f Ti with some algebra, allowing us to reconstruct the Wigner distribution of a crystal in presence of an electric field or a thermal gradient.
Having found the Wigner distribution, transport properties are readily obtained. For example, the charge current density is where V is the crystal unit cell volume and N k the number of wavevectors used to integrate the Brillouin zone. Additionally, the heat flux is ðε bs ðkÞ À mÞfvðkÞ; f ðkÞg bbss : These definitions readily allow us to compute transport coefficients using an approach typical of transport theory. After inserting the solution to the aiWTE in the definition of j and q, it is readily seen that charge current and heat flux can be written in the form j ¼ L EE E þ L ET VT and q ¼ L TE E þ L TT VT, where L denotes the Onsager coefficients. Transport coefficients can be expressed in terms of these Onsager coefficients: the electrical conductivity s ¼ EE L ET , and the thermal conductivity k ¼ À ðL TT À L TE L À1 EE L ET Þ. We now inspect the electrical conductivity in more details. Since j is linear in f bb 0 ss 0 ðkÞ, the total electrical conductivity s is a sum of diagonal and off-diagonal contributions s ij ¼ s aiBTE ij þ Ds ij , where i; j are cartesian labels. Here, s aiBTE ij comes from the diagonal terms of the Wigner distribution (b ¼ b 0 and s ¼ s 0 ) and is equal to the estimate from the aiBTE: where t bs ðkÞ ¼ Z=G bs ðkÞ is the carriers' lifetime. This semiclassical conductivity is corrected by an additional term Ds ij that can be found inserting f Ei bb 0 ss 0 ðkÞ in the definition of j and is We can now understand the effects of the off-diagonal corrections. First, the correction Ds is positive (note that f bs ðkÞ is a monotonic function of ε bs ðkÞ ), and therefore the aiWTE will always adjust the semiclassical prediction of conductivity to higher values. Second, the correction depends on a few quantities: the energy difference between electrons and holes, their linewidths, and velocity. One scenario where this correction is relevant, for example, is whenever the energy difference between an electron and a hole is comparable to their linewidth, so that the two carriers interact. The strength of such interaction is determined by the velocity matrix element v bb 0 ss 0 ðkÞ, i.e. the matrix element for the optical transition. In addition, a strong dipole coupling can also occur between states that are far from the conduction or valence band edge, as discussed below. Whenever this off-diagonal correction is large, the system allows for an additional transport mechanism, known as Zener tunneling, in which electrons propagate by tunneling across the band gap. In contrast, we stress that Zener tunneling, or similar spontaneous vertical electronic transitions, are entirely omitted in the aiBTE formalism, since the aiBTE doesn't include off-diagonal contributions from the velocity operator.

Computational methods
All quantities appearing in the aiWTE are available from firstprinciples codes and we can therefore apply this formalism using fully ab-initio parameters.
We use density functional theory as implemented in the planewave software suite Quantum-ESPRESSO [30,31]. To compute the ground state, we use ultrasoft pseudopotentials from the GBRV library [32], with the PBEsol functional. We use an energy cutoff of 80 Ry, and integrate the Brillouin zone with a 8 Â 8 Â 8 mesh of kpoints. We build the trigonal unit cell using experimental estimates of the crystal structure [33], i.e. with a lattice parameter of 9.839 Å, and an angle a such that cosa ¼ 0:91068. The Wannier functions are computed using p orbitals on both Bi and Se atoms as initial guesses.
Phonon properties, and the electron-phonon matrix elements are computed with density-functional perturbation theory [34] on a coarse grid of 4 Â 4 Â 4 q-points. Electron-phonon matrix elements are subsequently interpolated on a finer grid of 39 Â 39 Â 39 Brillouin zone wavevectors using a mixed Wannier and linear interpolation [12], while electronic energies and velocities are interpolated using Wannier90 [35]. Transport properties have been implemented in a custom-made software. The Dirac-delta ensuring energy conservation during an electron-phonon scattering event has been approximated using an adaptive-smearing scheme [36]. Transport properties have been converged with respect to the kpoints mesh used to integrate the Brillouin zone. The scattering operator, detailed in the Supplementary Information, is built considering only intrinsic electron-phonon scattering. The electron-phonon coupling includes the effect of long-range polar interaction, which contributes significantly to transport coefficients, and is evaluated according to the methodology of Ref. [37]. While the scattering operator may be treated beyond the relaxation time approximation with techniques similar to those of Refs. [8,9,12], all results of this study are obtained within the relaxation time approximation, as it allows for significantly faster simulations. We note that results beyond the relaxation time approximation would only affect the diagonal components of the Wigner distribution, but not the off-diagonal components since the action of the scattering operator on the off-diagonal components of the Wigner distribution involve only the carriers' lifetimes/linewidths (see Eq. 18 of Supplementary Information).

Results
We now apply the formalism to study the intrinsic phononlimited electronic transport of bulk Bi 2 Se 3 . In Fig. 1 we report the band structure [38,39] and the density of states (DOS) for this narrow-gap semiconductor. We estimate a quasiparticle gap of 0.2 eV, in agreement with experimental estimates [40]. We also mention that, while DFT may not always be accurate in estimating the band gap, the qualitative aspects of the results discussed below do not strongly depend on the precise band gap value. It is worth noting that the DOS increases away from the Fermi level (set at 0 eV at the middle of the band gap) and flattens at energies of approximately À0.8eV and 1.0eV for the valence and conduction bands, respectively, indicating that the subvalleys are separated by an energy of approximately 1.8 eV.
In Fig. 2a ( band are excited, the average carriers' group velocity is small, due to the quadratic nature of the band minimum. Therefore, the semiclassical contribution to electrical conductivity tends to be rather small. The aiWTE corrects this picture, including the Zener tunneling effect [2]. As carriers from valence and conduction band are close in energy, they can interact and contribute to the electrical transport through the tunneling effect, as discussed above. For small dopings, the aiWTE correction is significant, and can be much larger than the aiBTE conductivity value. For the smallest value of doping reported (10 16 cm À3 ), this correction is largest at lower temperatures. The doping of 10 18 cm À3 is an intermediate case, with metallic behavior at low temperatures (and thus smaller aiWTE correction) and semiconducting (with larger aiWTE correction) at higher temperatures as the chemical potential moves from the conduction band into the band gap. We can thus conclude that a substantial portion of electrical current at low doping is carried through the Zener tunneling included in the aiWTE formalism: the current is not only caused by the carriers traveling at a finite group velocity, but also by carriers' tunneling between single-particle Bloch states. In Fig. 2, panels c and d, we report the Seebeck coefficient S for n-type (p-type) doping concentrations, with aiBTE results in dashed lines and aiWTE in solid lines. The Seebeck coefficient is usually expected to have negative values for n-type doping and positive for p-type. However, there are deviations from this behavior, and the Seebeck coefficient can even change sign when temperature is varied at fixed doping concentrations. This complex behavior is not caused by the off-diagonal terms of the Wigner distributions, and the Seebeck coefficient sign changes are observed when using the aiBTE as well, e.g. due to the bipolar effect where minority carriers are thermally excited across the band gap. Such behavior has been found for Bi 2 Te 3 [41] and CoSb 3 [42], where the Seebeck coefficient exhibits marked temperature dependence.
The Wigner correction tends to make the Seebeck coefficient smaller in absolute value, which may be attributed to the electrical conductivity appearing at the denominator of the definition of the Seebeck coefficient. As for the case of electrical conductivity, the off-diagonal correction is more prominent at small values of doping concentrations and low temperatures. This may have a significant consequence when estimating thermoelectric properties, such as the power factor that depends on the square of the Seebeck coefficient. In this case, the inclusion of off-diagonal components of the Wigner distribution lead to a revision of power factor estimates to smaller values.
We now examine the relationship of electrical conductivity s and open-circuit electronic thermal conductivity k el . The Wiedemann-Franz (WF) law defines the Lorenz number L ¼ kel sT , which in the ideal metallic limit is a constant L 0 ¼ 2:44,10 À8 WUK À2 . Knowledge of L is necessary to decouple the electronic contribution k el and the lattice contribution k ph from measurements of the total k. In Fig. 2 e and f we plot the computed ratio L L0 Figure. 1. Electronic band structure and density of states of Bi2Se3. Electronic band structure and density of states (DOS) of Bi2Se3 near the Fermi energy, set at the center of the band gap. The bulk crystal is characterized by a small gap, opened by the spinorbit coupling. We also note from the DOS that subvalleys of valence and conduction bands are approximately 1.8 eV apart in energy.

A. Cepellotti and B. Kozinsky Materials Today Physics 19 (2021) 100412
for several temperatures and n-type or p-type kind of doping concentrations, using both aiBTE and aiWTE in dashed and solid lines respectively. At high doping the system has metallic character, and both predictions closely follow the WF law. In the case of small doping, in the bipolar regime, the semiclassical aiBTE predicts large deviations from the WF law, as has been discussed previously [13, 43e45]. Remarkably, in the aiWTE solution the Lorenz numbers are much closer to the expected range for semiconductors, indicating that quantum transport effects included in the aiWTE strongly suppress deviations and work towards restoring the validity of the WF law. In Fig. 3a, we analyze the contributions to the aiBTE electrical conductivity as a function of the carriers' energy at doping concentration of 10 18 cm À3 , and temperature of 700 K. This histogram is built such that the area under the curve integrates to the total electrical conductivity. Within the semiclassical relaxation time approximation, the quantity plotted is an energy-resolved histogram of e 2 V Nk vf bs ðkÞ vε v 2 bbss ðkÞt bs ðkÞ, i.e. the contribution of a single mode to the aiBTE electrical conductivity. As expected, the dominant contributions to electrical conductivity come from carriers whose energy is close to the chemical potential (set at 0 eV). The contributions of other carriers decay exponentially as their energy gets further from the chemical potential.
The aiWTE correction Ds cannot be resolved in terms of a single carrier's energy, since it involves the tunneling between two states at different energies. Therefore, in Fig. 3b, we plot the contributions to the electrical conductivity as a function of two interacting carriers energies. On the diagonal, we find again the aiBTE-like terms shown in Fig. 3b. In addition, we can see important off-diagonal contributions to the electrical conductivity that are not present in the aiBTE and are introduced with the aiWTE. In particular, there are two peaks of contributions to electrical conductivity, that couple electrons of energy 1.0 eV with holes at À0.8 eV. These two values correspond to the average energies of the valence and conduction bands, when the DOS reaches the corresponding maximum values. Therefore, in contrast to the typical intuition of the Zener tunneling, we find that the most significant coupling between carriers takes place far from the chemical potential, with carriers of energy much larger than the thermal energy. For this material, the dipole interaction between carriers in subvalleys of the valence and conduction is thus particularly strong, allowing for high-energy carriers to contribute to transport. As a result, we speculate that Zener tunneling may take place also in semiconductors with a wide gap and contribute significantly to electrical conductivity, provided that the inter-band dipole interaction is sufficiently strong, for instance in materials with high optical absorption character.
We can gain further insight by looking at how different points in the Brillouin zone contribute to the conductivity. In Fig. 4, we plot the contributions to the in-plane electrical conductivity from carriers with wavevector lying on high-symmetry lines of the Brillouin zone, the y-coordinate of the dot in figure is positioned according to the carrier's energy. The green dots' radius is proportional to the contribution of such carrier to the aiBTE electrical conductivity, i.e. the contributions from the diagonal part of the Wigner distribution, whereas the red dots radius is proportional to the electrical conductivity contribution from the off-diagonal components of the Wigner distribution. Quantities are computed for the electrical conductivity at temperature of 700K and doping concentration n ¼ 10 18 cm À3 . As one would expect from semiclassical arguments, the diagonal contributions are mostly originating from carriers of energy close to the band gap. The contributions from the off-diagonal components instead extend to carriers that are further from the band gap edge. States close to the band gap at the G point have a sizeable contribution to the off-diagonal conductivity. However, subvalleys overall contribute more significantly to the conductivity, especially due to the larger availability of states for the vertical transitions to take place. We further note that there is not a single transition that dominates the off-diagonal effects. Lastly, we stress that diagonal and off-diagonal contributions have not been drawn to scale with respect to one another. In fact, diagonal contributions are mostly determined by carriers close to the chemical potential and sum up to a conductivity of s aiBTE ¼ 2:9,10 4 S/m. The offdiagonal conductivity Ds ¼ 4:4,10 4 S/m receives contributions from carriers over a much larger energy range. As a result, if the diagonal contributions were drawn to the same scale of the offdiagonal ones, green circles would appear much larger than what is shown; therefore, green radii have been shrunk by a factor about 40 in order to fit in the figure.
So far, we have seen that the off-diagonal contributions to electrical conductivity at energies far from the chemical potential arise from a large presence of available states. However, Eq. (13) depends on several other factors. In Fig. 5a, we plot the average lifetime as a function of carrier's energy, at the same doping density and temperature of Fig. 3. Here we can see that lifetimes tend to be large for states close to the band edges and decrease away from them: to a first approximation, the average lifetime as a function of energy scales inversely proportional to the density of states. Therefore, the decrease of lifetimes is offset by the larger availability of states, which still allows for a sizeable off-diagonal contributions from states away of the band gap. In Fig. 5b we also plot the energy-resolved histogram of the velocity operator jv bb 0 ss 0 ðkÞj in the in-plane direction, where the two axis represent the energies of the carriers ε bs ðkÞ and ε b 0 s 0 ðkÞ respectively, and the color intensity represents the absolute value of the velocity operator matrix element. Clearly, in order to have a sizeable Wigner correction, it is essential for the velocity operator elements to be sufficiently large. However, one can see that the energy values where the velocity elements are largest do not always coincide with the energies of states that contribute most to the conductivity. Therefore, the offdiagonal effects of the Wigner distribution are not easily interpreted in terms of a single factor appearing in Eq. (13): the overall  correction arises due to a combination of large linewidth and velocity (dipoles), as well as high availability of states for the transitions.

Conclusions
First-principles simulations of electronic transport are often limited to the semiclassical approximation, that excludes contributions to transport coming from non-semiclassical particle dynamics. In this study, we have shown that the Moyal equation of motion leads to the Wigner transport equation, an equation of motion for the Wigner function that can be solved using firstprinciples techniques. The formalism has a broader range of applicability than the semiclassical Boltzmann transport equation, and includes off-diagonal contributions to charge and heat currents. In particular, the formalism captures contribution to transport coming from spontaneous vertical transitions between carriers' states, that are especially relevant for narrow-gap semiconductors. Moreover, the semiclassical Boltzmann transport equation is found as a simplified limit to the Wigner transport equation, provided that off-diagonal components are neglected. The formalism is well-suited to first-principles simulations and has been tested with a study of the transport properties of Bi 2 Se 3 . We have shown that, while at high doping concentrations the Boltzmann equation provides a fairly accurate description of transport, it fails at low doping concentrations. At low dopings, spontaneous vertical transitions across the band gap, also known as the Zener tunneling effect, contribute significantly to electronic transport, modifying both electrical conductivity and Seebeck coefficient. Surprisingly, Zener tunneling does not just take place across the states closest to the band gap, but involves states that are significantly further apart in energy, provided that the dipole interaction is sufficiently strong. As a result, we have extended the range of applicability of ab-initio transport simulations to materials where band tunneling couples carriers and a semiclassical description is no longer adequate.

Author contributions
A. C. derived results and performed numerical simulations. A.C. and B. K. conceived the project and wrote the paper.

Data availability
The datasets generated during and/or analyzed during the current study are available from the corresponding author on reasonable request.

Declaration of competing interest
The authors declare that they have no known competing financial interests or personal relationships that could have appeared to influence the work reported in this paper.