Breakdown of the Wiedemann-Franz law in AB-stacked bilayer graphene

TL;DR

This paper develops a simple theoretical model for thermoelectric transport in bilayer graphene, revealing how disorder influences the violation of the Wiedemann-Franz law near charge neutrality.

Contribution

It provides a detailed analysis of thermoelectric properties in bilayer graphene, highlighting the effects of disorder on Wiedemann-Franz law violations, which was less understood compared to single-layer graphene.

Findings

Thermal resistivity approaches zero at charge neutrality without disorder.

Electric resistivity jumps to a finite value at neutrality.

Disorder introduces a doping window with strong Wiedemann-Franz law violation.

Abstract

We present a simple theory of thermoelectric transport in bilayer graphene and report our results for the electrical resistivity, the thermal resistivity, the Seebeck coefficient, and the Wiedemann-Franz ratio as functions of doping density and temperature. In the absence of disorder, the thermal resistivity tends to zero as the charge neutrality point is approached; the electric resistivity jumps from zero to an intrinsic finite value, and the Seebeck coefficient diverges in the same limit. Even though these results are similar to those obtained for single-layer graphene, their derivation is considerably more delicate. The singularities are removed by the inclusion of a small amount of disorder, which leads to the appearance of a "window" of doping densities $0<n<n_c$ (with $n_c$ tending to zero in the zero-disorder limit) in which the Wiedemann-Franz law is severely violated.

Equations

$${\mathbf{j}}_n=\sum _{\mathbf{k},\gamma }{\mathbf{v}}_{\mathbf{k},\gamma }\delta f_{\mathbf{k},\gamma },\text{leavevmode}\text{leavevmode}{\mathbf{j}}_s=\sum _{\mathbf{k},\gamma }\beta {\widetilde{ϵ}}_{\mathbf{k},\gamma }{\mathbf{v}}_{\mathbf{k},\gamma }\delta f_{\mathbf{k},\gamma },$$

$$\mbox{\boldmath$\rho$}=\left(\begin{array}[]{cc}\bar{\rho}_{el}+\bar{Q}^{2}\bar{\rho}_{th}&\bar{Q}\bar{\rho}_{th}\\ \bar{Q}\bar{\rho}_{th}&\bar{\rho}_{th}\end{array}\right)\,,$$

$$D_{ij}=\frac{1}{2}\sum _{\mathbf{k},\gamma }{f}_{0\mathbf{k}\gamma }^{\prime }{\mathbf{u}}_{\mathbf{k},\gamma }^i\cdot {\mathbf{u}}_{\mathbf{k},\gamma }^j\text{leavevmode},$$

$$-f^{\prime}_{0{\bf k},\gamma}{\bf v}_{{\bf k},\gamma}\cdot\left[e{\bf E}+\beta\tilde{\epsilon}_{{\bf k},\gamma}k_{B}\mbox{\boldmath$\nabla$}T\right]=I_{{\bf k},\gamma}\leavevmode\nobreak\ .$$

$${\tilde{I}}_{ij}=-\frac{\beta }{4\pi }\sum _{\mathbf{q}}{\int }_{-\infty }^{\infty }d\omega \frac{|V(q){|}^2(\Im {\Pi }_i^1\Im {\Pi }_j^1-\Im {\Pi }^0\Im {\Pi }_{ij}^2)}{{\sinh }^2(\beta \hslash \omega /2)},$$

$${\Pi }^{\alpha }=4\sum _{\gamma ,{\gamma }^{\prime }}\sum _{\mathbf{k}}\frac{F_{\mathbf{k},\mathbf{k}+\mathbf{q}}^{\gamma {\gamma }^{\prime }}(f_{0\mathbf{k},\gamma }-f_{0\mathbf{k}+\mathbf{q},{\gamma }^{\prime }})}{{ϵ}_{\mathbf{k},\gamma }-{ϵ}_{\mathbf{k}+\mathbf{q},{\gamma }^{\prime }}+\hslash \omega +i{0}^{+}}{M}_{\mathbf{k},\mathbf{k}+\mathbf{q},\gamma ,{\gamma }^{\prime }}^{\alpha }.$$

$$\delta f_{\mathbf{k},\gamma }=f_{0\mathbf{k},\gamma }^{\prime }\left[\hslash v^2\mathbf{k}\cdot {\mathbf{p}}_k/t+{\mathbf{v}}_{\mathbf{k},\gamma }\cdot ({\mathbf{p}}_n+\beta {\widetilde{ϵ}}_{\mathbf{k},\gamma }{\mathbf{p}}_s)\right],$$

$$I_{jj^{\prime}}=\left(\begin{array}[]{ccc}\eta&0&0\\ 0&\tilde{I}_{11}&\tilde{I}_{12}\\ 0&\tilde{I}_{12}&\tilde{I}_{22}\end{array}\right)_{jj^{\prime}}$$

$${\mbox{\boldmath$\rho$}}=\frac{\left(\begin{array}[]{cc}[D_{02}]^{2}&-D_{01}D_{02}\\ -D_{01}D_{02}&[D_{01}]^{2}\end{array}\right)}{[D_{01}]^{2}\tilde{\sigma}_{22}+[D_{02}]^{2}\tilde{\sigma}_{11}-2D_{01}D_{02}[\tilde{\sigma}_{12}]^{2}}$$

$${\bar{\rho }}_{th}(n\ne 0)={\rho }_{22},\text{leavevmode}{\bar{\rho }}_{el}(n\ne 0)=0,\text{leavevmode}\bar{Q}=-D_{02}/D_{01}\text{leavevmode}.$$

$$WF={\left(\frac{\Gamma }{(D_{01}/D_{02}{)}^2+{\Gamma }^2}\right)}^2,\text{leavevmode}\text{leavevmode}{\Gamma }^2=\frac{{\bar{\rho }}_{th,D}}{{\bar{\rho }}_{el,C}}.$$

$${\tilde{I}}_{\mathbf{k},\gamma }$$

$$\delta f_{\mathbf{k},\gamma }=f_{0\mathbf{k},\gamma }^{\prime }{\mathbf{v}}_{\mathbf{k},\gamma }\cdot ({\mathbf{p}}_n+\beta {\widetilde{ϵ}}_{\mathbf{k},\gamma }{\mathbf{p}}_s)$$

