+ Affiliations1Department of Physics, University of Tokyo, Bunkyo, Tokyo 113-0033, Japan2Trans-scale Quantum Science Institute, University of Tokyo, Bunkyo, Tokyo 113-0033, Japan
Received November 9, 2021; Accepted March 23, 2022; Published April 22, 2022
We investigate the Seebeck effect in a three-dimensional Dirac electron system based on the linear response theory and Luttinger’s gravitational potential. The Seebeck coefficient S is defined by S = L12/L11T, where T is the temperature, and L11 and L12 are the longitudinal response coefficients of the charge current to the electric field and to the temperature gradient, respectively; L11 is the electric conductivity and L12 is the thermoelectric conductivity. We consider randomly-distributed impurity potentials as the source of the momentum relaxation of electrons and microscopically calculate the relaxation rate and the vertex corrections of L11 and L12 due to the impurities. It is confirmed that L11 and L12 are related through Mott’s formula in low temperatures when the chemical potential lies above the gap (|μ| > Δ), irrespective of the linear dispersion of the Dirac electrons and unconventional energy dependence of the lifetime of electrons. Alternatively, when the chemical potential lies in the bandgap (|μ| < Δ), Seebeck coefficient functions similar to that of conventional semiconductors; its dependences on the chemical potential, μ, and the temperature, T, are partially captured by S ∝ (Δ − μ)/kBT for μ > 0. The Seebeck coefficient takes a relatively large value |S| \( \simeq \) 1.7 mV/K at T \( \simeq \) 4.3 K for Δ = 8 meV by assuming doped bismuth.
©2022 The Physical Society of Japan
1. Introduction Heat flow is accompanied by an electric current in metallic materials. This phenomenon can be understood since electrons or holes, which have charge, are heat carriers in materials. This effect is called the Seebeck effect, which is a thermoelectric effect. When a material is in an open circuit, an electric voltage is generated in the longitudinal direction to the temperature gradient. This thermoelectomotive force, \(\boldsymbol{E}_{\text{emf}}\), is characterized by the Seebeck coefficient S as \(\boldsymbol{E}_{\text{emf}} = S {\boldsymbol{\nabla}} T\), where \({\boldsymbol{\nabla}} T\) is the temperature gradient. It is critical to discover materials with a high Seebeck coefficient for harvesting waste heat and converting it into usable electrical power.
Bismuth was the first material in which the Seebeck effect was observed for the first time.1) Despite its long history, the Seebeck effect in bismuth has yet to be microscopically calculated. According to current knowledge of the bismuth crystal, electron and hole pockets exist in the L and T points in the Brillouin zone, respectively, and electrons near the L point can be described by a three-dimensional (3D) Dirac Hamiltonian2–4) (see Fig. 1), in which the energy dispersion is linear for the momentum measured from the L point. Thus, it is important to study whether this kind of linear dispersion results in unusual temperature- or chemical-potential-dependences of the Seebeck coefficient compared with conventional metals with quadratic momentum dispersion. Further, it is well known that the impurity scattering in Dirac systems, e.g., in graphene, causes unconventional energy dependence of the relaxation rate.5) Therefore, it is necessary to develop a microscopic theory on the relaxation rate for the massive 3D Dirac electrons using the self energy of the Green's functions and the vertex corrections in calculating the electronic conductivities and Seebeck coefficients.
Figure 1. (Color online) Band structure of bismuth along various symmetry lines (a) for the entire energy scale and (b) near the Fermi level. The band structure is based on the tight-binding Hamiltonian provided in Ref. 9, and the symmetry points are provided in Ref. 10.
In this study, we investigate the Seebeck effect in a 3D Dirac electron system. The Seebeck coefficient is defined using longitudinal response coefficients of the charge current to the electric field and to the temperature gradient corresponding to the electric conductivity and thermoelectric conductivity,6,7) respectively. We calculate the response coefficients based on the linear response theory and Luttinger's gravitational potential.8) We consider randomly-distributed impurity potentials as the source of the momentum relaxation of electrons, and the relaxation rate and the vertex corrections of \(L_{11}\) and \(L_{12}\) due to the impurities are microscopically calculated, demonstrating different energy dependences from the quadratic momentum dispersion. We will confirm that the two response coefficients are related through Mott's formula at low temperatures. For the chemical potential dependence of the Seebeck coefficient, S, we find that S has a peak structure when the chemical potential lies in the bandgap at sufficiently low temperatures, whereas it has a monotonic behavior when the system is metallic. We also find a similar peak structure for the temperature dependence of S, as observed in the chemical potential dependence. We demonstrate that these peak structures are partially understood by the phenomenological theory used in semiconductors, which suggests that the Seebeck coefficient is proportional to \((\Delta -\mu)/k_{\text{B}} T\) for \(0 <\mu <\Delta\), where \(2\Delta\) is the bandgap, μ is the chemical potential, \(k_{\text{B}}\) is the Boltzmann constant, and T is the temperature. This behavior of the Seebeck coefficient indicates that the thermoelectric property of the 3D Dirac electrons in the bandgap is the same as that of conventional semiconductors.
The Fermi level may be tuned without changing the band structure by doping to bismuth. Appropriate doping results in the chemical potential inside the bandgap in the L point, where the Seebeck coefficient has the peak structure in the temperature dependence as aforementioned. We evaluated the Seebeck coefficient near the peak and found the value \(|S|\simeq 1.7\) mV/K at \(T\simeq 4.3\) K for \(\Delta = 8\) meV.
Luttinger's gravitational potential approach is mentioned here. In the linear response theory, the Kubo formula is formulated based on external mechanical forces. A mechanical force F couples to the physical quantity \(\hat{A}\) through its Hamiltonian; \(H_{\text{ext}} =\hat{A} F\). The Kubo formula indicates that the response function of a physical quantity \(\hat{B}\) to the external force, F, is given by the correlation function between \(\hat{B}\) and \(\hat{A}\). On the other hand, the temperature gradient is a statistical force that cannot be written by Hamiltonian. The validity of the Kubo formula for responses to a statistical force is not trivial. Then, Luttinger introduced a fictional gravitational potential, which is a mechanical force that couples to the Hamiltonian density.8) Note that the thermal current is coupled by the gradient of the fictional potential. The response to the temperature gradient may then be obtained by applying the Kubo formula to the gravitational potential and using the Einstein's relation, which states that the responses to the mechanical force and to the statistical force are equivalent for the nonequilibrium components.
2. Model and Green's Function Following Ref. 4, we consider the effective (isotropic) Dirac Hamiltonian, \begin{equation} \mathcal{H}_{\text{D}} = \begin{pmatrix} \Delta &i \hbar v \boldsymbol{k}\cdot {\boldsymbol{\sigma}}\\ - i \hbar v \boldsymbol{k}\cdot {\boldsymbol{\sigma}} &-\Delta \end{pmatrix}= - \hbar v \rho_{2} \boldsymbol{k}\cdot {\boldsymbol{\sigma}} + \Delta \rho_{3}, \end{equation} (1) where v is the velocity, \({\boldsymbol{\sigma}} = (\sigma^{x},\sigma^{y},\sigma^{z})\) is the Pauli matrix in spin space, and \(\rho_{i}\) with \(i = 1,2,3\) represents the Pauli matrix in particle–hole space. We use \(\rho_{0}\) and \(\sigma^{0}\) as the unit matrices when emphasizing them. We also consider the point-like impurity potential \(V = u\sum_{i}\rho_{0}\sigma^{0}\delta (\boldsymbol{r}-\boldsymbol{R}_{i})\), where u is the potential strength, and \(\boldsymbol{R}_{i}\) represents the positions of the impurities. The total Hamiltonian is given by \(\mathcal{H} =\mathcal{H}_{\text{D}} + V\).
By taking the average of the positions of impurities,11,12) the retarded self energy in the Born approximation is given by \begin{equation} \varSigma^{\text{R}} (\epsilon) = \frac{n_{\text{i}} u^{2}}{\Omega} \sum_{\boldsymbol{k}} G^{(0)}_{\boldsymbol{k}} (\epsilon + i 0), \end{equation} (2) where \(n_{\text{i}}\) is the impurity concentration (we assume \(n_{\text{i}}\ll 1\)), Ω is the system volume, and the bare Green's function is defined by \(G^{(0)}_{\boldsymbol{k}} (\epsilon + i 0) = (\epsilon +\mu -\mathcal{H}_{\text{D}} + i 0)^{-1}\). The imaginary part is evaluated as follows: \begin{equation} \mathop{\text{Im}}\nolimits \varSigma^{\text{R}} (\epsilon) = - \gamma_{0} (\epsilon+\mu) \rho_{0} - \gamma_{3} (\epsilon+\mu) \rho_{3} \end{equation} (3) with the damping constants \begin{align} \gamma_{0} (\epsilon) &= \frac{\pi}{2} n_{\text{i}} u^{2} \nu (\epsilon), \end{align} (4a) \begin{align} \gamma_{3} (\epsilon) &= \frac{\pi}{2} n_{\text{i}} u^{2} \frac{\Delta}{\epsilon} \nu (\epsilon), \end{align} (4b) where \(\nu (\epsilon)\) is the density of states (DOS), \begin{align} \nu (\epsilon) &= \frac{1}{\Omega} \sum_{\boldsymbol{k}, \eta = \pm} \delta (\epsilon - \eta \varepsilon_{k})\notag\\ &= \frac{|\epsilon|}{2 \pi^{2} \hbar^{3} v^{3}} \sqrt{\epsilon^{2} - \Delta^{2}} \sum_{\eta} \Theta (\eta \epsilon - \Delta) \end{align} (5) with \(\Theta (x)\) being the Heaviside step function.
As mentioned in Ref. 13, it is unclear how the impurity potential is expressed in the basis of the Dirac Hamiltonian. For simplicity, the impurity potential is here assumed to be proportional to \(\rho_{0}\sigma^{0}\) in this study. Further, bismuth has three L points in the Brillouin zone, and the inter-valley scatterings for the point-like impurity potential should be considered.5) However, we neglect the intervalley scatterings at the first step.
Then, the retarded Green's function of this system is given as \begin{equation} G^{\text{R}}_{\boldsymbol{k}} (\epsilon) = \frac{1}{D^{\text{R}}_{\boldsymbol{k}} (\epsilon+\mu)}(g^{\text{R}}_{0} (\epsilon+\mu) + \rho_{2} \boldsymbol{g}^{\text{R}}_{2} (\boldsymbol{k}) \cdot {\boldsymbol{\sigma}} + \rho_{3} g^{\text{R}}_{3} (\epsilon+\mu)) \end{equation} (6) with \begin{align} D^{\text{R}}_{\boldsymbol{k}} (\epsilon) &= (\epsilon + i \gamma_{0} (\epsilon))^{2} - \hbar^{2} v^{2} k^{2} - (\Delta - i \gamma_{3} (\epsilon))^{2}, \end{align} (7a) \begin{align} g^{\text{R}}_{0} (\epsilon) &= \epsilon + i \gamma_{0} (\epsilon), \end{align} (7b) \begin{align} \boldsymbol{g}^{\text{R}}_{2} (\boldsymbol{k}) &= - \hbar v \boldsymbol{k}, \end{align} (7c) \begin{align} g^{\text{R}}_{3} (\epsilon) &= \Delta - i \gamma_{3} (\epsilon). \end{align} (7d)
The denominator can be written as \begin{equation} D^{\text{R}}_{\boldsymbol{k}} (\epsilon) = \prod_{\eta = \pm} (\epsilon - \eta \epsilon_{k} + i \Gamma (\epsilon)), \end{equation} (8) where the eigenenergies are \(\pm\varepsilon_{k} =\pm \sqrt{\hbar^{2} v^{2} k^{2} +\Delta^{2}}\) and \(\Gamma (\epsilon)\) represents the damping of the electron \begin{equation} \Gamma (\epsilon) = \frac{\pi}{2} n_{\text{i}} u^{2} \left(1 + \frac{\Delta^{2}}{\epsilon^{2}} \right) \nu (\epsilon), \end{equation} (9) where we have used the on-shell condition \(\epsilon =\eta \epsilon_{k}\).
Here, we give some comments on the damping constants, \(\gamma_{0} (\epsilon)\) and \(\gamma_{3} (\epsilon)\), the damping \(\Gamma (\epsilon)\), and their energy dependences. There is one kind of damping constant for parabolic dispersion in the absence of magnetization, which is the \(\gamma_{0}\)-type. The two-dimensional Dirac electron system without the mass gap5) has also only the \(\gamma_{0}\)-type damping constant. However, the 3D Dirac electron system with mass gap has the two kinds of damping constants as shown in Eq. (3). This is not because of 3D, but because the system has the mass gap. In general, the damping constants in the Born approximation depends on the physical quantities in the Hamiltonian, such as energy, mass gap, and magnetization.14,15) Hence, the damping of the electron, \(\Gamma (\epsilon)\), has a different dependence on the energy from that of DOS as shown in Fig. 2. Here the lifetime is defined by \(\hbar/2\tau (\epsilon) =\Gamma (\epsilon)\).
Figure 2. (Color online) Energy dependences of DOS, \(\nu (\epsilon)\), the damping, \(\Gamma (\epsilon)\), and the lifetime, \(\tau (\epsilon)\), in the Born approximation.
The velocity operator is given by \begin{equation} v_{i} = \frac{1}{\hbar} \frac{\partial \mathcal{H}_{\text{D}}}{\partial k_{i}} = - v \rho_{2} \sigma^{i} \end{equation} (10) and the electric current operator in the second quantization is obtained as \begin{equation} J_{i} = e v \sum_{\boldsymbol{k}} c^{\dagger}_{\boldsymbol{k}} \rho_{2} \sigma^{i} c_{\boldsymbol{k}}, \end{equation} (11) where \(c_{\boldsymbol{k}}^{(\dagger)}\) is the Fourier component of the field operator, and −e is the electron charge. The thermal current operator in the imaginary time domain is calculated as follows:7,16,17) \begin{equation} \boldsymbol{J}_{\text{Q}} = - \frac{v}{2} \sum_{\boldsymbol{k}} (\dot{c}^{\dagger}_{\boldsymbol{k}} \rho_{2} {\boldsymbol{\sigma}} c_{\boldsymbol{k}} - c^{\dagger}_{\boldsymbol{k}} \rho_{2} {\boldsymbol{\sigma}} \dot{c}_{\boldsymbol{k}}), \end{equation} (12) where \(\dot{c}_{\boldsymbol{k}}^{(\dagger)} = d c_{\boldsymbol{k}}^{(\dagger)}/d\tau\) is the imaginary time derivative (see Appendix A).
3. Seebeck Effect In this study, we consider the longitudinal charge-current response to the electric field, \(E_{x}\), and to the temperature gradient, \(-\nabla_{x} T/T\), which is presented as follows: \begin{equation} \langle J_{x} \rangle = L_{11} E_{x} + L_{12} \left(- \frac{\nabla_{x} T}{T} \right) \end{equation} (13) with the response coefficients, \(L_{11}\) and \(L_{12}\). The Seebeck coefficient, S, is given as12) \begin{equation} S = \frac{L_{12}}{L_{11} T}. \end{equation} (14)
From the linear response theory, the response coefficients are calculated from the following: \begin{align} L_{11} &= \lim_{\omega \to 0} \frac{\chi_{c} (\omega) - \chi_{c} (0)}{i \omega}, \end{align} (15a) \begin{align} L_{12} &= \lim_{\omega \to 0} \frac{\chi_{Q} (\omega) - \chi_{Q} (0)}{i \omega}, \end{align} (15b) where \(\chi_{i}\) with \(i = c, Q\) is evaluated from the corresponding thermal correlation functions \begin{align} \chi_{c} (i \omega_{\lambda}) &= \frac{1}{\Omega} \int_{0}^{\beta} \mathrm{d}\tau\,e^{i \omega_{\lambda} \tau}\langle \mathrm{T}_{\tau} J_{x} (\tau) J_{x} \rangle, \end{align} (16a) \begin{align} \chi_{Q} (i \omega_{\lambda}) &= \frac{1}{\Omega} \int_{0}^{\beta} \mathrm{d}\tau\,e^{i \omega_{\lambda} \tau}\langle \mathrm{T}_{\tau} J_{x} (\tau) J_{\text{Q}, x} \rangle, \end{align} (16b) by taking the analytic continuation \(i\omega_{\lambda}\to \hbar\omega + i 0\), where \(\beta = 1/k_{\text{B}} T\) is the inverse temperature, \(\omega_{\lambda} = 2\pi \lambda k_{\text{B}} T\) with an integer λ is the Matsubara frequency of bosons, T\(_{\tau}\) is the imaginary time ordering operator, and \(J_{x} (\tau) = e^{\tau\mathcal{H}} J_{x} e^{-\tau \mathcal{H}}\) is the Heisenberg representation in the imaginary time.
Precisely, we first calculate the charge current response to Luttinger's gravitational potential based on the linear response theory, and then use the Einstein's relation, which results in the aforementioned formulation.7,8,18–20) We do not need worry about the local equilibrium correction for the longitudinal component,18–20) since it is zero.
Rewriting Eq. (16) using thermal Green's function \(\mathcal{G}_{\boldsymbol{k}} (\tau) = -\langle \mathrm{T}_{\tau} c_{\boldsymbol{k}} (\tau) c^{\dagger}_{\boldsymbol{k}}\rangle\) with the following relation (see Appendix B for the derivation) \begin{align} \langle \mathrm{T}_{\tau} c_{\boldsymbol{k}} (\tau) \dot{c}^{\dagger}_{\boldsymbol{k}} \rangle &= \frac{\mathrm{d}}{\mathrm{d}\tau} \mathcal{G}_{\boldsymbol{k}} (\tau) + \delta (\tau), \end{align} (17a) \begin{align} \langle \mathrm{T}_{\tau} \dot{c}_{\boldsymbol{k}} c^{\dagger}_{\boldsymbol{k}} (\tau) \rangle &= \frac{\mathrm{d}}{\mathrm{d}\tau} \mathcal{G}_{\boldsymbol{k}} (-\tau) - \delta (\tau), \end{align} (17b) with \(\dot{c}_{\boldsymbol{k}}^{(\dagger)} = d c_{\boldsymbol{k}}^{(\dagger)}/d\tau\), and then taking the analytic continuation \(i\omega_{\lambda}\to \hbar\omega + i 0\), we obtain \begin{align} L_{11} &= \int_{-\infty}^{\infty} \mathrm{d}\epsilon\left(- \frac{\partial f}{\partial \epsilon} \right) \sigma (\epsilon + \mu), \end{align} (18) \begin{align} L_{12} &= \frac{1}{- e} \int_{-\infty}^{\infty} \mathrm{d}\epsilon\left(- \frac{\partial f}{\partial \epsilon} \right) \epsilon \sigma (\epsilon + \mu), \end{align} (19) where \(f (\epsilon) = (e^{\beta\epsilon} + 1)^{-1}\), and \begin{equation} \sigma (\epsilon + \mu) = \frac{\hbar e^{2} v^{2}}{4 \pi \Omega}\sum_{\boldsymbol{k}} \mathop{\text{tr}}\nolimits [\rho_{2} \sigma^{x} G^{\text{R}}_{\boldsymbol{k}} (\epsilon) \Lambda_{2, x} G^{\text{A}}_{\boldsymbol{k}} (\epsilon)]. \end{equation} (20) The derivations of Eqs. (18) and (19) with Eq. (20) are given in Appendix C. The vertex \(\Lambda_{2, x}\) represents the full velocity vertex including the ladder-type vertex corrections, \begin{equation} \Lambda_{2, x} = \rho_{2} \sigma^{x} + \frac{n_{\text{i}} u^{2}}{\Omega} \sum_{\boldsymbol{k}} G^{\text{R}}_{\boldsymbol{k}} (\epsilon) \Lambda_{2, x} G^{\text{A}}_{\boldsymbol{k}} (\epsilon), \end{equation} (21) which can be solved as \begin{equation} \Lambda_{2, x} = \frac{1}{1 - U (\epsilon + \mu)} \rho_{2} \sigma^{x} + \mathcal{O} (n_{\text{i}}) \end{equation} (22) with \begin{equation} U (\epsilon) = \frac{1}{3} \frac{\epsilon^{2} - \Delta^{2}}{\epsilon^{2} + \Delta^{2}} \end{equation} (23) as presented in a previous study.13) The detailed calculation from Eqs. (21) to (22) is shown in Appendix D. Note that the vertex corrections have \(\mathcal{O} (n_{\text{i}}^{0})\)-contributions, since \(\sum_{\boldsymbol{k}} G^{\text{R}}_{\boldsymbol{k}} (\epsilon)\rho_{2}\sigma^{x} G^{\text{A}}_{\boldsymbol{k}} (\epsilon)\propto \tau\propto 1/n_{\text{i}} u^{2}\). We then obtain \begin{align} \sigma (\epsilon + \mu) &= \frac{\hbar e^{2} v^{2}}{\pi} \frac{1}{4 n_{\text{i}} u^{2}} \mathop{\text{tr}}\nolimits \left[\rho_{2} \sigma^{x}\left(\frac{n_{\text{i}} u^{2}}{\Omega} \sum_{\boldsymbol{k}} G^{\text{R}}_{\boldsymbol{k}} \Lambda_{2, x} G^{\text{A}}_{\boldsymbol{k}} \right)\right]\notag\\ &= \frac{\hbar e^{2} v^{2}}{\pi} \frac{1}{4 n_{\text{i}} u^{2}} \mathop{\text{tr}}\nolimits [\rho_{2} \sigma^{x} (\Lambda_{2, x} - \rho_{2} \sigma^{x})]\notag\\ &= \frac{\hbar e^{2} v^{2}}{\pi} \frac{1}{n_{\text{i}} u^{2}} \frac{U (\epsilon + \mu)}{1 - U (\epsilon + \mu)}. \end{align} (24) Here, we note that we only calculated the leading order in \(\sigma (\epsilon)\) with respect to \(\mu\tau/\hbar\gg 1\). When calculating the higher order contributions, the terms which include only \(G^{\text{R}}_{\boldsymbol{k}}\) or \(G^{\text{A}}_{\boldsymbol{k}}\) should also be calculated (see Appendix C).
4. Results and Discussion The result is summarized below. The longitudinal charge current response coefficient to the electric field, i.e., the electric conductivity, is calculated as \begin{equation} L_{11} = \int_{-\infty}^{\infty} \mathrm{d}\epsilon \left(- \frac{\partial f_{\text{FD}}}{\partial \epsilon} \right)\sigma (\epsilon), \end{equation} (25) and the longitudinal charge-current response coefficient to the temperature gradient, i.e., the thermoelectric conductivity, is calculated as \begin{equation} L_{12} = \frac{1}{- e} \int_{-\infty}^{\infty} \mathrm{d}\epsilon \left(- \frac{\partial f_{\text{FD}}}{\partial \epsilon} \right) (\epsilon - \mu)\sigma (\epsilon), \end{equation} (26) where \(f_{\text{FD}} (\epsilon) = \{e^{(\epsilon -\mu)/k_{\text{B}} T} + 1 \}^{-1}\) is the Fermi–Dirac distribution function, \begin{equation} \sigma (\epsilon) = \frac{e^{2} v^{2} \tau (\epsilon)}{2} \frac{(\epsilon^{2} - \Delta^{2})(\epsilon^{2} + \Delta^{2})}{\epsilon^{2} (\epsilon^{2} + 2 \Delta^{2})}\nu (\epsilon), \end{equation} (27) \(\tau (\epsilon)\) is the lifetime of electron, and \(\nu (\epsilon)\) is DOS given by Eq. (5). Notably, Eq. (25) reads \(L_{11} =\sigma (\epsilon_{\text{F}})\) at \(T = 0\), implying that \(\sigma (\epsilon)\) describes the electric conductivity at zero temperature.
At the low temperatures, we use the Sommerfeld expansion \begin{equation} \int_{-\infty}^{\infty} \mathrm{d}\epsilon\,H (\epsilon) \left(- \frac{\partial f_{\text{FD}}}{\partial \epsilon} \right) = H (\mu) +\frac{\pi^{2}}{6} H''(\mu) (k_{\text{B}} T)^{2}+ \cdots, \end{equation} (28) which leads to Mott's formula \begin{equation} L_{12} = \frac{\pi^{2}}{3} (k_{\text{B}} T)^{2}\frac{\partial}{\partial \epsilon_{\text{F}}} \left(\frac{1}{- e} \sigma (\epsilon_{\text{F}}) \right), \end{equation} (29) where \(\epsilon_{\text{F}}\) is the Fermi energy (μ is equivalent to the Fermi energy at absolute zero). Therefore, the Seebeck coefficient is expressed as follows: \begin{equation} S = \frac{\pi^{2}}{3} \frac{k_{\text{B}}^{2} T}{- e}\frac{\partial \ln \sigma (\epsilon_{\text{F}})}{\partial \epsilon_{\text{F}}}. \end{equation} (30)
Figure 3 depicts the chemical potential dependences of the response coefficients, \(L_{11}\) and \(L_{12}\), and the Seebeck coefficient, S. The electric conductivity, \(L_{11}\), does not change significantly at the lower temperature (\(k_{\text{B}} T/\Delta\lesssim 0.1\)) as shown in Fig. 3(a), whereas the thermoelectric conductivity, \(L_{12}\), has the relatively large dependence on the temperature [Fig. 3(b)]; thus the Seebeck coefficient defined by Eq. (14) has a peak structure [Fig. 3(c)]. Figures 3(d)–3(f) show the chemical potential dependences in a wider temperature range. For \(k_{\text{B}} T/\Delta\gtrsim 1\), the conductivities \(L_{11}\) and \(L_{12}\) have finite values even in the bandgap, \(|\mu|/\Delta < 1\) [Figs. 3(d) and 3(e)]. This phenomenon implies that thermal fluctuation excite electron and hole pairs, which contribute to the conductivities, as is widely understood in semiconductor physics and superconducting junctions.21) The Seebeck coefficient's peak structure becomes broader as the temperature increases, and the structure disappears for \(k_{\text{B}} T/\Delta\gtrsim 1\), as shown in Fig. 3(f). We discuss the peak structure using a phenomenological analysis in semiconductors after showing the temperature dependences.
Figure 3. (Color online) Chemical potential dependences of the response coefficients, \(L_{11}\), \(L_{12}\), and the Seebeck coefficient, S. (a)–(c) \(L_{11}\), \(L_{12}\), and S at relatively low temperatures, (d)–(f) those in the wider range. The units are given as \(\sigma_{0} = e^{2} v^{2}\tau_{0}\nu_{0}\) and \(\kappa_{0} = -\sigma_{0}\Delta/e\) with \(\tau_{0} =\hbar/\pi n_{\text{i}} u^{2}\nu_{0}\) and \(\nu_{0} =\Delta^{2}/2\pi^{2}\hbar^{3} v^{3}\). The difference between the temperature dependences of \(L_{11}\) and \(L_{12}\) results in the peak structure of S in the bandgap at low temperatures, which disappears above \(k_{\text{B}} T/\Delta\gtrsim 1\). The Seebeck coefficient takes \(S\simeq - 20\,k_{\text{B}}/e\simeq - 1.7\) meV/K at \(k_{\text{B}} T/\Delta = 0.05\) for \(\Delta = 8\) meV, which corresponds to \(T = 4.3\) K.
Here, we evaluate the Seebeck coefficient, S, by assuming the rigid band model for doped bismuth. Assume that one can prepare sufficient doping such that the chemical potential becomes \(\mu\simeq 0.1\Delta\) that lies in the bandgap. For this scenario, the Seebeck coefficient takes \(S\simeq - 20\,k_{\text{B}}/e\simeq - 1.7\) meV/K at \(k_{\text{B}} T/\Delta = 0.05\) with \(\Delta = 8\) meV, corresponding to \(T = 4.3\) K. Note that we do not use the specific values of the velocity and the impurity potential for the quantitative estimation of the Seebeck coefficient. Whereas the electric and thermoelectric conductivities are affected by the velocity and the impurity potential, the Seebeck coefficient in the good metal region is unaffected because the Seebeck coefficient is defined by the ratio of the electric and thermoelectric conductivities, and the velocity and the impurity potential cancel out in the numerator and denominator. We also note that, in actual instances, doping could shorten the electron lifetime, but it is expected to change slightly the Seebeck coefficient, since the Seebeck coefficient is defined by the ratio of the electric and thermoelectric conductivities, both of which are proportional to the lifetime.
Figure 4 shows the temperature dependence of the response coefficients, \(L_{11}\) and \(L_{12}\), and the Seebeck coefficient, S. The electric conductivity, \(L_{11}\), for \(\mu/\Delta = 0\) and 0.9 are zero for sufficiently low temperature since the system has no Fermi surface, whereas \(L_{11}\) for \(\mu/\Delta = 2\) and 5 are metallic even in the absolute zero, as shown in Fig. 4(a). The thermoelectric conductivity, \(L_{12}\), and the Seebeck coefficient S become zero when \(T\to 0\) for all the chemical potentials [Figs. 4(b) and 4(c)]. However, the Seebeck coefficient, when the chemical potential lies in the bandgap (and \(\mu\neq 0\)), has different dependence on the temperature from those when the chemical potential lies out of the bandgap, and has a peak structure [see the blue line (\(\mu/\Delta = 0.9\)) in the inset of Fig. 4(c)]. We briefly discuss the peak structure using the phenomenological analysis. Note that, in the metallic case (\(\mu/\Delta > 1\)), there is no unexpected dependence of S on temperature, and \(L_{12}\) and S in this case can be approximately described by Eqs. (29) and (30), respectively.
Figure 4. (Color online) Temperature dependences of the response coefficients, \(L_{11}\) and \(L_{12}\), and the Seebeck coefficient, S, at various chemical potentials. The units are given in the caption of Fig. 3. At \(\mu = 0\), the response coefficient \(L_{12}\) is zero, and hence the Seebeck coefficient, S, is also zero. The response coefficient \(L_{12}\) and the Seebeck coefficient S become zero when \(T\to 0\) for all the chemical potentials.
Now, we discuss the peak structures in the Seebeck coefficient's chemical potential and temperature dependences. These peak structures appear when the chemical potential lies in the bandgap [Figs. 3(c) and 4(c)]. Here, we regard the system as a semiconductor and apply the phenomenological theory for semiconductors.6,22) The theory says that \(L_{11}\propto \exp[(\Delta -\mu)/k_{\text{B}} T]\) and \(- e L_{12}\propto (\Delta -\mu)\exp[(\Delta -\mu)/k_{\text{B}} T]\), which yields \begin{equation} S_{\text{semi}} = \frac{k_{\text{B}}}{-e} \frac{\Delta - \mu}{k_{\text{B}} T}. \end{equation} (31) Fitting the μ- and T-dependences of the Seebeck coefficient with \(S_{\text{semi}}\) yields Figs. 5(a) and 5(b), respectively. The red lines describing \(S_{\text{semi}}\) are consistent with the μ- and T-dependences of S. This consistency indicates that the Seebeck effect in the 3D Dirac electron system when the chemical potential lies in the bandgap has the similar behavior to that exhibited by semiconductors. This is mainly because the Seebeck effect is related to the charge degree of freedom and not to the spin degree of freedom. It would be different from usual semiconductors for the phenomena related to spin, such as the spin Nernst effect. Note that the deviations for the much lower \(\mu/\Delta\) and \(k_{\text{B}} T/\Delta\) are visible, but the mechanism is unknown.
Figure 5. (Color online) Fittings of the Seebeck coefficient, S, by the phenomenological function \(S_{\text{semi}}\) (a) for the chemical potential dependence and (b) for the temperature dependence.
As briefly described in Ref. 13, we discuss the case of another type of impurity potential which is proportional to \(V'\propto \rho_{0} +\rho_{1}\) obtained from the \(\boldsymbol{k}\cdot \boldsymbol{p}\) theory,23) while assuming that the impurity potential is proportional to \(V\propto \rho_{0}\) in this paper. The self energy for \(V'\) contains only the \(\gamma_{0}\)-term. Moreover, the ladder-type vertex correction for \(V'\) is found to be \(\mathcal{O} (n_{\text{i}})\). However, Eqs. (25) and (26) do not change even in the case of \(V'\), while \(\sigma (\epsilon)\) changes that replacing \(U/(1 - U)\) with U. Thus, only the quantitative difference in \(L_{11}\) and \(L_{12}\) increases in the cases of V and \(V'\). The Seebeck coefficient is more robust when the type of impurity potential is changed than with \(L_{11}\) and \(L_{12}\) since the ratio only contributes to the Seebeck coefficient.
The results are compared with the experiments,24,25) where Ref. 25 is the experiment on the thin film, whereas our calculation is for a bulk system. (The experimental value of the Seebeck coefficient in thin film becomes smaller than that in the bulk system.) The Seebeck coefficient in experiments are in the order of \(-10{\text{–}}{-}100\) µV/K near the room temperature, and our results indicate that the theoretical value is in the order of \(- k_{\text{B}}/e\simeq - 8.6\) µV/K, which is small compared with the experimental value. This may be because we did not consider the contribution from the phonon effects. This remains as a future problem.
Let us discuss the hole contribution near the T point in the Brillouin zone to the Seebeck coefficient. It is believed that the holes can be treated as conventional carriers. The total Seebeck coefficient, \(S_{\text{total}}\), is given by \begin{equation} S_{\text{total}} = \frac{1}{T} \frac{L_{12}^{e} + L_{12}^{h}}{L_{11}^{e} + L_{11}^{h}}, \end{equation} (32) where \(L_{11}^{e}\) and \(L_{12}^{e}\) (\(L_{11}^{h}\) and \(L_{12}^{h}\)) are the electric and thermoelectric conductivities of the electrons near the L point (those of the holes near the T point). The mass of the holes is estimated as \(m_{h}\simeq 0.06 m_{0}{\text{–}}0.22 m_{0}\) and that of the electrons is \(m_{e}\simeq 0.0019 m_{0}{\text{–}}0.027 m_{0}\),26) where \(m_{0}\) is the bare electron mass. The carrier densities of electrons and holes are almost same as \(n_{e}\simeq n_{h}\sim 3\times 10^{17}\) cm−3.9) Therefore, if we assume the same magnitude of the transport lifetime for the electrons and holes, the simple estimation by the Drude conductivity leads to \begin{equation} \frac{L_{11}^{h}}{L_{11}^{e}} = \frac{m_{e}}{m_{h}} \simeq 0.0316, \end{equation} (33) where we use the smallest mass values since the lightest carriers are considered to contribute to the transport. Note that, if we consider the difference of the transport lifetime, the ratio, \(L_{11}^{h}/L_{11}^{e}\), will be further reduced, because the lifetime of Dirac electrons will be longer than that of holes owing to the reduced density of states for Dirac electrons. Here, the intervalley (between L and T points) scattering is neglected assuming that the large momentum exchange reduces the scattering potential. Similarly, for the thermoelectric conductivities, if we assume the Boltzmann-type estimation: \begin{equation} \frac{L_{12}^{h}}{L_{12}^{e}} = \cfrac{\cfrac{1}{+ e} \displaystyle\int \mathrm{d}\epsilon\,f_{\text{FD}} (\epsilon) (\epsilon - \mu) \sigma_{h} (\epsilon)}{\cfrac{1}{- e} \displaystyle\int \mathrm{d}\epsilon\,f_{\text{FD}} (\epsilon) (\epsilon - \mu) \sigma_{e} (\epsilon)} \end{equation} (34) for metallic region, where \(\sigma_{e (h)}\) is electric conductivity of electrons (holes) at zero temperature, the ratio \(L_{12}^{h}/L_{12}^{e}\) is also roughly proportional to \(m_{e}/m_{h}\). Note that the ϵ-dependences of \(\sigma_{e}\) and \(\sigma_{h}\) can differ, but it may be less significant than the mass difference. Thus, Eq. (34) shows that \(S_{\text{total}}\) is approximated by \(S_{\text{total}}\simeq S = L_{12}^{e}/T L_{11}^{e}\).
Since our model is isotropic, we cannot discuss the anisotropy of the Seebeck coefficient, which is thus a future work. Because the experimental data are near the room temperature, phonon effects may be present, although we do not consider them. Hence, to make a quantitative comparison with the experiments, the calculation of phonon effects may be required, which is also a future study.
5. Conclusion We have considered the 3D Dirac electron system and calculated the Seebeck effect based on the linear response theory and Luttinger's gravitational potential. The Seebeck coefficient, S, is defined using the longitudinal responses of the charge current to the electric field and to the temperature gradient, whose response coefficients are, \(L_{11}\) and \(L_{12}\), respectively. We confirm that \(L_{11}\) and \(L_{12}\) for low temperatures are related through Mott's formula. We discuss the dependences of \(L_{11}\), \(L_{12}\), and S on the chemical potential μ and the temperature T. The Seebeck coefficient S, when \(|\mu| <\Delta\) with the bandgap \(2\Delta\), can be fitted with the function \(S_{\text{semi}}\propto (\Delta -\mu)/k_{\text{B}} T\) derived from the phenomenological theory for semiconductors. We also evaluate the Seebeck coefficient by assuming the doped bismuth, which has a relatively large value \(S\simeq 1.7\) mV/K at \(T\simeq 4.3\) K.
Acknowledgments
We thank H. Matsuura for valuable information. This work was supported by Grants-in-Aid for Scientific Research from the Japan Society for the Promotion of Science (Grant No. JP18H01162), and by JST-Mirai Program (Grant No. JPMJMI19A1).
Appendix A: Heat Current Operator In this Appendix, we derive the heat current operator in the Dirac system. Consider the system given by \begin{align} \tilde{\mathscr{H}} &= \mathscr{H} - \mu \mathscr{N}\notag\\ &= \int \mathrm{d}\boldsymbol{r}\,\psi^{\dagger} (\boldsymbol{r}) (- v \rho_{2} \boldsymbol{p}\cdot {\boldsymbol{\sigma}}+ \Delta \rho_{3}- \mu \rho_{0}) \psi (\boldsymbol{r}) \end{align} (A·1) \begin{align} &= \sum_{\boldsymbol{k}} c_{\boldsymbol{k}}^{\dagger} (- \hbar v \rho_{2} \boldsymbol{k}\cdot {\boldsymbol{\sigma}} + \Delta \rho_{3} - \mu \rho_{0}) c_{\boldsymbol{k}}, \end{align} (A·2) where \(\psi^{(\dagger)} (\boldsymbol{r})\) is the field operator of the electron, μ is the chemical potential, and \(\mathscr{N}\) is the electron number. We make the Hamiltonian symmetric as \begin{equation} \tilde{\mathscr{H}} = \int \mathrm{d}\boldsymbol{r}\,Q (\boldsymbol{r}) \end{equation} (A·3) with the heat operator \begin{align} Q (\boldsymbol{r}) &= \psi^{\dagger} (\boldsymbol{r}) \left(\frac{- \hbar v}{2 i} \rho_{2} {\boldsymbol{\sigma}} \cdot (- \overleftarrow{{\boldsymbol{\nabla}}} + \overrightarrow{{\boldsymbol{\nabla}}}) + \Delta \rho_{3} - \mu \rho_{0}\right) \psi (\boldsymbol{r})\notag\\ &\equiv \psi^{\dagger} (\boldsymbol{r}) \bar{\mathcal{H}} \psi (\boldsymbol{r}). \end{align} (A·4)
The continuity equation for the heat operator is calculated as \begin{align} \frac{\mathrm{d}Q (\boldsymbol{r})}{\mathrm{d}t} &= \dot{\psi}^{\dagger} \bar{\mathcal{H}} \psi + \psi^{\dagger} \bar{\mathcal{H}} \dot{\psi}\notag\\ &= - {\boldsymbol{\nabla}} \cdot \left\{- \frac{\hbar v}{2 i} (\dot{\psi}^{\dagger} \rho_{2} {\boldsymbol{\sigma}} \psi - \psi^{\dagger} \rho_{2} {\boldsymbol{\sigma}} \dot{\psi})\right\}\notag\\ &\quad + \dot{\psi}^{\dagger} \left(- \frac{\hbar v}{i} \rho_{2} {\boldsymbol{\sigma}} \cdot \overrightarrow{{\boldsymbol{\nabla}}} + \Delta \rho_{3} - \mu \rho_{0}\right) \psi\notag\\ &\quad + \psi^{\dagger} \left(+ \frac{\hbar v}{i} \rho_{2} {\boldsymbol{\sigma}} \cdot \overleftarrow{{\boldsymbol{\nabla}}} + \Delta \rho_{3}- \mu \rho_{0}\right) \dot{\psi}. \end{align} (A·5) The last two terms cancel out since the following relation \begin{align} i \hbar \frac{\mathrm{d}\psi (\boldsymbol{r})}{\mathrm{d}t} &= \left(- \frac{\hbar v}{i} \rho_{2} {\boldsymbol{\sigma}} \cdot \overrightarrow{{\boldsymbol{\nabla}}}+ \Delta \rho_{3}- \mu \rho_{0}\right) \psi (\boldsymbol{r}), \end{align} (A·6a) \begin{align} - i \hbar \frac{\mathrm{d}\psi^{\dagger} (\boldsymbol{r})}{\mathrm{d}t} &= \psi^{\dagger} (\boldsymbol{r}) \left(+ \frac{\hbar v}{i} \rho_{2} {\boldsymbol{\sigma}} \cdot \overleftarrow{{\boldsymbol{\nabla}}}+ \Delta \rho_{3}- \mu \rho_{0}\right). \end{align} (A·6b) Hence, the continuity equation is obtained as \(\mathrm{d}Q (\boldsymbol{r})/\mathrm{d}t = - {\boldsymbol{\nabla}}\cdot \boldsymbol{j}_{\text{Q}} (\boldsymbol{r})\) with the heat current operator \begin{equation} \boldsymbol{j}_{\text{Q}} (\boldsymbol{r}) = - \frac{\hbar v}{2 i} (\dot{\psi}^{\dagger} \rho_{2} {\boldsymbol{\sigma}} \psi - \psi^{\dagger} \rho_{2} {\boldsymbol{\sigma}} \dot{\psi}) \end{equation} (A·7) in real time domain. Using \(t = - i\hbar \tau\), we obtain the heat current operator in the imaginary time domain as \begin{equation} \boldsymbol{j}_{\text{Q}} (\boldsymbol{r}) = - \frac{v}{2} (\dot{\psi}^{\dagger} \rho_{2} {\boldsymbol{\sigma}} \psi - \psi^{\dagger} \rho_{2} {\boldsymbol{\sigma}} \dot{\psi}). \end{equation} (A·8)
Appendix B: Derivation of Eq. (17) In this Appendix, we derive the relation (17). First, we see the τ-derivative of the thermal Green's function \(\mathcal{G}_{\boldsymbol{k}} (\tau -\tau') = -\langle \mathrm{T}_{\tau} c_{\boldsymbol{k}} (\tau) c_{\boldsymbol{k}}^{\dagger} (\tau')\rangle\). Using the definition of the time-ordering operator in the thermal Green's function, we have \begin{align} \mathcal{G}_{\boldsymbol{k}} (\tau - \tau') &= - \Theta (\tau - \tau') \langle c_{\boldsymbol{k}} (\tau) c_{\boldsymbol{k}}^{\dagger} (\tau') \rangle\notag\\ &\quad + \Theta (\tau' - \tau) \langle c_{\boldsymbol{k}}^{\dagger} (\tau') c_{\boldsymbol{k}} (\tau) \rangle, \end{align} (B·1) and its τ-derivation is obtained as \begin{align} \frac{\partial}{\partial \tau} \mathcal{G}_{\boldsymbol{k}} (\tau - \tau') &= - \delta (\tau - \tau') \langle c_{\boldsymbol{k}} (\tau) c_{\boldsymbol{k}}^{\dagger} (\tau') \rangle\notag\\ &\quad - \delta (\tau' - \tau) \langle c_{\boldsymbol{k}}^{\dagger} (\tau') c_{\boldsymbol{k}} (\tau) \rangle\notag\\ &\quad - \left\langle \mathrm{T}_{\tau} \frac{\partial c_{\boldsymbol{k}} (\tau)}{\partial \tau} c_{\boldsymbol{k}}^{\dagger} (\tau') \right\rangle. \end{align} (B·2) The first two terms in the right hand side read \begin{align} &{-} \delta (\tau - \tau') \langle c_{\boldsymbol{k}} (\tau) c_{\boldsymbol{k}}^{\dagger} (\tau') \rangle - \delta (\tau' - \tau) \langle c_{\boldsymbol{k}}^{\dagger} (\tau') c_{\boldsymbol{k}} (\tau) \rangle\notag\\ &\quad = - \delta (\tau - \tau') \end{align} (B·3) from the anticommutator \(\langle c_{\boldsymbol{k}} (\tau) c_{\boldsymbol{k}}^{\dagger} (\tau) + c_{\boldsymbol{k}}^{\dagger} (\tau) c_{\boldsymbol{k}} (\tau)\rangle = 1\). Hence, \begin{equation} \frac{\partial}{\partial \tau} \mathcal{G}_{\boldsymbol{k}} (\tau - \tau') = - \delta (\tau - \tau') - \langle \mathrm{T}_{\tau} \dot{c}_{\boldsymbol{k}} (\tau) c_{\boldsymbol{k}}^{\dagger} (\tau') \rangle. \end{equation} (B·4) Then, we see the Heisenberg equation for the field operator \(c_{\boldsymbol{k}} (\tau)\), which is given as \begin{equation} \dot{c}_{\boldsymbol{k}} (\tau) = e^{\tau \tilde{\mathscr{H}}} (\tilde{\mathscr{H}} c_{\boldsymbol{k}} - c_{\boldsymbol{k}} \tilde{\mathscr{H}}) e^{- \tau \tilde{\mathscr{H}}}= e^{\tau \tilde{\mathscr{H}}} \dot{c}_{\boldsymbol{k}} e^{- \tau \tilde{\mathscr{H}}}, \end{equation} (B·5) where \(\tilde{\mathscr{H}}\) is defined by Eq. (A·2). Substituting Eq. (B·5) into Eq. (B·4), we have \begin{equation} \frac{\partial}{\partial \tau} \mathcal{G}_{\boldsymbol{k}} (\tau - \tau') = - \delta (\tau - \tau') - \langle \mathrm{T}_{\tau} \dot{c}_{\boldsymbol{k}} c_{\boldsymbol{k}}^{\dagger} (\tau' - \tau) \rangle, \end{equation} (B·6) which leads to Eq. (17b) by replacing \(\tau' -\tau\) with τ.
Similarly, taking the \(\tau'\)-derivative for Eq. (B·1), we have \begin{align} \frac{\partial}{\partial \tau'} \mathcal{G}_{\boldsymbol{k}} (\tau - \tau') &= \delta (\tau - \tau') \langle c_{\boldsymbol{k}} (\tau) c_{\boldsymbol{k}}^{\dagger} (\tau') \rangle\notag\\ &\quad + \delta (\tau' - \tau) \langle c_{\boldsymbol{k}}^{\dagger} (\tau') c_{\boldsymbol{k}} (\tau) \rangle\notag\\ &\quad - \left\langle \mathrm{T}_{\tau} c_{\boldsymbol{k}} (\tau) \frac{\partial c_{\boldsymbol{k}}^{\dagger} (\tau')}{\partial \tau'} \right\rangle\notag\\ &= \delta (\tau - \tau') - \langle \mathrm{T}_{\tau} c_{\boldsymbol{k}} (\tau - \tau') \dot{c}_{\boldsymbol{k}}^{\dagger} \rangle, \end{align} (B·7) which leads to Eq. (17a) by replacing \(\tau -\tau'\) with τ.
Appendix C: Calculation Detail of Correlation Functions Here, we show the detail of the calculation procedures for obtaining Eqs. (18) and (19). Using the thermal Green's function and taking the average on the impurity positions, Eq. (16a) is rewritten as \begin{align} &\chi_{c} (i \omega_{\lambda}) \notag\\ &\quad= - \frac{e^{2} v^{2}}{\beta \Omega} \sum_{n} \sum_{\boldsymbol{k}}\mathop{\text{tr}}\nolimits [\rho_{2} \sigma^{x} \mathcal{G}_{\boldsymbol{k}} (i \epsilon_{n}^{+})\Lambda_{2,x} (i \epsilon_{n}^{+}, i \epsilon_{n}) \mathcal{G}_{\boldsymbol{k}} (i \epsilon_{n})], \end{align} (C·1) where \(i\epsilon_{n}^{+} = i\epsilon_{n} + i\omega_{\lambda}\), and \(\Lambda_{2,x} (i\epsilon_{n}^{+}, i\epsilon_{n})\) is the full velocity vertex including the ladder-type vertex corrections, which is defined by \begin{equation} \Lambda_{2,x} (i \epsilon_{n}^{+}, i \epsilon_{n}) = \rho_{2} \sigma^{x} + \frac{n_{\text{i}} u^{2}}{\Omega} \sum_{\boldsymbol{k}}\mathcal{G}_{\boldsymbol{k}} (i \epsilon_{n}^{+})\Lambda_{2,x} (i \epsilon_{n}^{+}, i \epsilon_{n})\mathcal{G}_{\boldsymbol{k}} (i \epsilon_{n}). \end{equation} (C·2)
The correlation function \(\chi_{Q} (i\omega_{\lambda})\) is also rewritten as follows. For the case without the ladder-type vertex corrections, it reads \begin{equation} \chi_{Q} (i \omega_{\lambda}) = - \frac{e v^{2}}{2 \Omega} \int_{0}^{\beta} \mathrm{d}\tau\,e^{i \omega_{\lambda} \tau} \sum_{\boldsymbol{k}} \mathop{\text{tr}}\nolimits [\rho_{2} \sigma^{x} \langle \mathrm{T}_{\tau} c_{\boldsymbol{k}} (\tau) \dot{c}_{\boldsymbol{k}}^{\dagger} \rangle\rho_{2} \sigma^{x} \mathcal{G}_{\boldsymbol{k}} (- \tau) - \rho_{2} \sigma^{x} \mathcal{G}_{\boldsymbol{k}} (\tau)\rho_{2} \sigma^{x} \langle \mathrm{T}_{\tau} \dot{c}_{\boldsymbol{k}} c_{\boldsymbol{k}}^{\dagger} (\tau) \rangle], \end{equation} (C·3) and using the relation (17), we have \begin{align} \chi_{Q} (i \omega_{\lambda}) &= - \frac{e v^{2}}{2 \Omega} \int_{0}^{\beta} \mathrm{d}\tau\,e^{i \omega_{\lambda} \tau} \sum_{\boldsymbol{k}}\mathop{\text{tr}}\nolimits \left[\rho_{2} \sigma^{x} \left(\frac{\mathrm{d}}{\mathrm{d}\tau} \mathcal{G}_{\boldsymbol{k}} (\tau) \right)\rho_{2} \sigma^{x} \mathcal{G}_{\boldsymbol{k}} (- \tau) - \rho_{2} \sigma^{x} \mathcal{G}_{\boldsymbol{k}} (\tau)\rho_{2} \sigma^{x} \left(\frac{\mathrm{d}}{\mathrm{d}\tau} \mathcal{G}_{\boldsymbol{k}} (- \tau) \right)\right]\notag\\ &= \frac{e v^{2}}{\beta \Omega} \sum_{n} \sum_{\boldsymbol{k}} \mathop{\text{tr}}\nolimits \left[\rho_{2} \sigma^{x} \mathcal{G}_{\boldsymbol{k}} (i \epsilon_{n} + i \omega_{\lambda})\left(i \epsilon_{n} + \frac{i \omega_{\lambda}}{2} \right) \rho_{2} \sigma^{x} \mathcal{G}_{\boldsymbol{k}} (i \epsilon_{n})\right]. \end{align} (C·4) Considering the vertex corrections, we have \begin{equation} \chi_{Q} (i \omega_{\lambda}) = \frac{e v^{2}}{\beta \Omega} \sum_{n} \sum_{\boldsymbol{k}} \mathop{\text{tr}}\nolimits \left[\rho_{2} \sigma^{x} \mathcal{G}_{\boldsymbol{k}} (i \epsilon_{n}^{+})\left(i \epsilon_{n} + \frac{i \omega_{\lambda}}{2} \right) \Lambda_{2, x} (i \epsilon_{n}^{+}, i \epsilon_{n}) \mathcal{G}_{\boldsymbol{k}} (i \epsilon_{n})\right]. \end{equation} (C·5)
Rewriting the Matsubara summation into the contour integral in \(\chi_{c} (i\omega_{\lambda})\) and \(\chi_{Q} (i\omega_{\lambda})\), and taking the analytic continuation \(i\omega_{\lambda}\to \hbar\omega + i 0\), we obtain \begin{align} \chi_{c} (\omega) &= - \frac{i e^{2} v^{2}}{2 \Omega} \sum_{\boldsymbol{k}} \int_{-\infty}^{\infty} \frac{\mathrm{d}\epsilon}{2 \pi}\biggl\{(f (\epsilon_{+}) - f (\epsilon_{-}))\mathop{\text{tr}}\nolimits [\rho_{2} \sigma^{x} G^{\text{R}}_{\boldsymbol{k}} (\epsilon_{+}) \Lambda_{2, x}^{\text{RA}} G^{\text{A}}_{\boldsymbol{k}} (\epsilon_{-})]\notag\\ &\quad + f (\epsilon_{-}) \mathop{\text{tr}}\nolimits [\rho_{2} \sigma^{x} G^{\text{R}}_{\boldsymbol{k}} (\epsilon_{+}) \Lambda_{2, x}^{\text{RR}} G^{\text{R}}_{\boldsymbol{k}} (\epsilon_{-})] - f (\epsilon_{+}) \mathop{\text{tr}}\nolimits [\rho_{2} \sigma^{x} G^{\text{A}}_{\boldsymbol{k}} (\epsilon_{+}) \Lambda_{2, x}^{\text{AA}} G^{\text{A}}_{\boldsymbol{k}} (\epsilon_{-})]\biggr\}, \end{align} (C·6) \begin{align} \chi_{Q} (\omega) &= \frac{i e v^{2}}{2 \Omega} \sum_{\boldsymbol{k}} \int_{-\infty}^{\infty} \frac{\mathrm{d}\epsilon}{2 \pi}\epsilon \biggl\{(f (\epsilon_{+}) - f (\epsilon_{-}))\mathop{\text{tr}}\nolimits [\rho_{2} \sigma^{x} G^{\text{R}}_{\boldsymbol{k}} (\epsilon_{+}) \Lambda_{2, x}^{\text{RA}} G^{\text{A}}_{\boldsymbol{k}} (\epsilon_{-})]\notag\\ &\quad + f (\epsilon_{-}) \mathop{\text{tr}}\nolimits [\rho_{2} \sigma^{x} G^{\text{R}}_{\boldsymbol{k}} (\epsilon_{+}) \Lambda_{2, x}^{\text{RR}} G^{\text{R}}_{\boldsymbol{k}} (\epsilon_{-})] - f (\epsilon_{+}) \mathop{\text{tr}}\nolimits [\rho_{2} \sigma^{x} G^{\text{A}}_{\boldsymbol{k}} (\epsilon_{+}) \Lambda_{2, x}^{\text{AA}} G^{\text{A}}_{\boldsymbol{k}} (\epsilon_{-})]\biggr\}, \end{align} (C·7) where \(\epsilon_{\pm} =\epsilon \pm\hbar \omega/2\), \(f (\epsilon) = (e^{\beta\epsilon} + 1)^{-1}\), \(G^{\text{R}}_{\boldsymbol{k}} (\epsilon)\) are the retarded Green's function given by Eq. (6), \(G^{\text{A}}_{\boldsymbol{k}} (\epsilon)\) is the advanced Green's function obtained by replacing \(i\gamma_{0}\) and \(i\gamma_{3}\) with \(- i\gamma_{0}\) and \(- i\gamma_{3}\) in \(G^{\text{R}}_{\boldsymbol{k}} (\epsilon)\), respectively, and the full velocity vertex \(\Lambda_{2, x}^{\text{RA}}\), \(\Lambda_{2, x}^{\text{RR}}\), and \(\Lambda_{2, x}^{\text{AA}}\) are given as \begin{equation} \Lambda_{2, x}^{\text{XY}} = \rho_{2} \sigma^{x} + \frac{n_{\text{i}} u^{2}}{\Omega} \sum_{\boldsymbol{k}} G^{\text{X}}_{\boldsymbol{k}} (\epsilon_{+}) \Lambda_{2, x}^{\text{XY}} G^{\text{Y}}_{\boldsymbol{k}} (\epsilon_{-}),\quad (\text{X}, \text{Y}\in \{\text{R}, \text{A}\}). \end{equation} (C·8) Note that \(\Lambda_{2, x}^{\text{RA}}\) in this Appendix is equivalent to \(\Lambda_{2, x}\) [Eq. (21)] in the main text. We also note that Eqs. (C·6) and (C·7) are different only in the factor \(-\epsilon/e\).
We extract the ω-linear terms in \(\chi_{c} (\omega)\) and \(\chi_{Q} (\omega)\), which are given as \begin{align} L_{11} &= \frac{\hbar e^{2} v^{2}}{2 \Omega} \sum_{\boldsymbol{k}} \int_{-\infty}^{\infty} \frac{\mathrm{d}\epsilon}{2 \pi}\biggl\{\left(- \frac{\partial f}{\partial \epsilon} \right) \mathop{\text{tr}}\nolimits [\rho_{2} \sigma^{x} G^{\text{R}}_{\boldsymbol{k}} (\epsilon) \Lambda_{2, x}^{\text{RA}} G^{\text{A}}_{\boldsymbol{k}} (\epsilon) - \mathbb{R}\mathrm{e}\{\rho_{2} \sigma^{x} G^{\text{R}}_{\boldsymbol{k}} (\epsilon) \Lambda_{2, x}^{\text{RR}} G^{\text{R}}_{\boldsymbol{k}} (\epsilon) \}]\notag\\ &\quad - i f (\epsilon) (\partial_{\epsilon} - \partial_{\epsilon'}) \mathop{\text{tr}}\nolimits [\mathbb{I}\mathrm{m}\{\rho_{2} \sigma^{x} G^{\text{R}}_{\boldsymbol{k}} (\epsilon) \Lambda_{2, x}^{\text{RR}} G^{\text{R}}_{\boldsymbol{k}} (\epsilon') \}] |_{\epsilon' \to \epsilon}\biggr\}, \end{align} (C·9) \begin{align} L_{12} &= - \frac{\hbar e v^{2}}{2 \Omega} \sum_{\boldsymbol{k}} \int_{-\infty}^{\infty} \frac{\mathrm{d}\epsilon}{2 \pi}\epsilon \biggl\{\left(- \frac{\partial f}{\partial \epsilon} \right) \mathop{\text{tr}}\nolimits [\rho_{2} \sigma^{x} G^{\text{R}}_{\boldsymbol{k}} (\epsilon) \Lambda_{2, x}^{\text{RA}} G^{\text{A}}_{\boldsymbol{k}} (\epsilon) - \mathbb{R}\mathrm{e}\{\rho_{2} \sigma^{x} G^{\text{R}}_{\boldsymbol{k}} (\epsilon) \Lambda_{2, x}^{\text{RR}} G^{\text{R}}_{\boldsymbol{k}} (\epsilon) \}]\notag\\ &\quad - i f (\epsilon) (\partial_{\epsilon} - \partial_{\epsilon'}) \mathop{\text{tr}}\nolimits [\mathbb{I}\mathrm{m}\{\rho_{2} \sigma^{x} G^{\text{R}}_{\boldsymbol{k}} (\epsilon) \Lambda_{2, x}^{\text{RR}} G^{\text{R}}_{\boldsymbol{k}} (\epsilon') \}] |_{\epsilon' \to \epsilon}\biggr\}, \end{align} (C·10) where \(\mathbb{R}\mathrm{e}[\cdots]\) and \(\mathbb{I}\mathrm{m}[\cdots]\) are defined as \begin{equation} \mathbb{R}\mathrm{e}\{P^{\text{R}} \} = \frac{1}{2} [P^{\text{R}} + P^{\text{A}}],\quad \mathbb{I}\mathrm{m}\{P^{\text{R}} \} = \frac{1}{2 i} [P^{\text{R}} - P^{\text{A}}]. \end{equation} (C·11) Then, we neglect the terms which include only \(G^{\text{R}}_{\boldsymbol{k}}\) or \(G^{\text{A}}_{\boldsymbol{k}}\), since they contribute only in the higher orders with respect to \(q =\mu \tau/\hbar\) for \(q\gg 1\). Finally, we have Eqs. (18) and (19) with Eq. (20) in the main text.
Appendix D: Calculation Detail of Vertex Corrections Here, we show the calculation of the ladder-type vertex corrections. Equation (21) reads \begin{equation} \Lambda_{2,x} = \rho_{2} \sigma^{x}+ \frac{n_{\text{i}} u^{2}}{\Omega} \sum_{\boldsymbol{k}} G^{\text{R}}_{\boldsymbol{k}} (\epsilon) \rho_{2} \sigma^{x} G^{\text{A}}_{\boldsymbol{k}} (\epsilon)+ \left(\frac{n_{\text{i}} u^{2}}{\Omega} \right)^{2} \sum_{\boldsymbol{k}, \boldsymbol{k}'} G^{\text{R}}_{\boldsymbol{k}} (\epsilon) G^{\text{R}}_{\boldsymbol{k}'} (\epsilon) \rho_{2} \sigma^{x} G^{\text{A}}_{\boldsymbol{k}'} (\epsilon) G^{\text{A}}_{\boldsymbol{k}} (\epsilon)+ \cdots. \end{equation} (D·1) Firstly, we calculate the second term in the right hand side, which is computed as \begin{equation} \frac{n_{\text{i}} u^{2}}{\Omega} \sum_{\boldsymbol{k}} G^{\text{R}}_{\boldsymbol{k}} (\epsilon) \rho_{2} \sigma^{x} G^{\text{A}}_{\boldsymbol{k}} (\epsilon) = U (\epsilon + \mu) \rho_{2} \sigma^{x}- V (\epsilon + \mu) \rho_{1} \sigma^{x}, \end{equation} (D·2) where \(U (\epsilon)\) and \(V (\epsilon)\) are given by \begin{align} U (\epsilon) &= \frac{1}{4} \frac{n_{\text{i}} u^{2}}{\Omega} \sum_{\boldsymbol{k}} \mathop{\text{tr}}\nolimits [G^{\text{R}}_{\boldsymbol{k}} (\epsilon - \mu) \rho_{2} \sigma^{x} G^{\text{A}}_{\boldsymbol{k}} (\epsilon - \mu) \rho_{2} \sigma^{x}], \end{align} (D·3) \begin{align} V (\epsilon) &= \frac{1}{4} \frac{n_{\text{i}} u^{2}}{\Omega} \sum_{\boldsymbol{k}} \mathop{\text{tr}}\nolimits [G^{\text{R}}_{\boldsymbol{k}} (\epsilon - \mu) \rho_{2} \sigma^{x} G^{\text{A}}_{\boldsymbol{k}} (\epsilon - \mu) \rho_{1} \sigma^{x}]. \end{align} (D·4) Using \(U (\epsilon)\) and \(V (\epsilon)\), the third term in Eq. (D·1) is calculated as \begin{align} &\left(\frac{n_{\text{i}} u^{2}}{\Omega} \right)^{2} \sum_{\boldsymbol{k}, \boldsymbol{k}'}G^{\text{R}}_{\boldsymbol{k}} (\epsilon) G^{\text{R}}_{\boldsymbol{k}'} (\epsilon) \rho_{2} \sigma^{x} G^{\text{A}}_{\boldsymbol{k}'} (\epsilon) G^{\text{A}}_{\boldsymbol{k}} (\epsilon)\notag\\ &\quad = \frac{n_{\text{i}} u^{2}}{\Omega} \sum_{\boldsymbol{k}} G^{\text{R}}_{\boldsymbol{k}} (\epsilon) (U (\epsilon + \mu) \rho_{2} \sigma^{x}- V (\epsilon + \mu) \rho_{1} \sigma^{x}) G^{\text{A}}_{\boldsymbol{k}} (\epsilon)\notag\\ &\quad = U (\epsilon + \mu) \{U (\epsilon + \mu) \rho_{2} \sigma^{x}- V (\epsilon + \mu) \rho_{1} \sigma^{x}\}- V (\epsilon + \mu)\frac{n_{\text{i}} u^{2}}{\Omega} \sum_{\boldsymbol{k}} G^{\text{R}}_{\boldsymbol{k}} (\epsilon) \rho_{1} \sigma^{x} G^{\text{A}}_{\boldsymbol{k}} (\epsilon)\notag\\ &\quad = U (\epsilon + \mu) \{U (\epsilon + \mu) \rho_{2} \sigma^{x}- V (\epsilon + \mu) \rho_{1} \sigma^{x}\}- V (\epsilon + \mu)\{2 U (\epsilon + \mu) \rho_{2} \sigma^{x}+ V (\epsilon + \mu) \rho_{1} \sigma^{x}\}\notag\\ &\quad = \begin{pmatrix} 0 &1 \end{pmatrix}\begin{pmatrix} 2 U (\epsilon + \mu) &V (\epsilon + \mu)\\ - V (\epsilon + \mu) &U (\epsilon + \mu) \end{pmatrix}^{2} \begin{pmatrix} \rho_{1} \sigma^{x}\\ \rho_{2} \sigma^{x} \end{pmatrix}. \end{align} (D·5) Here, we used \begin{equation} \frac{n_{\text{i}} u^{2}}{\Omega} \sum_{\boldsymbol{k}} G^{\text{R}}_{\boldsymbol{k}} (\epsilon) \rho_{1} \sigma^{x} G^{\text{A}}_{\boldsymbol{k}} (\epsilon) = 2 U (\epsilon + \mu) \rho_{1} \sigma^{x}+ V (\epsilon + \mu) \rho_{2} \sigma^{x}. \end{equation} (D·6) Similarly, the first and second terms in Eq. (D·1) are also written as \begin{align} \rho_{2} \sigma^{x} &= \begin{pmatrix} 0 &1 \end{pmatrix}\begin{pmatrix} \rho_{1} \sigma^{x}\\ \rho_{2} \sigma^{x} \end{pmatrix}, \end{align} (D·7) \begin{align} \frac{n_{\text{i}} u^{2}}{\Omega} \sum_{\boldsymbol{k}} G^{\text{R}}_{\boldsymbol{k}} (\epsilon) \rho_{2} \sigma^{x} G^{\text{A}}_{\boldsymbol{k}} (\epsilon) &= \begin{pmatrix} 0 &1 \end{pmatrix}\begin{pmatrix} 2 U (\epsilon + \mu) &V (\epsilon + \mu)\\ - V (\epsilon + \mu) &U (\epsilon + \mu) \end{pmatrix}\begin{pmatrix} \rho_{1} \sigma^{x}\\ \rho_{2} \sigma^{x} \end{pmatrix}. \end{align} (D·8) From these, we have \begin{align} \Lambda_{2, x} &= \begin{pmatrix} 0 &1 \end{pmatrix}\left\{ \begin{pmatrix} 1 &0\\ 0 &1 \end{pmatrix}+ \begin{pmatrix} 2 U (\epsilon + \mu) &V (\epsilon + \mu)\\ - V (\epsilon + \mu) &U (\epsilon + \mu) \end{pmatrix}+ \begin{pmatrix} 2 U (\epsilon + \mu) &V (\epsilon + \mu)\\ - V (\epsilon + \mu) &U (\epsilon + \mu) \end{pmatrix}^{2}{} + \cdots\right\} \begin{pmatrix} \rho_{1} \sigma^{x}\\ \rho_{2} \sigma^{x} \end{pmatrix}\notag\\ &= \begin{pmatrix} 0 &1 \end{pmatrix}\{A^{0}+ A+ A^{2}+ \cdots\} \begin{pmatrix} \rho_{1} \sigma^{x}\\ \rho_{2} \sigma^{x} \end{pmatrix}, \end{align} (D·9) where \begin{equation} A = \begin{pmatrix} 2 U (\epsilon + \mu) &V (\epsilon + \mu)\\ - V (\epsilon + \mu) &U (\epsilon + \mu) \end{pmatrix}. \end{equation} (D·10) Introducing the matrix P which diagonalize A and the eigenvalue matrix \(E =\mathop{\text{diag}}\nolimits (\lambda_{1},\lambda_{2})\); \(A = P (P^{-1} A P) P^{-1} = P E P^{-1}\), we find \begin{align} A^{0} + A + A^{2} + \cdots &= P (1 + E + E^{2} + \cdots) P^{-1} \notag\\ &= P \begin{pmatrix} \dfrac{1}{1 - \lambda_{1}} &0\\ 0 &\dfrac{1}{1 - \lambda_{2}} \end{pmatrix}P^{-1}. \end{align} (D·11) Hence, we obtain \begin{equation} \Lambda_{2, x} = \begin{pmatrix} 0 &1 \end{pmatrix}P \begin{pmatrix} \dfrac{1}{1 - \lambda_{1}} &0\\ 0 &\dfrac{1}{1 - \lambda_{2}} \end{pmatrix}P^{-1} \begin{pmatrix} \rho_{1} \sigma^{x}\\ \rho_{2} \sigma^{x} \end{pmatrix}. \end{equation} (D·12)
Before showing the forms of P and E, we calculate \(U (\epsilon)\) and \(V (\epsilon)\). \begin{align} U (\epsilon) &= \frac{n_{\text{i}} u^{2}}{\Omega} \sum_{\boldsymbol{k}} \frac{1}{| D^{\text{R}}_{\boldsymbol{k}} (\epsilon) |^{2}}| g^{\text{R}}_{2,x} (\boldsymbol{k}) |^{2}, \end{align} (D·13) \begin{align} V (\epsilon) &= \frac{n_{\text{i}} u^{2}}{\Omega} \sum_{\boldsymbol{k}} \frac{1}{| D^{\text{R}}_{\boldsymbol{k}} (\epsilon) |^{2}}(g^{\text{R}}_{3} (\epsilon) g^{\text{A}}_{0} (\epsilon) - g^{\text{R}}_{0} (\epsilon) g^{\text{A}}_{3} (\epsilon)). \end{align} (D·14) From the symmetry of the system, \(| g^{\text{R}}_{2,x} (\boldsymbol{k}) |^{2} =\hbar^{2} v^{2} k^{2}/3 = (\epsilon_{k}^{2} -\Delta^{2})/3\). We express \(D^{\text{R}}_{\boldsymbol{k}} (\epsilon) = D' + i D''\) with \begin{align} D' &= (\epsilon - \epsilon_{k}) (\epsilon + \epsilon_{k}) + \mathcal{O} (n_{\text{i}}^{2}), \end{align} (D·15) \begin{align} D''&= 2 (\epsilon \gamma_{0} (\epsilon) + \Delta \gamma_{3} (\epsilon)), \end{align} (D·16) and use the following approximations \begin{align} \frac{1}{| D^{\text{R}}_{\boldsymbol{k}} (\epsilon) |^{2}} &\simeq \frac{\pi}{| D''|} \delta (D')\notag\\ &\simeq \frac{\pi}{4} \frac{1}{| \epsilon \gamma_{0} (\epsilon) + \Delta \gamma_{3} (\epsilon) |} \frac{1}{|\epsilon|} \sum_{\eta = \pm} \delta (\epsilon - \eta \epsilon_{k}), \end{align} (D·17) where, we assumed small \(n_{\text{i}}\) and approximated \(D''/((D')^{2} + (D'')^{2})\) by the δ-function in the first line, and dropped \(\gamma_{0}^{2}\) and \(\gamma_{3}^{2}\) in the second line. From these, we have \begin{align} U (\epsilon) &= \frac{1}{3} \frac{\epsilon^{2} - \Delta^{2}}{\epsilon^{2} + \Delta^{2}} + \mathcal{O} (n_{\text{i}}^{2}), \end{align} (D·18) \begin{align} V (\epsilon) &= n_{\text{i}} u^{2} \frac{\pi \Delta}{\epsilon^{2} + \Delta^{2}} + \mathcal{O} (n_{\text{i}}^{3}). \end{align} (D·19)
For small \(n_{\text{i}}\), which leads to \(\{U (\epsilon +\mu) \}^{2} - 4 \{V (\epsilon +\mu) \}^{2}\simeq \{U (\epsilon +\mu) \}^{2}\), the eigenvalue is obtained as \(\lambda_{1} = 2 U (\epsilon +\mu)\) and \(\lambda_{2} = U (\epsilon +\mu)\). The diagonalization matrix P is calculated as \begin{align} &P \simeq \begin{pmatrix} U (\epsilon + \mu) &-V (\epsilon + \mu)\\ - V (\epsilon + \mu) &U (\epsilon + \mu) \end{pmatrix},\notag\\ &P^{-1} \simeq \frac{1}{\{U (\epsilon + \mu) \}^{2}} \begin{pmatrix} U (\epsilon + \mu) &V (\epsilon + \mu)\\ V (\epsilon + \mu) &U (\epsilon + \mu) \end{pmatrix}. \end{align} (D·20) Hence, \begin{align} \Lambda_{2, x} &\simeq \begin{pmatrix} 0 &1 \end{pmatrix}\begin{pmatrix} \dfrac{1}{1 - 2 U (\epsilon + \mu)} &\dfrac{1}{\{1 - 2 U (\epsilon + \mu)\} \{1 - U (\epsilon + \mu) \}}\\ -\dfrac{1}{\{1 - 2 U (\epsilon + \mu)\} \{1 - U (\epsilon + \mu) \}} &\dfrac{1}{1 - U (\epsilon + \mu)} \end{pmatrix}\begin{pmatrix} \rho_{1} \sigma^{x}\\ \rho_{2} \sigma^{x} \end{pmatrix}\notag\\ &= \frac{1}{1 - U (\epsilon + \mu)} \rho_{2} \sigma^{x}. \end{align} (D·21)
References
- 1 T. J. Seebeck, Abh. Akad. Wiss. Berlin 289, 1820 (1822). Google Scholar
- 2 M. H. Cohen and E. I. Blount, Philos. Mag. 5, 115 (1960). 10.1080/14786436008243294 Crossref, Google Scholar
- 3 P. A. Wolff, J. Phys. Chem. Solids 25, 1057 (1964). 10.1016/0022-3697(64)90128-3 Crossref, Google Scholar
- 4 Y. Fuseya, M. Ogata, and H. Fukuyama, J. Phys. Soc. Jpn. 81, 093704 (2012). 10.1143/JPSJ.81.093704 Link, Google Scholar
- 5 N. H. Shon and T. Ando, J. Phys. Soc. Jpn. 67, 2421 (1998). 10.1143/JPSJ.67.2421 Link, Google Scholar
- 6 K. Behnia, Fundamentals of Thermoelectricity (Oxford University Press, New York, 2015). Crossref, Google Scholar
- 7 M. Ogata and H. Fukuyama, J. Phys. Soc. Jpn. 88, 074703 (2019). 10.7566/JPSJ.88.074703 Link, Google Scholar
- 8 J. Luttinger, Phys. Rev. 135, A1505 (1964). 10.1103/PhysRev.135.A1505 Crossref, Google Scholar
- 9 Y. Liu and R. E. Allen, Phys. Rev. B 52, 1566 (1995). 10.1103/PhysRevB.52.1566 Crossref, Google Scholar
- 10 L. M. Falicov and S. Golin, Phys. Rev. 137, A871 (1965). 10.1103/PhysRev.137.A871 Crossref, Google Scholar
- 11 W. Kohn and J. M. Luttinger, Phys. Rev. 108, 590 (1957). 10.1103/PhysRev.108.590 Crossref, Google Scholar
- 12 G. D. Mahan, Many-Particle Physics (Springer US, Boston, MA, 2000). Crossref, Google Scholar
- 13 T. Fukazawa, H. Kohno, and J. Fujimoto, J. Phys. Soc. Jpn. 86, 094704 (2017). 10.7566/JPSJ.86.094704 Link, Google Scholar
- 14 J. Fujimoto and H. Kohno, Phys. Rev. B 90, 214418 (2014). 10.1103/PhysRevB.90.214418 Crossref, Google Scholar
- 15 J. Fujimoto, Phys. Rev. B 97, 104421 (2018). 10.1103/PhysRevB.97.104421 Crossref, Google Scholar
- 16 M. Jonson and G. D. Mahan, Phys. Rev. B 42, 9350 (1990). 10.1103/PhysRevB.42.9350 Crossref, Google Scholar
- 17 H. Kontani, Phys. Rev. B 67, 014408 (2003). 10.1103/PhysRevB.67.014408 Crossref, Google Scholar
- 18 L. Smrcka and P. Streda, J. Phys. C 10, 2153 (1977). 10.1088/0022-3719/10/12/021 Crossref, Google Scholar
- 19 N. R. Cooper, B. I. Halperin, and I. M. Ruzin, Phys. Rev. B 55, 2344 (1997). 10.1103/PhysRevB.55.2344 Crossref, Google Scholar
- 20 H. Kohno, Y. Hiraoka, M. Hatami, and G. E. W. Bauer, Phys. Rev. B 94, 104417 (2016). 10.1103/PhysRevB.94.104417 Crossref, Google Scholar
- 21 M. Tinkham, Introduction to Superconductivity (Dover, New York, 2004). Google Scholar
- 22 N. F. Mott and E. A. Davis, Electronic Processes in Non-Crystalline Materials, Oxford Classic Texts in the Physical Sciences (Oxford University Press, Oxford, New York, 2012). Google Scholar
- 23 K. Sakai, C. Ishii, and H. Fukuyama, J. Phys. Soc. Jpn. 50, 3590 (1981). 10.1143/JPSJ.50.3590 Link, Google Scholar
- 24 B. S. Chandrasekhar, J. Phys. Chem. Solids 11, 268 (1959). 10.1016/0022-3697(59)90225-2 Crossref, Google Scholar
- 25 V. D. Das and N. Soundararajan, Phys. Rev. B 35, 5990 (1987). 10.1103/PhysRevB.35.5990 Crossref, Google Scholar
- 26 Y. Fuseya, Z. Zhu, B. Fauqué, W. Kang, B. Lenoir, and K. Behnia, Phys. Rev. Lett. 115, 216401 (2015). 10.1103/PhysRevLett.115.216401 Crossref, Google Scholar
, and 