$$\delta ({ϵ}_{\mathbf{k},\gamma }+{ϵ}_{{\mathbf{k}}^{\prime },\eta }-{ϵ}_{\mathbf{k}-\mathbf{q},{\gamma }^{\prime }}-{ϵ}_{\mathbf{k}+\mathbf{q},{\eta }^{\prime }})=\hslash {\int }_{-\infty }^{\infty }d\omega \delta ({ϵ}_{{\mathbf{k}}^{\prime },\eta }-{ϵ}_{{\mathbf{k}}^{\prime }+\mathbf{q},{\eta }^{\prime }}+\hslash \omega )\delta ({ϵ}_{\mathbf{k},\gamma }-{ϵ}_{\mathbf{k}-\mathbf{q},{\gamma }^{\prime }}-\hslash \omega )$$

$$f_{0\mathbf{k},\gamma }(1-f_{0\mathbf{k}\pm \mathbf{q},{\gamma }^{\prime }})\delta ({ϵ}_{\mathbf{k},\gamma }-{ϵ}_{\mathbf{k}\pm \mathbf{q},{\gamma }^{\prime }}\pm \hslash \omega )=\frac{f_{0\mathbf{k},\gamma }-f_{0\mathbf{k}\pm \mathbf{q},{\gamma }^{\prime }}}{\mp 2e^{\pm \beta \hslash \omega /2}\sinh (\beta \hslash \omega /2)}\delta ({ϵ}_{\mathbf{k},\gamma }-{ϵ}_{\mathbf{k}\pm \mathbf{q},{\gamma }^{\prime }}\pm \hslash \omega ),$$

$${\tilde{I}}_{\mathbf{k},\gamma }=\frac{2\pi \beta }{4}\sum _{{\gamma }^{\prime },\eta ,{\eta }^{\prime }}\sum _{{\mathbf{k}}^{\prime }}\sum _{\mathbf{q}}{\int }_{-\infty }^{\infty }d\omega \frac{|V(\mathbf{q}){|}^2}{{\sinh }^2(\beta \hslash \omega /2)}\times$$

$$F_{\mathbf{k},\mathbf{k}-\mathbf{q}}^{\gamma {\gamma }^{\prime }}(f_{0\mathbf{k},\gamma }-f_{0\mathbf{k}-\mathbf{q},{\gamma }^{\prime }})\delta ({ϵ}_{\mathbf{k},\gamma }-{ϵ}_{\mathbf{k}-\mathbf{q},{\gamma }^{\prime }}-\hslash \omega )F_{{\mathbf{k}}^{\prime },{\mathbf{k}}^{\prime }+\mathbf{q}}^{\eta {\eta }^{\prime }}(f_{0{\mathbf{k}}^{\prime },\eta }-f_{0{\mathbf{k}}^{\prime }+\mathbf{q},{\eta }^{\prime }})\delta ({ϵ}_{{\mathbf{k}}^{\prime },\eta }-{ϵ}_{{\mathbf{k}}^{\prime }+\mathbf{q},{\eta }^{\prime }}+\hslash \omega )\times$$

$$[({\mathbf{v}}_{\mathbf{k},\gamma }-{\mathbf{v}}_{\mathbf{k}-\mathbf{q},{\gamma }^{\prime }}+{\mathbf{v}}_{{\mathbf{k}}^{\prime },\eta }-{\mathbf{v}}_{{\mathbf{k}}^{\prime }+\mathbf{q},{\eta }^{\prime }}){\mathbf{p}}_n+\beta ({\mathbf{v}}_{\mathbf{k},\gamma }{\widetilde{ϵ}}_{\mathbf{k},\gamma }-{\mathbf{v}}_{\mathbf{k}-\mathbf{q},{\gamma }^{\prime }}{\widetilde{ϵ}}_{\mathbf{k}-\mathbf{q},{\gamma }^{\prime }}+{\mathbf{v}}_{{\mathbf{k}}^{\prime },\eta }{\widetilde{ϵ}}_{{\mathbf{k}}^{\prime },\eta }-{\mathbf{v}}_{{\mathbf{k}}^{\prime }+\mathbf{q},{\eta }^{\prime }}{\widetilde{ϵ}}_{{\mathbf{k}}^{\prime }+\mathbf{q},{\eta }^{\prime }})]{\mathbf{p}}_s,$$

$$\left(\begin{array}[]{c}\sum_{{\bf k},\gamma}{\bf v}_{{\bf k},\gamma}\tilde{I}_{{\bf k},\gamma}\\ \sum_{{\bf k},\gamma}\beta\tilde{\epsilon}_{{\bf k},\gamma}{\bf v}_{{\bf k},\gamma}\tilde{I}_{{\bf k},\gamma}\end{array}\right)=\tilde{\bf I}\cdot\left(\begin{array}[]{c}{\bf p}_{n}\\ {\bf p}_{s}\end{array}\right).$$

$$\left(\begin{array}[]{c}\sum_{{\bf k},\gamma}{\bf v}_{{\bf k},\gamma}\tilde{I}_{{\bf k},\gamma}\\ \sum_{{\bf k},\gamma}\beta\tilde{\epsilon}_{{\bf k},\gamma}{\bf v}_{{\bf k},\gamma}\tilde{I}_{{\bf k},\gamma}\end{array}\right)=\frac{\pi\beta}{2}\sum_{{\bf q}}\int_{-\infty}^{\infty}d\omega\frac{|V({\bf q})|^{2}}{\sinh^{2}(\beta\hbar\omega/2)}\left(\begin{array}[]{cc}G_{11}&G_{12}\\ G_{21}&G_{22}\end{array}\right)\cdot\left(\begin{array}[]{c}{\bf p}_{n}\\ {\bf p}_{s}\end{array}\right),$$

$$G_{ij}=\left(\sum _{\mathbf{k},\gamma ,{\gamma }^{\prime }}{u}_{\mathbf{k},\gamma }^i{F}_{\mathbf{k},\mathbf{k}-\mathbf{q}}^{\gamma {\gamma }^{\prime }}(f_{0\mathbf{k},\gamma }-f_{0\mathbf{k}-\mathbf{q},{\gamma }^{\prime }})\delta ({ϵ}_{\mathbf{k},\gamma }-{ϵ}_{\mathbf{k}-\mathbf{q},{\gamma }^{\prime }}-\hslash \omega )\right)\times$$

$$\left(\sum _{{\mathbf{k}}^{\prime },\eta ,{\eta }^{\prime }}(u_{{\mathbf{k}}^{\prime },\eta }^j-u_{{\mathbf{k}}^{\prime }+\mathbf{q},{\eta }^{\prime }}^j)F_{{\mathbf{k}}^{\prime },{\mathbf{k}}^{\prime }+\mathbf{q}}^{\eta {\eta }^{\prime }}(f_{0{\mathbf{k}}^{\prime },\eta }-f_{0{\mathbf{k}}^{\prime }+\mathbf{q},{\eta }^{\prime }})\delta ({ϵ}_{{\mathbf{k}}^{\prime },\eta }-{ϵ}_{{\mathbf{k}}^{\prime }+\mathbf{q},{\eta }^{\prime }}+\hslash \omega )\right)+$$

$$\left(\sum _{\mathbf{k},\gamma ,{\gamma }^{\prime }}{u}_{\mathbf{k},\gamma }^i(u_{\mathbf{k},\gamma }^j-u_{\mathbf{k}-\mathbf{q},{\gamma }^{\prime }}^j)F_{\mathbf{k},\mathbf{k}-\mathbf{q}}^{\gamma {\gamma }^{\prime }}(f_{0\mathbf{k},\gamma }-f_{0\mathbf{k}-\mathbf{q},{\gamma }^{\prime }})\delta ({ϵ}_{\mathbf{k},\gamma }-{ϵ}_{\mathbf{k}-\mathbf{q},{\gamma }^{\prime }}-\hslash \omega )\right)\times$$

$$\left(\sum _{{\mathbf{k}}^{\prime },\eta ,{\eta }^{\prime }}{F}_{{\mathbf{k}}^{\prime },{\mathbf{k}}^{\prime }+\mathbf{q}}^{\eta {\eta }^{\prime }}(f_{0{\mathbf{k}}^{\prime },\eta }-f_{0{\mathbf{k}}^{\prime }+\mathbf{q},{\eta }^{\prime }})\delta ({ϵ}_{{\mathbf{k}}^{\prime },\eta }-{ϵ}_{{\mathbf{k}}^{\prime }+\mathbf{q},{\eta }^{\prime }}+\hslash \omega )\right).$$

$$\sum _{\mathbf{k},\gamma ,{\gamma }^{\prime }}{u}_{\mathbf{k},\gamma }^i{F}_{\mathbf{k},\mathbf{k}-\mathbf{q}}^{\gamma {\gamma }^{\prime }}(f_{\mathbf{k},\gamma }^0-f_{\mathbf{k}-\mathbf{q},{\gamma }^{\prime }}^0)\delta ({ϵ}_{\mathbf{k},\gamma }-{ϵ}_{\mathbf{k}-\mathbf{q},{\gamma }^{\prime }}-\hslash \omega )=$$

$$\frac{1}{2}\sum _{\mathbf{k},\gamma ,{\gamma }^{\prime }}{F}_{\mathbf{k},\mathbf{k}+\mathbf{q}}^{\gamma {\gamma }^{\prime }}(u_{\mathbf{k},\gamma }^i-u_{\mathbf{k}+\mathbf{q},{\gamma }^{\prime }}^i)(f_{\mathbf{k},\gamma }^0-f_{\mathbf{k}+\mathbf{q},{\gamma }^{\prime }}^0)\delta ({ϵ}_{\mathbf{k},\gamma }-{ϵ}_{\mathbf{k}+\mathbf{q},{\gamma }^{\prime }}+\hslash \omega )$$

$$\sum _{\mathbf{k},\gamma ,{\gamma }^{\prime }}{u}_{\mathbf{k},\gamma }^i(u_{\mathbf{k},\gamma }^j-u_{\mathbf{k}-\mathbf{q},{\gamma }^{\prime }}^j)F_{\mathbf{k},\mathbf{k}-\mathbf{q}}^{\gamma {\gamma }^{\prime }}(f_{\mathbf{k},\gamma }^0-f_{\mathbf{k}-\mathbf{q},{\gamma }^{\prime }}^0)\delta ({ϵ}_{\mathbf{k},\gamma }-{ϵ}_{\mathbf{k}-\mathbf{q},{\gamma }^{\prime }}-\hslash \omega )=$$

$$-\frac{1}{2}\sum _{\mathbf{k},\gamma ,{\gamma }^{\prime }}(u_{\mathbf{k},\gamma }^i-u_{\mathbf{k}+\mathbf{q},{\gamma }^{\prime }}^i)(u_{\mathbf{k},\gamma }^j-u_{\mathbf{k}+\mathbf{q},{\gamma }^{\prime }}^j)F_{\mathbf{k},\mathbf{k}+\mathbf{q}}^{\gamma {\gamma }^{\prime }}(f_{\mathbf{k},\gamma }^0-f_{\mathbf{k}+\mathbf{q},{\gamma }^{\prime }}^0)\delta ({ϵ}_{\mathbf{k},\gamma }-{ϵ}_{\mathbf{k}+\mathbf{q},{\gamma }^{\prime }}+\hslash \omega ).$$

$$G_{ij}=\frac{1}{2}\left(\sum _{\mathbf{k},\gamma ,{\gamma }^{\prime }}(u_{\mathbf{k},\gamma }^i-u_{\mathbf{k}+\mathbf{q},{\gamma }^{\prime }}^i)F_{\mathbf{k},\mathbf{k}+\mathbf{q}}^{\gamma {\gamma }^{\prime }}(f_{0\mathbf{k},\gamma }-f_{0\mathbf{k}+\mathbf{q},{\gamma }^{\prime }})\delta ({ϵ}_{\mathbf{k},\gamma }-{ϵ}_{\mathbf{k}+\mathbf{q},{\gamma }^{\prime }}+\hslash \omega )\right)\times$$

$$\left(\sum _{{\mathbf{k}}^{\prime },\eta ,{\eta }^{\prime }}(u_{{\mathbf{k}}^{\prime },\eta }^j-u_{{\mathbf{k}}^{\prime }+\mathbf{q},{\eta }^{\prime }}^j)F_{{\mathbf{k}}^{\prime },{\mathbf{k}}^{\prime }+\mathbf{q}}^{\eta {\eta }^{\prime }}(f_{0{\mathbf{k}}^{\prime },\eta }-f_{0{\mathbf{k}}^{\prime }+\mathbf{q},{\eta }^{\prime }})\delta ({ϵ}_{{\mathbf{k}}^{\prime },\eta }-{ϵ}_{{\mathbf{k}}^{\prime }+\mathbf{q},{\eta }^{\prime }}+\hslash \omega )\right)-$$

$$\frac{1}{2}\left(\sum _{\mathbf{k},\gamma ,{\gamma }^{\prime }}(u_{\mathbf{k},\gamma }^i-u_{\mathbf{k}+\mathbf{q},{\gamma }^{\prime }}^i)(u_{\mathbf{k},\gamma }^j-u_{\mathbf{k}+\mathbf{q},{\gamma }^{\prime }}^j)F_{\mathbf{k},\mathbf{k}+\mathbf{q}}^{\gamma {\gamma }^{\prime }}(f_{\mathbf{k},\gamma }^0-f_{\mathbf{k}+\mathbf{q},{\gamma }^{\prime }}^0)\delta ({ϵ}_{\mathbf{k},\gamma }-{ϵ}_{\mathbf{k}+\mathbf{q},{\gamma }^{\prime }}+\hslash \omega )\right)\times$$

$$\left(\sum _{{\mathbf{k}}^{\prime },\eta ,{\eta }^{\prime }}{F}_{{\mathbf{k}}^{\prime },{\mathbf{k}}^{\prime }+\mathbf{q}}^{\eta {\eta }^{\prime }}(f_{0{\mathbf{k}}^{\prime },\eta }-f_{0{\mathbf{k}}^{\prime }+\mathbf{q},{\eta }^{\prime }})\delta ({ϵ}_{{\mathbf{k}}^{\prime },\eta }-{ϵ}_{{\mathbf{k}}^{\prime }+\mathbf{q},{\eta }^{\prime }}+\hslash \omega )\right).$$

$$\tilde{\bf I}=-\frac{\beta}{4\pi}\sum_{{\bf q}}\int_{-\infty}^{\infty}d\omega\frac{|V(q)|^{2}}{\sinh^{2}(\beta\hbar\omega/2)}\left(\begin{array}[]{cc}(\Im\Pi_{10})^{2}-\Im\Pi_{00}\Im\Pi_{20}&\Im\Pi_{10}\Im\Pi_{01}-\Im\Pi_{00}\Im\Pi_{11}\\ \Im\Pi_{10}\Im\Pi_{01}-\Im\Pi_{00}\Im\Pi_{11}&\leavevmode\nobreak\ \leavevmode\nobreak\ (\Im\Pi_{01})^{2}-\Im\Pi_{00}\Im\Pi_{02}\\ \end{array}\right)\,$$

$${\Pi }_{nm}=4\sum _{\gamma ,{\gamma }^{\prime }}\sum _{\mathbf{k}}\frac{F_{\mathbf{k},\mathbf{k}+\mathbf{q}}^{\gamma {\gamma }^{\prime }}(v_{\mathbf{k},\gamma }-v_{\mathbf{k}+\mathbf{q},{\gamma }^{\prime }}{)}^n({\widetilde{ϵ}}_{k,\gamma }{v}_{\mathbf{k},\gamma }-{\widetilde{ϵ}}_{|k+q|,{\gamma }^{\prime }}{v}_{\mathbf{k}+\mathbf{q},{\gamma }^{\prime }}{)}^m(f_{0\mathbf{k},\gamma }-f_{0\mathbf{k}+\mathbf{q},{\gamma }^{\prime }})}{{ϵ}_{\mathbf{k},\gamma }-{ϵ}_{\mathbf{k}+\mathbf{q},{\gamma }^{\prime }}+\hslash \omega +i{0}^{+}}.$$

$${\mathbf{I}}_{D\mathbf{k},\gamma }=\frac{2\pi g}{\hslash }\sum _{{\gamma }^{\prime }}\sum _{{\mathbf{k}}^{\prime }}|U_{\mathbf{k}-{\mathbf{k}}^{\prime }}{|}^2{F}_{\mathbf{k},{\mathbf{k}}^{\prime }}^{\gamma {\gamma }^{\prime }}(f_{\mathbf{k},\gamma }-f_{{\mathbf{k}}^{\prime },{\gamma }^{\prime }})\delta ({ϵ}_{\mathbf{k},\gamma }-{ϵ}_{{\mathbf{k}}^{\prime },{\gamma }^{\prime }})$$

Full text

*Introduction – * The electric and thermal transport properties of graphene-based devices are a topic of great interest. Even setting aside their great potential for real-world applications, these systems have already offered unprecedented opportunities to study new modalities of transport, in which hydrodynamic flow patterns, governed by global conservation laws, supersede the conventional diffusive dynamics of individual carriers Narozhny et al. (2017); Principi et al. (2016); Narozhny et al. (2015); Briskot et al. (2015); Lucas and Fong (2018). The simultaneous presence of carriers of opposite polarities – electrons and holes – whose density can be tuned by chemical doping, electrostatic gating, or simply by changing the temperature, creates a rich scenario of transport behaviors Bandurin et al. (2016); Ghahari et al. (2016); Krishna Kumar et al. (2017); Gurzhi (1963). A particularly interesting one has recently been observed in a single layer of ultra-clean graphene near the charge neutrality point (CNP), where the chemical potential $\mu=0$. The system is a zero-gap semiconductor with linearly dispersing conduction and valence bands and equal numbers of electrons and holes arising from thermal fluctuations at finite temperature. Because $\mu=0$ the thermal (entropy) current coincides with the energy current, and because of the linear dispersion the latter coincides with the total momentum density, which is a constant of the motion as long as impurities, lattice vibrations and umklapp effects are negligible. Thus we have an interesting situation in which the thermal resistivity vanishes while the electric resistivity remains finite because electrons and holes, moving in opposite directions under the action of an electric field, exert mutual friction on each other Zarenia et al. (2018); Fritz et al. (2008); Müller et al. (2008); Svintsov et al. (2012). The result is that the Wiedemann-Franz (WF) ratio between the electric resistivity and the thermal resistivity is enhanced well above the standard value of $L_{0}\equiv\pi^{2}(k_{\rm B}/e)^{2}/3$ (the so-called Wiedemann-Franz law) – an effect that is indeed observed experimentally Crossno et al. (2016), but only in a narrow window of doping densities around the CNP – a window whose width shrinks to zero as the system is made less and less disordered. It should also be noted that this behavior is diametrically opposite to what one expects and observes in heavily doped graphene: in that case, only one polarity of carriers contributes to both electric and thermal transport and the electric resistivity plummets (barring electron-impurity and electron-phonon scattering) while the thermal resistivity rises to a finite value: in this case, the WF ratio drops below $L_{0}$ Principi and Vignale (2015); Lucas and Das Sarma (2018); Lucas and Hartnoll (2018); Xie and Foster (2016).

Motivated by these interesting findings, in this paper we investigate thermoelectric transport in a more complex system, AB-stacked bilayer graphene (BLG), which is a zero-gap semiconductor with electron-hole symmetric parabolic bands touching at the Dirac point McCann and Koshino (2013a). BLG is a quite distinct system from its single layer with the possibility of having a gate-induced tunable band gap Zhang et al. (2009). We ask in particular, whether a large violation of the WF law will still be present in the double-layer system near CNP. At first glance, a major qualitative difference exists between the two systems, because neither the energy current nor the particle current are conserved in the bilayer. There seems to be no reason why the thermal resistivity would vanish in a bilayer at CNP when momentum-non-conserving processes are negligible. Indeed, a naive calculation, based on the textbook theory of thermoelectric transport, would lead precisely to this conclusion: that the thermo-electric transport coefficients are free of singularities and qualitatively similar to what was obtained in single-layer graphene after the inclusion of disorder.

One of the main purposes of this paper is to show that the naive conclusion is, in fact, incorrect. After introducing a more careful treatment of the Boltzmann equation, which includes a conserved mode in which electrons and holes travel in the same direction, we are able to show that the transport coefficients remain singular as long as total momentum is conserved: that is to say, we find that, just as in single-layer graphene Zarenia et al. (2018), the electrical resistivity jumps from zero to a finite universal value (controlled by the strength of the Coulomb interaction) at the CNP, the thermal resistivity tends to zero, and the Seebeck coefficient diverges. This happens because the conserved mode – electrons and holes traveling in the same direction – strongly overlaps the energy current mode and provides a dissipation-free channel of energy transport at CNP. The inclusion of disorder, even in the smallest amount, “cures” the singularities and creates a region of “disorder-enabled hydrodynamics” in the immediate vicinity of the CNP. Thus, the second purpose of this paper is to determine the qualitative behavior of the thermo-electric transport coefficients of bilayer graphene in this regime. We find that the WF ratio, plotted as a function of doping density, follows a squared Lorentzian behavior, whose quarter-maximum occurs at doping density $n_{c}$ proportional to the strength of disorder, such that the WF ratio increases with decreasing disorder for $n<n_{c}$ and decreases with decreasing disorder for $n>n_{c}$. At the same time, the Seebeck coefficient exhibits an interesting non-monotonic behavior, vanishishing at CNP and peaking in absolute value at $n\sim n_{c}$. This behavior, admittedly very similar to what has been predicted and observed in single-layer graphene, should be promptly comparable with the results of experimental measurements, as soon as they become available.

Within the framework of quasi-classical transport theory Ashcroft and Mermin (2011), the state of the carriers is described by a non-equilibrium distribution function $f_{{\bf k},\gamma}$, where ${\bf k}$ is the Bloch wave vector and $\gamma=\pm 1$ is the band index. The deviation from equilibrium is $\delta f_{{\bf k},\gamma}=f_{{\bf k},\gamma}-f_{0{\bf k},\gamma}$, where $f_{0{\bf k},\gamma}$ is the equilibrium distribution function at chemical potential $\mu$ and temperature $T$. The quantities of interest are the electric current ${\bf j}_{e}$ and the thermal current ${\bf j}_{q}$, however, in order to homogenize the dimensions we will be working with the particle current ${\bf j}_{n}={\bf j}_{e}/(-e)$ ($-e$ is the charge of the electron) and the entropy current in units of the Boltzmann constant $k_{B}$, ${\bf j}_{s}=\beta{\bf j}_{q}$, where $\beta=(k_{B}T)^{-1}$. With this choice, the thermoelectric matrix [see below Eq. (2)] manifestly satisfies Onsager reciprocity. The currents are related to the non-equilibrium distribution function by Pines and Nozières (1966)

[TABLE]

where $\epsilon_{{\bf k},\gamma}=\gamma[\sqrt{(t/2)^{2}+(\hbar vk)^{2}}-t/2]$ and ${\bf v}_{{\bf k},\gamma}=2\gamma\hbar v^{2}{\bf k}/\sqrt{t^{2}+(2\hbar vk)^{2}}$ are, respectively, the energy and the velocity of band $\gamma=\pm$ at the wave vector ${\bf k}$, while $\tilde{\epsilon}_{{\bf k},\gamma}\equiv\epsilon_{{\bf k},\gamma}-\mu$. Here $v=10^{8}\leavevmode\nobreak\ {\rm cm/s}$ is the Fermi velocity and $t=0.4$ eV is the vertical interlayer hopping between the two layers McCann et al. (2007); McCann and Koshino (2013b). Our two-band approximation is justified up to $T\sim t/k_{\rm B}\sim 4,600\leavevmode\nobreak\ {\rm K}$.

The currents are connected to the electric field ${\bf E}$ and to the temperature gradient \mbox{\boldmath\nabla}T by the thermoelectric resistivity matrix, $\rho$, which we define as -\left(e{\bf E}\leavevmode\nobreak\ ,\leavevmode\nobreak\ k_{B}\mbox{\boldmath\nabla}T\right)^{\mathrm{T}}=\mbox{\boldmath\rho}\cdot\left({\bf j}_{n}\leavevmode\nobreak\ ,\leavevmode\nobreak\ {\bf j}_{s}\right)^{\mathrm{T}} (here “T” stands for the vector transposition). The elements of $\rho$ are expressed in terms of three transport coefficients: the reduced electric resistivity $\bar{\rho}_{el}$, i.e. the ordinary electric resistivity multiplied by $e^{2}$, the reduced thermal resistivity $\bar{\rho}_{th}$, i.e. the usual thermal resistivity multiplied by $k_{B}^{2}T$, and the dimensionless Seebeck coefficient $\bar{Q}$, i.e. the ordinary Seebeck coefficient expressed in units of $k_{B}/e$, in the following form (see Ref. Ashcroft and Mermin, 2011)

[TABLE]

with {\rm det}\mbox{\boldmath\rho}=\bar{\rho}_{el}\bar{\rho}_{th}. The (dimensionless) Wiedemann-Franz ratio is defined as WF\equiv\bar{\rho}_{el}/\bar{\rho}_{th}={\rm det}\mbox{\boldmath\rho}/\bar{\rho}_{th}^{2},. $WF=\pi^{2}/3$ when the Wiedemann-Franz law is satisfied in its standard form, e.g., for a parabolic band electron gas in the presence of quenched short-range disorder.

Two-mode Ansatz. The standard textbook calculation of thermoelectric coefficients Ashcroft and Mermin (2011) starts with the introduction of a 2-parameter Ansatz $\delta f_{{\bf k},\gamma}=f^{\prime}_{0{\bf k},\gamma}{\bf v}_{{\bf k},\gamma}\cdot({\bf p}_{n}+\beta\tilde{\epsilon}_{{\bf k},\gamma}{\bf p}_{s})$, where ${\bf v}_{{\bf k},\gamma}$ and $\beta\tilde{\epsilon}_{{\bf k},\gamma}{\bf v}_{{\bf k},\gamma}$ are the “modes” used to expand $\delta f_{{\bf k},\gamma}$. The two parameters ${\bf p}_{n}$ and ${\bf p}_{s}$ correspond to shifts of the particle momentum associated with the ${\bf j}_{n}$ and ${\bf j}_{s}$, respectively. Finally, the factor $f^{\prime}_{0{\bf k},\gamma}$, which denotes the derivative of the Fermi distribution with respect to energy, accounts for the fact that only electrons around the Fermi surface are mobile.

We stress that the choice of the Ansatz is the delicate point of the entire calculation. (The rest of the section is completely general and valid independently of such choice.) The two-parameter Ansatz above is, for BLG, incomplete and leads to the wrong results for $\bar{\rho}_{el}$ and $\bar{\rho}_{th}$. In the following section we will amend it by including a third parameter, corresponding to the current conserved by the collision integral (i.e. the momentum density).

Using Eq. (1) we obtain $\left({\bf j}_{n}\leavevmode\nobreak\ ,\leavevmode\nobreak\ {\bf j}_{s}\right)^{\mathrm{T}}={\bf D}\cdot\left({\bf p}_{n}\leavevmode\nobreak\ ,\leavevmode\nobreak\ {\bf p}_{s}\right)^{\mathrm{T}}$, where ${\bf D}$ is the $2\times 2$ matrix of “Drude weights”

[TABLE]

where ${\bf u}^{i}_{{\bf k},\gamma}=({\bf v}_{{\bf k},\gamma},\beta\tilde{\epsilon}_{{\bf k},\gamma}{\bf v}_{{\bf k},\gamma})_{i}$ and $i,j=1,2$. $D_{ij}$ quantifies the “overlap” between the modes ${\bf u}^{i}$ and ${\bf u}^{j}$ (it can in fact be interpreted as a scalar product in the mode space). To determine ${\bf p}_{n}$ and ${\bf p}_{s}$, we substitute the Ansatz for $\delta f_{{\bf k},\gamma}$ into the Boltzmann equation for the steady state response in the presence of fields ${\bf E}$ and \mbox{\boldmath\nabla}T,

[TABLE]

Here $I_{{\bf k},\gamma}$ is the collision integral, which depends on the details of the microscopic scattering mechanism. Eq. (4) is projected over the same set of modes which are used to expand $\delta f_{{\bf k},\gamma}$, i.e. it is multiplied by one of the modes ${\bf u}^{i}_{{\bf k},\gamma}$, integrated over ${\bf k}$ and summed over bands. In this way, the differential Eq. (4) is transformed into an algebraic one, and is easily solved for ${\bf p}_{n}$ and ${\bf p}_{s}$. The key inputs are the moments of the collision integral, which to linear order in ${\bf p}_{n}$ and ${\bf p}_{s}$ are given by $\sum_{{\bf k},\gamma}{\bf u}^{i}_{{\bf k},\gamma}I_{{\bf k},\gamma}=\tilde{I}_{ij}({\bf p}_{n}\leavevmode\nobreak\ ,\leavevmode\nobreak\ {\bf p}_{s})_{j}$. Hereafter summation of repeated latin indices is understood. Such equation defines the $2\times 2$ matrix $\tilde{\bf I}$, which we refer to as “collision kernel”, and whose matrix elements are the ${\tilde{I}}_{ij}$ ($i,j=1,2$). In the supplementary online material sup we make use of a standard approximation for the Coulomb collision integral (screened interaction plus Fermi golden rule) to find

[TABLE]

where the response functions $\Pi^{\alpha}({\bf q},\omega)$ are defined as

[TABLE]

Here $M^{0}_{{\bf k},{\bf k}+{\bf q},\gamma,\gamma^{\prime}}=1$, $M^{1}_{{\bf k},{\bf k}+{\bf q},\gamma,\gamma^{\prime}}=\hat{\bf x}\cdot({\bf u}^{i}_{{\bf k},\gamma}-{\bf u}^{i}_{{\bf k}+{\bf q},\gamma^{\prime}})$, and $M^{2}_{{\bf k},{\bf k}+{\bf q},\gamma,\gamma^{\prime}}=({\bf u}^{i}_{{\bf k},\gamma}-{\bf u}^{i}_{{\bf k}+{\bf q},\gamma^{\prime}})\cdot({\bf u}^{j}_{{\bf k},\gamma}-{\bf u}^{j}_{{\bf k}+{\bf q},\gamma^{\prime}})$. Furthermore, $V(q)=2\pi e^{2}/[\kappa(q+q_{\mathrm{TF}})]$ is the screened Coulomb interaction, $q_{\mathrm{TF}}=4e^{2}\varepsilon_{F}/(\kappa\hbar^{2}v^{2})$ is the Thomas-Fermi screening wave vector ($\varepsilon_{F}$ is the Fermi energy). The factor 4 accounts for the spin and valley degeneracy and $\kappa$ is the dielectric constant of the substrate (within our calculations, we set $\kappa=4$ as for an h-BN substrate). The form factors $F_{{\bf k},{\bf k}+{\bf q}}^{\gamma\gamma^{\prime}}=[1+\gamma\gamma^{\prime}\cos(2\theta_{{\bf k}}-2\theta_{{\bf k}+{\bf q}})]/2$, where $\theta_{\bf k}$ is the angle formed by the ${\bf k}$ vector with the $x$-axis, come from the overlap of the wave functions at wave vectors ${\bf k}$ and ${\bf k}+{\bf q}$.

Simple algebraic manipulations lead to the final expression for the two-mode thermoelectric resistivity matrix, \tilde{}\mbox{\boldmath\rho}={\bf D}^{-1}\cdot\tilde{\bf I}\cdot{\bf D}^{-1}. For future purposes we define the conductivity matrix \tilde{}\mbox{\boldmath\sigma}\equiv\tilde{}\mbox{\boldmath\rho}^{-1}={\bf D}\cdot\tilde{\bf I}^{-1}\cdot{\bf D}.

This procedure works remarkably well for both parabolic-band electron gases and monolayer graphene in the clean limit. Zarenia et al. (2018) In such limit, i.e. when the only collision mechanism is the electron-electron interaction, the thermoelectric transport is strongly influenced by the exact conservation of the momentum density. For example, when one of the two currents overlaps sufficiently with the momentum density [where the overlap is defined as in Eq. (3)], then it is also conserved. This in turn implies that either $\bar{\rho}_{el}$ or $\bar{\rho}_{th}$ vanishes. For parabolic-band electron gases and (massive or massless) Dirac systems, such result is readily obtained with a two-mode Ansatz. In such special cases, the overlap between the two currents and momentum (which coincides, in fact, with one of them at all times) is automatically taken into account.

It is then clear what fails when the two-mode Ansatz is applied to BLG: the overlap between the currents and the momentum density, the conserved mode of the collision integral, is not explicitly taken into account. Hence, all the coefficients of the matrix \tilde{}\mbox{\boldmath\rho} are non-zero, and so is its determinant. Therefore, all thermoelectric transport coefficients are finite, see Fig. 1(a). Such result is wrong as we proceed to show in the next section.

Three-mode Ansatz. To account for the conserved mode of the collision integral (the momentum density) we express the deviation of $f_{{\bf k},\gamma}$ from equilibrium as

[TABLE]

where ${\bf p}_{k}$ is the shift associated with the momentum density scaled with $t/\hbar v^{2}$, in such a way that it has the same dimension as the other modes.

A more general situation (away from the CNP) might require to include more modes, in order to achieve a better quantitative agreement with the true solution of the Boltzmann equation. However, as we will shown later, the Ansatz (7) is able to capture all the main qualitative features (vanishing of the thermal/charge conductivities, divergence of the Seebeck coefficient, etc.). In addition, there is no real locking imposed a priori between the electron and hole velocities. As such, the three modes of the anstaz describe different situations, in which particles and holes may co-propagate, counter-propagate, or any combination in between. The coefficients of the separate modes are allowed to vanish, describing situations in which one particular configuration is realised.

We now define the three-component mode vector as ${\bf u}^{j}_{{\bf k},\gamma}=(\hbar v^{2}{\bf k}/t,{\bf v}_{{\bf k},\gamma},\beta\tilde{\epsilon}_{{\bf k},\gamma}{\bf v}_{{\bf k},\gamma})_{j}$. Using Eq. (1) and (3), we get $({\bf j}_{n},{\bf j}_{s})_{i}=D_{ij}({\bf p}_{k},{\bf p}_{n},{\bf p}_{s})_{j}$, where $i=1,2$ and $j=0,1,2$,. Note that ${\bf D}$ is now a $2\times 3$ matrix. Similarly, for the Coulomb collision kernel we find

[TABLE]

where $j,j^{\prime}=0,1,2$. The $\tilde{I}_{ij}$ are defined in Eq. (A.12). The conductivity matrix is easily seen to be \mbox{\boldmath\sigma}={\bf D}\cdot{\bf I}^{-1}\cdot{\bf D}. Inverting it, in the limit of $\eta\to 0$, we obtain

[TABLE]

where $\tilde{\sigma}_{ii^{\prime}}=\lim_{\eta\to 0}\sigma_{ii^{\prime}}$ coincide with the conductivities obtained with the two-mode Ansatz. Note that, since {\rm det}\mbox{\boldmath\rho}=\bar{\rho}_{el}\bar{\rho}_{th}=0, either the electrical or thermal resistivity must necessarily be zero. At charge neutrality $D_{01}=0$. Therefore, ${\bf\rho}_{11}(n=0)=\tilde{\sigma}_{11}^{-1}$, while all the other components of $\rho$ are zero. Thus, at the CNP, the thermal resistivity vanishes and the electrical resistivity is finite. Away from the CNP we find

[TABLE]

Numerical results for $\bar{\rho}_{el}$ and $\bar{\rho}_{th}$ are shown in Fig. 1. Within the two-mode approximation [Fig. 1(a)] both the electrical and thermal resistivities are always finite. In contrast, with the three mode Ansatz [Fig. 1(b)], the electrical resistivity $\bar{\rho}_{el}(n\neq 0)=0$ and exhibits a discontinuity at the CNP in the absence of disorder (solid curves). Interestingly, both models yield the same intrinsic resistivity at $n=0$ [full red dots in Fig. 1(a) and (b)]. At this point the only carriers in the system are thermally excited electrons and holes, in the conduction and valence bands, respectively. The two types of carries drift in opposite directions under the action of an electric field. Because of the transfer of momentum (known as Coulomb drag Schütt et al. (2013)) between them, the resistivity becomes finite.

Disorder-enabled hydrodynamics. An infinitesimal amount of disorder, which breaks the exact conservation of momentum, regularizes the singular results of the clean limit. We assume a momentum-non-conserving kernel ${\bf I}_{D}$ (D for disorder) proportional to a dimensionless momentum relaxation rate $\lambda$, which we take to be $\ll 1$. The precise form of ${\bf I}_{D}$ is not important for our purposes. However, for the sake of illustration, we will later make use of a simple model of electrons and holes scattering against randomly distributed impurities of density $n_{d}$ with short-range potential $V_{0}\delta({\bf r})$: for this model $\lambda=n_{d}V_{0}^{2}/(\hbar v)^{2}$, as detailed in the supplementary online material sup .

The resistivity matrix is now the sum of two term \mbox{\boldmath\rho}=\mbox{\boldmath\rho}_{C}+\mbox{\boldmath\rho}_{D}, where \mbox{\boldmath\rho}_{C} is given by Eq. (9), whereas \mbox{\boldmath\rho}_{D}={\bf D}^{-1}\cdot{\bf I}_{D}\cdot{\bf D}^{-1}. The electrical and thermal resistivities now read Zarenia et al. (2018) $\bar{\rho}_{el}^{-1}\simeq\bar{\rho}_{el,C}^{-1}+\tilde{\rho}_{D}^{-1}$ and $\bar{\rho}_{th}=\bar{\rho}_{th,C}+\bar{\rho}_{th,D}$, where $\tilde{\rho}_{D}=\bar{\rho}_{th,D}(D_{02}/D_{01})^{2}$. In Fig. 1(b) we plot such $\bar{\rho}_{el}$ and $\bar{\rho}_{th}$ (dotted curves) and compare them with the results of the previous section (solid curves). The effect of disorder on the thermal resistivity is just a small shift. Conversely, it regularizes the electric resistivity. From the results above we obtain the WF ratio

[TABLE]

which is a square of a Lorentzian. This formula shows that, at the CNP, $WF\to 1/\Gamma^{2}$, i.e. it is greatly enhanced relative to its standard noninteracting value $\pi^{2}/3$ and diverges as the strength of disorder tends to zero.

In Figs. 2(a) and (b) we show numerical results for $WF$ and the Seebeck coefficient $\bar{Q}$, respectively. We define a crossover density $n_{c}$, below which the enhancement of the WF persists, from the condition $\tilde{\rho}_{D}=\bar{\rho}_{el,C}$ [the dots in Fig. 2(a) indicate its position]. We call such regime ($n<n_{c}$) “disorder-enabled hydrodynamics” Zarenia et al. (2018). In it, $WF$ remains much larger than $\pi^{2}/3$, in fact larger than $1/\Gamma^{2}\gg\pi^{2}/3$. (A more lenient crossover density could be defined Zarenia et al. (2018) as that at which $WF$ first drops below $\pi^{2}/3$.)

The disorder-regularized Seebeck coefficient [Figs. 2(b)] exhibits a large swing about the CNP and goes to zero at $n=0$, as expected from particle-hole symmetry. The swing region, in which the derivative of $\bar{Q}$ vs density reverses its sign, is yet another incarnation of the disorder-enabled hydrodynamic regime Zarenia et al. (2018). Its width is defined by the same condition $|n|<n_{c}$.

In Summary, we have presented the theory of thermoelectric transport in clean bilayer graphene. We have shown that the conventional semi-classical (Boltzmann) textbook approach to the calculation of thermoelectric coefficients Ashcroft and Mermin (2011), which works remarkably well for both parabolic-band electron gases and Dirac systems Zarenia et al. (2018), fails for such system. This is attributed to the fact that neither the particle current nor the energy current are conserved quantities in bilayer graphene. The correct results are found by explicitly including the current conserved by the electron-electron collision integral, i.e. the momentum density, in the “augmented” Ansatz for the non-equilibrium distribution function. Note that in previously studied examples the explicit addition of the conserved momentum mode was not required. This fortunate situation occurred because the momentum mode was automatically subsumed under either the particle current or the energy current mode. Bilayer graphene is the first system studied in this context in which the two-mode Ansatz fails and the momentum mode must be introduced explicitly.

We find that, at the charge neutrality point: (i) the thermal resistivity vanishes; (ii) the electric resistivity jumps from zero to a finite value; (iii) the Seebeck coefficient diverges. These singularities are cured by the inclusion of a small amount of disorder. Breaking the exact momentum conservation, it opens a “window” of doping densities $0<n<n_{c}$ (with $n_{c}$ tending to zero in the zero-disorder limit) in which the Wiedemann-Franz law is largely violated. There, the WF ratio greatly exceeds the standard value and the Seebeck coefficient exhibits a non-monotonic behavior (as a function of doping density). Such predictions can be tested in experiments in sufficiently clean samples of bilayer graphene.

Supplemental Material

Breakdown of the Wiedemann-Franz law in AB-stacked bilayer graphene

In this supplemental material, we present more details on our calculations for the electron-electron Coulomb (Sec. A) as well as our short-range model for the disorder collision moments (Sec. B).

A. Electron-electron collision moments

The electron-electron collision integral for the band $\gamma$ and at wave vector ${\bf k}$, $I_{{\bf k},\gamma}$ is given by Ashcroft and Mermin (2011); Pines and Nozières (1966),

[TABLE]

where the momentum conservation appears naturally when doing the second quantization of Coulomb interaction in ${\bf k}$-space, and the energy conservation stems from the Fermi golden rule. $f_{{\bf k},\gamma}=f_{0{\bf k},\gamma}+\delta f_{{\bf k},\gamma}$ is the non-equilibrium distribution function and $W(q)=(2\pi/\hbar)|V(q)|^{2}$ defines the collision probability where $V({\bf q})$ is the statistic screened Coulomb interaction. Having the two-mode Ansatz

[TABLE]

and using

[TABLE]

and

[TABLE]

we obtain

[TABLE]

where $\{\gamma,\gamma^{\prime},\eta,\eta^{\prime}\}=\pm 1$ denote the band index and $F_{{\bf k},{\bf k}\pm{\bf q}}^{\gamma\gamma^{\prime}}$ is the form factor comes from the overlap of the wave functions at vector ${\bf k}$ and ${\bf k}\pm{\bf q}$. Now the moments of the electron-electron collision integrals are calculated using

[TABLE]

Using Eq. (A. Electron-electron collision moments) we write,

[TABLE]

where the matrix elements $G_{ij}({\bf q},\omega,T)$ are

[TABLE]

As defined in the main text, the two-component vector ${\bf u}^{i}_{{\bf k},\gamma}=({\bf v}_{{\bf k},\gamma},\beta\tilde{\epsilon}_{{\bf k},\gamma}{\bf v}_{{\bf k},\gamma})_{i}$ ($i=1,2$) . By simple algebraic manipulations, one can show that the first and the third lines in Eq. (A. Electron-electron collision moments) are equivalent to

[TABLE]

and

[TABLE]

Inserting Eqs. (A. Electron-electron collision moments) and (A. Electron-electron collision moments) into Eq. (A. Electron-electron collision moments) leads to

[TABLE]

Using this equation, $\tilde{\bf I}$ becomes,

[TABLE]

where the $\delta$-functions are written in terms of the imaginary part of the response functions $\Pi_{nm}(q,\omega)$ ($n,m=0,1,2$)

[TABLE]

Equations (A.12) and (A.13) are respectively represented in a more compact forms of Eqs. (5) and (6) in the main text.

I B. Disorder collision moments

The non-momentum-conserving disorder collision integral is given by,

[TABLE]

where the factor $g=4$ accounts for the spin and valley degeneracy. Considering the short-range disorder characterized by an effective strength of $|U_{{\bf k}-\boldsymbol{k^{\prime}}}|\approx(n_{d}V_{0}^{2})^{1/2}$ ($n_{d}$ is the disorder density and $V_{0}\delta(\bf r)$ is the disorder potential), the linearized collision integral becomes

[TABLE]

in which we set $F_{{\bf k},{\bf k}^{\prime}}^{\gamma\gamma^{\prime}}=1$ for the sake of simplicity and without loss of qualitative behavior of disorder. Using $\delta(\epsilon_{{\bf k},\gamma}-\epsilon_{{\bf k}^{\prime},\gamma^{\prime}})=\delta({\bf k}-{\bf k}^{\prime})/\hbar|{\bf v}_{{\bf k},\gamma}|$, we easily construct the moments of the disorder collision integral

[TABLE]

where,

[TABLE]

Taking $k_{F}$, $v$, and $\beta^{-1}$, respectively as units of wave vector, velocity and energy, ${\bf I}_{D}$ can be further simplified to,

[TABLE]

where the dimensionless disorder strength $\lambda=n_{d}V_{0}^{2}/(\hbar v)^{2}$. We solve the integrals in Eq. (B.4) numerically for the full energy bands of bilayer graphene.