Quantum Boltzmann equation for bilayer graphene

TL;DR

This paper derives a quantum Boltzmann equation for bilayer graphene near charge neutrality, incorporating various scattering effects, and computes transport coefficients, providing a comprehensive framework for understanding its non-equilibrium transport properties.

Contribution

It develops a quantum Boltzmann formalism for bilayer graphene including phonons, disorder, and magnetic fields, and introduces a simplified hydrodynamic model for this system.

Findings

Calculated electrical and thermal conductivities near charge neutrality.

Extended formalism to include magnetic field effects.

Derived a hydrodynamic model matching numerical solutions.

Abstract

A-B stacked bilayer graphene has massive electron and hole-like excitations with zero gap in the nearest-neighbor hopping approximation. In equilibrium, the quasiparticle occupation approximately follows the usual Fermi-Dirac distribution. In this paper we consider perturbing this equilibrium distribution so as to determine DC transport coefficients near charge neutrality. We consider the regime $\beta |\mu| \lesssim 1$ (with $\beta$ the inverse temperature and $\mu$ the chemical potential) where there is not a well formed Fermi surface. Starting from the Kadanoff-Baym equations, we obtain the quantum Boltzmann equation of the electron and hole distribution functions when the system is weakly perturbed out of equilibrium. The effect of phonons, disorder, and boundary scattering for finite sized systems are incorporated through a generalized collision integral. The transport coefficients, including the electrical and thermal conductivity, thermopower, and shear viscosity, are calculated in the linear response regime. We also extend the formalism to include an external magnetic field. We present results from numerical solutions of the quantum Boltzmann equation. Finally, we derive a simplified two-fluid hydrodynamic model appropriate for this system, which reproduces the salient results of the full numerical calculations.

Figures

Equations

$$H_{\xi }=\xi \left(\begin{array}{cccc}0& v_3\pi & 0& v_F{\pi }^{†}\\ v_3{\pi }^{†}& 0& v_F\pi & 0\\ 0& v_F{\pi }^{†}& 0& \xi {\gamma }_1\\ v_F\pi & 0& \xi {\gamma }_1& 0\end{array}\right)-e{A}_{0}I,$$

$${\psi }_K(\mathbf{x})$$

$$p_i=-i\hslash {\partial }_i-e{A}_i$$

$$\pi =p_x+i{p}_y,{\pi }^{†}=p_x-i{p}_y.$$

$$H_K=-\frac{1}{2m}\left(\begin{array}{cc}0& ({\pi }^{†}{)}^2\\ {\pi }^2& 0\end{array}\right),{\psi }_K=\left(\begin{array}{c}{\phi }_{A_1}(\mathbf{x})\\ {\phi }_{B_2}(\mathbf{x})\end{array}\right),$$

$${\psi }_K^{\lambda }=\frac{1}{\sqrt{2}}\left(\begin{array}{c}-\lambda e^{-2i{\theta }_{\mathbf{k}}}\\ 1\end{array}\right),{ϵ}_{\mathbf{k}}^{\lambda }=\lambda \frac{k^2}{2m},$$

$${\mathcal{L}}_{\xi }=i{\hat{\Psi }}_{\xi }^{†}\overleftrightarrow{{\partial }_t}{\hat{\Psi }}_{\xi }-{\hat{\Psi }}_{\xi }^{†}{H}_{\xi }{\hat{\Psi }}_{\xi },$$

$${\hat{\Psi }}_K(\mathbf{x})=\int \frac{d^{2}k}{(2\pi {)}^2}\left(\begin{array}{c}{c}_{K;A_1}(\mathbf{k})\\ c_{K;B_2}(\mathbf{k})\\ c_{K;A_2}(\mathbf{k})\\ c_{K;B_1}(\mathbf{k})\end{array}\right)e^{i\mathbf{kx}},$$

$${\mathcal{S}}_{\xi }=\int d^{3}x{\mathcal{L}}_{\xi }.$$

$$\rho (\mathbf{x})=\sum _{\xi }{\rho }_{\xi }(\mathbf{x})=\sum _{\xi ,s}\frac{1}{e}\frac{\delta {\mathcal{S}}_{\xi }}{\delta A_0}(\mathbf{x})=\sum _{\xi ,s}{\hat{\Psi }}_{\xi }^{†}(\mathbf{x}){\hat{\Psi }}_{\xi }(\mathbf{x}).$$

$$c_{\pm K}(\mathbf{k})=\frac{1}{\sqrt{2}}\left(\mp e^{-2i{\theta }_{\mathbf{k}}}{c}_{K;A_1}(\mathbf{k})+c_{K;B_2}(\mathbf{k})\right),$$

$$\begin{gathered} {\rho }^{\text{eff}}(\mathbf{q})=\sum _f\int \frac{d^{2}k}{(2\pi {)}^2}\frac{1}{2} \\ \times \sum _{\lambda ,{\lambda }^{\prime }}{c}_{\lambda f}^{†}(\mathbf{k}-\mathbf{q})c_{{\lambda }^{\prime }f}(\mathbf{k})\left(1+\lambda {\lambda }^{\prime }{e}^{-2i({\theta }_{\mathbf{k}-\mathbf{q}}-{\theta }_{\mathbf{k}})}\right). \end{gathered}$$

$${\hat{V}}_C=\frac{1}{2}\int \frac{d^{2}q}{(2\pi {)}^2}{V}_C(\mathbf{q}){\rho }^{\text{eff}}(-\mathbf{q}){\rho }^{\text{eff}}(\mathbf{q}).$$

$$V(q)=\frac{2\pi \alpha }{q}$$

$$V_C(q,\omega )=\frac{V(q,\omega )}{1-{\Pi }^0(q,\omega )V(q,\omega )},$$

$${\Pi }^0(\mathbf{q},0)\approx -\frac{m{N}_f}{2\pi }(1+\frac{\beta q^2}{12m}).$$

$$V_C(q)=\frac{2\pi \alpha }{q+q_{TF}(q)}$$

$$q_{TF}(q)=-{\Pi }^0(\mathbf{q},0)2\pi \alpha .$$

$$V_C(q)=\frac{2\pi \alpha }{q_{TF}(q)}.$$

$$\tilde{\mathcal{V}}(\mathbf{q},0)=\frac{2\pi \alpha }{\sqrt{g}}\frac{1}{|\mathbf{q}|+q_{TF}(|\mathbf{q}|)},$$

$$J_{\xi }^i(\mathbf{x})=\frac{\delta {\mathcal{S}}_{\xi }}{\delta A_i(\mathbf{x})}=-v_{F}e\sum _s\xi {\hat{\Psi }}_{\xi }^{†}(\mathbf{x})\left(\begin{array}{cc}0& ({\sigma }^i{)}^{†}\\ ({\sigma }^i{)}^{†}& 0\end{array}\right){\hat{\Psi }}_{\xi }(\mathbf{x})$$

$$\mathbf{J}(\mathbf{x})=\sum _{\xi }{\mathbf{J}}_{\xi }^{\text{I}}(\mathbf{x})+{\mathbf{J}}_{\xi }^{\text{II}}(\mathbf{x})$$

$${\mathbf{J}}^{\text{I}}(\mathbf{q}=0)=\frac{e}{m}\sum _f\sum _{\lambda }\int \frac{d^{2}k}{(2\pi {)}^2}\lambda \mathbf{k}{c}_{\lambda f}^{†}(\mathbf{k})c_{\lambda f}(\mathbf{k}),$$

$$\begin{gathered} {\mathbf{J}}^{\text{II}}(\mathbf{q}=0)=-i\frac{e}{m}\sum _f\int \frac{d^{2}k}{(2\pi {)}^2}(\hat{z}\times \mathbf{k}) \\ \times \left[c_{+f}^{†}(\mathbf{k})c_{-f}(\mathbf{k})-c_{-f}^{†}(\mathbf{k})c_{+f}(\mathbf{k})\right]. \end{gathered}$$

$${\mathbf{J}}^Q(\mathbf{q}=0)={\mathbf{J}}^E(\mathbf{q}=0)-\frac{\mu }{e}\mathbf{J}(\mathbf{q}=0),$$

$$\partial \cdot {\mathbf{J}}^E(\mathbf{x},t)+\dot{\mathcal{E}}(\mathbf{x},t)=0,$$

$${\mathcal{E}}_{kin}(\mathbf{x})=\sum _{\xi }{\mathbf{H}}_{\xi }(\mathbf{x})=\sum _{\xi }{\hat{\Psi }}_{\xi }^{†}(\mathbf{x}){\overleftrightarrow{H}}_{\xi }{\hat{\Psi }}_{\xi }(\mathbf{x})$$

$${\mathbf{H}}_{\xi }(\mathbf{q})$$

$$=\int \frac{d^{2}k}{(2\pi {)}^2}{\hat{\Psi }}_{\xi }^{†}(\mathbf{k}-\mathbf{q}){\hat{\mathbf{H}}}_{\xi }(\mathbf{k},\mathbf{q}){\hat{\Psi }}_{\xi }(\mathbf{k}),$$

$$\begin{gathered} {\hat{\mathbf{H}}}_{\xi }(\mathbf{k},\mathbf{q})= \\ \xi \left(\begin{array}{cccc}0& 0& 0& v_F(\bar{k}-\frac{1}{2}\bar{q})\\ 0& 0& v_F(k-\frac{1}{2}q)& 0\\ 0& v_F(\bar{k}-\frac{1}{2}\bar{q})& 0& \xi {\gamma }_1\\ v_F(k-\frac{1}{2}q)& 0& \xi {\gamma }_1& 0\end{array}\right), \end{gathered}$$

$$\mathbf{q}\cdot {\mathbf{J}}_{kin}^E(\mathbf{q})=[{\mathcal{E}}_{kin}(\mathbf{q}),\mathcal{H}],$$

Full text

I Introduction

Graphene is a remarkable material that has generated enormous interest in both the theoretical and experimental community since its discovery in 2004 [1]. While there are many reasons for interest in this unusual material, one particularly exciting feature is that it is possible to produce graphene samples with extreme purity and thus observe an electronic regime that had not previously been explored — the hydrodynamic regime [2, 3, 4, 5, 6]. In this regime the shortest scattering time is that for electron-electron collisions, and collisions with phonons and impurities are subdominant. As a result, semi-classical kinetic theory reduces to a form of quantum hydrodynamics [7].

In general, quantum Hydrodynamics (QH) describes the dynamics of systems that vary slowly in space and time. The foundation of QH is the ability to obtain a set of conservation laws for the electron liquid. These conservation laws are derived from the quantum Boltzmann equation (QBE), which is the equation of motion for the electron fluid’s phase space distribution function [8, 9]. The QH approach has been extremely successful in studying electron plasmas, fractional quantum Hall fluids, [10, 11, 12], and now also the electron fluid of graphene [13, 2]. In the current work we will derive QH equations from the QBE for the case of bilayer graphene (BLG) near charge neutrality (CN) where there is no well defined Fermi surface. We note, however, that our kinetic theory formalism is more general than QH.

Whereas QBE and QH have been studied extensively in monolayer graphene [13, 14, 15, 16], they have not been as well studied in the case of bilayer graphene. There are, however, several reasons why the BLG case should be of substantial interest. Firstly, the band structure of BLG is fundamentally different from that of monolayer graphene. In the nearest-neighbor hopping approximation, A-B stacked bilayer graphene has quadratic bands of electron and hole-like excitations at low energies [17] which touch at zero energy. This interesting band structure, which was confirmed experimentally in Ref [18], provides the unique quantum transport properties of BLG [19].

A second reason that BLG is now of interest is due to recent experimental advances that have allowed measurements with unprecedented precision — in particular [18] reports measurements of the electrical conductivity of BLG. These advances have been possible due to the development of suspended BLG devices [18, 20, 21, 22, 23, 24, 25, 26, 27]. As with monolayer graphene, BLG on a substrate suffers an inhomogeneous potential, which can lead to charge-puddle physics and superlattice effects. Suspended samples, in comparison, are far cleaner. The current limitation on suspended samples is on the size of these devices, however recent BLG devices have achieved sizes longer than the disorder scattering length [28].

Further, due to the low impurity scattering rates in clean samples, there has also been recent interest in studying the viscosity of the electron fluid in materials such as graphene [5, 29]. Some signatures of electron viscosity, such as negative non-local resistance, have already been measured experimentally in monolayer graphene[30]. The extension of these experiments to BLG seems natural.

Finally, measurements of the thermal conductivity of suspended single-layer graphene have been performed [31] and it may be possible to extend this study to BLG as well.

In this paper, we develop the QBE formalism for calculating transport properties of BLG which can in principle be compared against existing and future transport experiments. The analogous formalism for the electrical conductivity of single-layer graphene was worked out in [13]. This work was later extended to study Coulomb drag between two monolayers of graphene [16] and BLG exactly at charge neutrality (CN) [32]. We extend this work away from CN. To do so, we start with the derivation of the quantum transport operators including the charge current, the heat current and the stress tensor in terms of electron and hole fields. We include the Zitterbewegung contribution to these operators. Taking the expectation value of these operators in the DC limit, we obtain the main results for the electrical current (50), heat current (51) and stress tensor (52). We then calculate the collision integral, which in general will have four contributions. A central result of our work is the collision integral (69), which contains contributions from the Coulomb interactions (70), impurity scattering (76), scattering off the boundary (78) and phonon scattering (81).

Once we have derived the QBE formalism, we can study the linear response of the system to perturbations. In particular, we study the behavior under an applied external electric field, thermal gradient, and straining motion in order to calculate the electrical conductivity, thermal conductivity, and viscosity respectively. Finally, in anticipation of future experiments, we also consider the behavior of the transport properties in an applied magnetic field which leads to, e.g., nonzero Hall-conductivity. The inclusion of a magnetic field is possible for both electrical conductivity and thermal conductivity calculations.

Since the dominant scattering mechanism is electron-electron, and hole-hole, we should be able to represent the transport with a two-fluid hydrodynamics — where the electron fluid and the hole fluid are individually in a thermal equilibrium, each having a well defined temperature, chemical potential, and velocity. As such we use the QBE to derive a two-fluid model which describes the evolution of the mean fluid velocities of the electron and hole fluids on timescales long compared to the electron-electron collision time. These are given by equations (89) and (90). The two-fluid model includes Coulomb drag between the electron and hole fluids and the momentum-relaxing scattering from scattering with phonons. This then allows us to derive simple analytical expressions for the transport properties.

The plan of the paper is as follows. In section II we review the electronic structure of bilayer graphene and the low-energy effective Hamiltonian. Section III deals with the Coulomb interaction between electrons, in particular the screening thereof in the RPA approximation. We perform the calculation in flat space first and then generalize to curved space so that we can later calculate the stress tensor, which gives the response to a change in the metric. In section IV we calculate the conserved currents. We then derive the kinetic equation in section V. Section V.1 explores the effect of a constant magnetic field on the kinetic equation. In VI we write down the collision integral. Then in section VII we discuss the symmetries of the collision integral. We introduce the two-fluid model in VIII and derive analytical expressions for the transport properties in terms of the fundamental parameters of the problem. In section IX we evaluate the collision integral numerically and show the results for some transport properties of interest. We compare to the results from the two-fluid model and find good agreement. Detailed calculations are left for the Appendices.

II Review of the electronic structure of bilayer graphene

To start off, we briefly review the derivation of the band structure and the explicit expression of the wave function in BLG. This part serves to make this work self-contained and to introduce notation.

II.1 Hamiltonian

The tight-binding Hamiltonian of A-B stacked bilayer graphene with the coupling defined in Fig. 1 and external gauge field $A_{\mu}$ has the explicit form [33, 34, 19]

[TABLE]

where the velocity $v_{3}$ is given by $v_{3}=\frac{\sqrt{3}}{2}a\gamma_{3}/\hbar$, where $a$ is the lattice constant, and the Fermi-velocity is given by $v_{F}=\frac{\sqrt{3}}{2}a\gamma_{0}/\hbar$ 111We omit the spin indices here and will take the spin degeneracy into account later on, when we introduce the fermion flavor. The Zeeman splitting is negligible for the fields that we consider here ($B<100$Gauss).. Here, $\xi=1$ corresponds to $K$ valley with corresponding wave function 222Here we use the integral over $\mathbf{k}$ as the summation over $\mathbf{k}$ about each Dirac cone $K$ or $K^{\prime}$ with a cut off $\Lambda$ of order $2\pi/a$. The physical results in our paper do not depend on the cut off since we consider the physics in the long wavelength limit where $k\ll\Lambda$.

[TABLE]

and $\xi=-1$ corresponds to the $K^{\prime}$ valley with corresponding wave function $\psi_{K^{\prime}}(\mathbf{x})=(\varphi_{B_{2}}(\mathbf{x}),\varphi_{A_{1}}(\mathbf{x}),\varphi_{B_{1}}(\mathbf{x}),\varphi_{A_{2}}(\mathbf{x}))$. We have defined the momentum operator and its holomorphic and anti-holomorphic notation

[TABLE]

$e<0$ is the electron charge. In the following, we will set $v_{3}=0$, since we are only interested in the quadratic bands (see Ref. [19] for details).

II.2 Effective Hamiltonian

Since the Hamiltonian (1) provides information about both high energy and low energy states, it will be useful to create a low-energy effective Hamiltonian. To simplify our model, we consider only the low-energy bands near the $K$ valley. In the long wavelength limit $v_{F}k\ll\gamma_{1}$, one can derive for $\xi=1$

[TABLE]

where $m=\frac{\gamma_{1}}{2v_{F}^{2}}$. The wave function is given by

[TABLE]

where $\theta_{\mathbf{k}}$ is the angle between the vector $\mathbf{k}$ and the $x$-axis. $\lambda=-1$ denotes electrons in the valence band and $\lambda=1$ denotes electrons in the conduction band. In the low energy limit, we only consider the electrons appearing at sites $A_{1}$ and $B_{2}$. In the following sections, we will omit the spin indices for simplicity and consider them back in the counting factors. Similarly, we can derive the effective Hamiltonian near the $K^{\prime}$ valley, the only difference is that the wave function is now $\psi_{K^{\prime}}=(\varphi_{B_{2}}(\mathbf{x}),\varphi_{A_{1}}(\mathbf{x}))$

III Coulomb interaction and screening

Coulomb screening is the damping of the electric field due to mobile charge carriers which are quasi-particles and quasi-holes. As a result of screening, the long-range Coulomb interaction becomes short-range. In this section, we will calculate the screening effect of the Coulomb interactions in BLG, or in other words we will calculate the screening momentum $q_{TF}$.

III.1 Charge density operator

The Hamiltonian (1) shows explicitly the coupling of BLG with an external gauge field. The free Lagrangian density is given by

[TABLE]

where $\overleftrightarrow{\partial_{t}}=\frac{1}{2}(\overrightarrow{\partial_{t}}-\overleftarrow{\partial_{t}})$. The field operator in the $K$ valley in second quantization language is given by

[TABLE]

where the operator $c_{K;a}(\mathbf{k})$ is the annihilation operator of an electron on the sublattice $a$ at momentum $\mathbf{K+k}$. Similarly, the field operator in the $K^{\prime}$ valley in second quantization language can be derived. The free action in flat space is given by

[TABLE]

The total number density operator in both valleys is given by the definition

[TABLE]

From the wave function of the electron and hole bands at low energies, we can derive the transformation of the field operator

[TABLE]

and similarly for the $K^{\prime}$ bands.

Combining the spin index and the valley index to a flavor index, we obtain the effective low-energy density

[TABLE]

We can separate the effective charge density (12) into a normal part where $\lambda=\lambda^{\prime}$ and the Zitterbewegung part where $\lambda=-\lambda^{\prime}$. The homogeneous contribution to the charge density $\rho^{\text{eff}}(\mathbf{0})$ is only due to the normal part. The effective Coulomb interaction of the effective theory is given by

[TABLE]

The screened Coulomb interaction $V_{C}(q)$ will be calculated in the next subsection. The bare Coulomb interaction is

[TABLE]

and $\alpha=\frac{e^{2}}{4\pi\varepsilon_{0}}$.

III.2 Screening in flat background metric

We need to account for the screening of the long-range Coulomb interaction by the mobile charge carriers. In the random phase approximation (RPA) the dressed interaction is given by

[TABLE]

where $\varPi^{0}(q,\omega)$ is the bare susceptibility. In order to calculate $\varPi^{0}(q,\omega)$, we need to calculate the fermion loop of BLG at a finite temperature and at a given chemical potential $\mu$. This is the textbook Lindhart calculation and in the regime $\beta\mu\lesssim 1$ and $\beta q^{2}/m\lesssim 1$ we use the approximate result from Ref [37]

[TABLE]

This equation starts deviating from the full result in Ref [37] for $\beta q^{2}/m>1$, however large momentum transfer is suppressed by the Fermi occupation factor. We have that the screened potential is given according to (15) by

[TABLE]

with the screening momentum

[TABLE]

For $\beta\mu\lesssim 1$, the typical momentum is $k_{T}=\sqrt{2k_{B}Tm/\hbar^{2}}\ll q_{TF}$ for any realistic temperature, so we can safely approximate

[TABLE]

III.3 Screening in a homogeneous metric

In order to calculate the stress tensor, we need to generalize the formalism to curved space. For a homogeneous metric $g_{ij}=\delta_{ij}+\delta g_{ij}$, we can follow the steps in the last subsection and obtain the screened Coulomb interaction as

[TABLE]

where $|\mathbf{q}|=\sqrt{g^{ij}q_{i}q_{j}}$ and $g=\textrm{det}(g_{ij})$ and where $q_{TF}(|\mathbf{q}|)$ takes the form of (18). The detailed derivation of equation (20) is given in Appendix A. Equation (20) is a new result of this paper and is required in order to calculate the stress tensor operator in the next section.

IV Conserved current operators

In order to calculate the transport coefficients, we need to start by deriving the conserved current operators in the effective theory in second quantization language. The detailed derivations for the energy current and the stress tensor operators of BLG are new contributions of this paper. As has been shown in Refs [13] and [32], there are Zitterbewegung contributions to the charge current operator of graphene as well as BLG. To obtain the DC transport coefficients of BLG, one can neglect the Zitterbewegung part which is just the contribution from the off-diagonal component of the Green’s function in the generalized Kadanoff-Baym formalism 333The detailed discussions of the Kadanoff-Baym formalism are left for Appendix C.. However, this contribution will be necessary for studying the quantum transport at finite frequency and momentum using the QBE and we therefore include it for future extensions of this work. In this paper, we are only interested in spatially-independent current operators.

IV.1 Charge current operator

The current operator is by definition

[TABLE]

where $s$ stands for spin. The current density is given by

[TABLE]

where $\mathbf{J}^{\text{I}}_{\xi}(\mathbf{x})$ is the contribution of quasi-particle and quasi-hole flow, and the operator $\mathbf{J}^{\text{II}}_{\xi}(\mathbf{x})$ creates a quasi-particle-quasi-hole pair. Using the explicit wave function of low-energy modes (131), (132) and equation (21), one can derive the spatially independent current operator

[TABLE]

where we combined the spin index $s$ and valley index $\xi$ to flavor index $f$. Similarly, the Zitterbewegung current operator is given by

[TABLE]

IV.2 Energy current operator

The heat current is related to the energy current via

[TABLE]

where $\mu$ is the chemical potential and so we now calculate the energy current, which has contributions from both the kinetic and interaction energy terms in the Hamiltonian. We will follow Ref [39] in deriving the energy current operator. The conservation of energy gives us the continuity equation

[TABLE]

where $\mathbf{J}^{E}$ is the total energy current, which includes both kinetic and interaction contributions, $\mathcal{E}$ is the energy density. We will use equation (26) as the definition of the energy current.

IV.2.1 Kinetic contribution

The kinetic energy density operator is given by

[TABLE]

where $\overleftrightarrow{H}_{\xi}$ means we replace $\partial_{i}$ in $H_{\xi}$ by $\overleftrightarrow{\partial}_{i}=\frac{1}{2}(\overrightarrow{\partial}_{i}-\overleftarrow{\partial}_{i})$. We can write down the kinetic energy density in momentum space by Fourier transformation

[TABLE]

where $\hat{\mathbf{H}}_{\xi}(\mathbf{k},\mathbf{q})$ is given explicitly as follows

[TABLE]

where we defined the holomorphic and anti-holomorphic vectors $X=X_{1}+iX_{2},\qquad\bar{X}=X_{1}-iX_{2}$. Using Heisenberg’s equation $\dot{\mathcal{E}}_{kin}=-i[\dot{\mathcal{E}}_{kin},\mathcal{H}]$ and the continuity equation of energy (26) in momentum space we obtain the formula to determine the kinetic contribution to the energy current

[TABLE]

where the total Hamiltonian $\mathcal{H}$ is defined in equation (133). We leave the detailed calculation, which is quite technical, to Appendix B. We only quote here the results after taking the limit $\mathbf{q}\rightarrow 0$

[TABLE]

In comparison with the charge current operator, there is no Zitterbewegung contribution to the kinetic part of the energy current. We also see that quasi-particle and quasi-hole bands contribute to the energy current equally. At each momentum $\mathbf{k}$, quasi-particle and quasi-hole bands have the same energy and velocity and hence make the same contribution to the energy current .

IV.2.2 Interaction contribution

In the linear response calculation, the contribution of the Coulomb interaction to the energy density is given by

[TABLE]

where $N_{0}$ is the total background charge number.

The contribution to the energy current from the Coulomb interaction is then given by

[TABLE]

where $\mathbf{J}(\mathbf{q}=0)$ is nothing but the charge current. However in the kinetic formalism, we consider $\delta\mathbf{H}^{C}(\mathbf{x})$ as the shift of the chemical potential due to the background charge (the Hartree diagram).

IV.3 Stress tensor operator

The effective Lagrangian in curved space is defined as

[TABLE]

where the free Lagrangian density $\mathcal{L}_{\xi}(\mathbf{x})$ is defined in equation (7) and $\hat{V}_{C}$ is the effective Coulomb interaction. The stress tensor is defined as the response of the system with respect to a perturbation of the local metric,

[TABLE]

where $g_{ij}(\mathbf{x})=\delta_{ij}+\delta g_{ij}(\mathbf{x})$ and the rescaled field is

[TABLE]

IV.3.1 Kinetic contribution

We calculate the stress tensor operator for the kinetic Hamiltonian (1) following Ref [40]. We leave the detailed calculation to Appendix B, where we derive the results directly using definition (36) and the explicit form of the kinetic Hamiltonian in curved space. Here we quote the result of the kinetic contribution to the stress tensor

[TABLE]

where the normal contribution to the kinetic part of stress tensor is

[TABLE]

Equation (IV.3.1) has a similar form as the stress tensor operator for a quadratic semimetal [41] like HgTe.

The Zitterbewegung contribution to the kinetic part of the stress tensor is given by

[TABLE]

[TABLE]

IV.3.2 Interaction contribution

Besides the kinetic contribution, the stress tensor also has a contribution from the interactions, which we will calculate now. We turn on the homogeneous metric perturbation $\delta g_{ij}$. The Coulomb interaction in curved space near charge neutrality is given by

[TABLE]

where $\tilde{\rho}^{\text{eff}}(\mathbf{q})$ is given by substituting

[TABLE]

in (12), and $q=|\mathbf{q}|=\sqrt{g^{ij}q_{i}q_{j}}$. The factor $g^{-1/2}$ appears in equation (42) as was shown in the previous section. The transformation (43) is equivalent to the transformation (37). We now can use the definition of the stress tensor (36) to derive the contribution of the Coulomb interaction to the stress tensor in flat space time by taking the derivative of $\hat{V}_{C}$ with respect to the homogeneous metric 444There is a metric dependence in the definition of $|\mathbf{p}|$.

[TABLE]

We leave the detailed derivation of (IV.3.2) to Appendix B. The contribution to $T^{ij}_{C}(\mathbf{q}=0)$ up to linear order in the perturbation is given by555In this paper, we are only interested in the linear transport calculation.

[TABLE]

where $N_{0}$ is the background charge. We can view this contribution simply as a shift in the chemical potential. In the calculation for the shear stress tensor in section V, we will calculate $T^{12}$ under a constant shear. The contribution from the interactions $T^{ij}_{C}\sim\delta^{ij}$ will hence not enter our calculation.

V Kinetic equation and quantum transport

After having derived the conserved currents, we are now ready to begin the derivation of the actual QBE. In this vein, we will set up the semi-classical problem of electron and hole transport in bilayer graphene at a finite temperature $T$. We define the retarded Green’s function as follows 666We ignore the flavor index $f$.

[TABLE]

where $\lambda$ and $\lambda^{\prime}$ are the band indices. The expectation value $\langle\,\,\rangle$ is evaluated at finite temperature as explained in more detail in Appendix C. In order to study the DC transport, we can ignore the off-diagonal part of the retarded Green’s function since this part depends explicitly on time as explained in Appendix C. Equation (46) now takes the explicit form

[TABLE]

where $f_{+(-)}(\mathbf{p},\mathbf{x},t)$ is the distribution function of electrons in the conduction (valence) band. We can write down formally the QBE for the distribution function as

[TABLE]

where the group velocity of band $\lambda$ is defined as

[TABLE]

$\mathbf{E}(\mathbf{x},t)$ is slowly varying applied electric field. The right-hand side of the equation is the collision integral, which can be derived explicitly from first principles. In section VI, we will discuss in detail the collision integral, which takes into account the scattering of quasi-particles off each other, on impurities as well as at the boundary. The microscopic derivation of equation (48) is left for Appendix C. In the subsequent subsections, we will employ the equation (48) to set up the calculation of the transport coefficients. In DC transport, we can ignore the contribution from the Zitterbewegung contribution which comes from the off-diagonal part of the Green’s function (46). From the results (23), (25) with (31), and (IV.3.1) in the previous section, we can obtain expressions for the expectation value of the normal contribution to the conserved currents in terms of the local distribution function as follows

[TABLE]

[TABLE]

and the kinetic contribution to the stress tensor has the expectation value

[TABLE]

We note that the results in Eqns. (50)-(52) look similar to the Fermi liquid results for two types of particles, although we are in a very different regime without a well-formed Fermi surface. If we replace the distribution function $f_{\lambda}(\mathbf{k})$ in equations (50), (51) and (52) by the unperturbed Fermi distribution $f^{0}_{\lambda}(\mathbf{k})$, we get zero. In this section, we will use the above equations to obtain the expectation value of the conserved currents in terms of the distribution function perturbations:

[TABLE]

V.1 Constant applied magnetic field

So far, the experiments performed on the electrical conductivity of suspended BLG have been performed in zero magnetic field. However, we believe it is eminently possible to extend the experiments in this direction and to this end, we will set up the calculation process to obtain the transport coefficients with an applied magnetic field $\mathbf{B}=B\hat{\mathbf{z}}$. In order to use the kinetic equation with a magnetic field, we need to consider a weak magnetic field. In a Fermi liquid at zero temperature, the requirement is $k_{F}\ell_{B}\gg 1$ where the magnetic length is given by $\ell_{B}=\sqrt{\frac{\hbar c}{eB}}$. For neutral BLG, at finite temperature, one may guess that the valid limit of the kinetic equation is $k_{T}\ell_{B}\gg 1$, where the thermal momentum is defined as $k_{T}=\sqrt{2k_{B}Tm/\hbar^{2}}$. For temperature $T=10K$, the appropriate magnetic field is $B<100$ Gauss. Such a small magnetic field also guarantees that the Zeeman energy term is small enough, that we can neglect the energy different between the two spin species.

With the appearance of a magnetic field, we need to add one more term in the left-hand side of the kinetic equation to take into account the Lorentz force [45]777We ignore the Zeeman effect, since it is small in comparison with the experimental temperature [18].

[TABLE]

where the group velocity is given by (49). In this section, we only consider the charge conductivity and thermal conductivity in the appearance of a magnetic field. We can rewrite (54) as

[TABLE]

V.2 Thermoelectric coefficients

We define the electrical conductivity $\sigma$, the thermal conductivity $K$ and the thermopower $\Theta$ by

[TABLE]

where each of the thermoelectric coefficients is a $2\times 2$ matrix. The fact that $\Theta$ appears twice in (56) is due to the Onsager reciprocity relation [47]. Using these definitions, the Seebeck coefficient is $S=\sigma^{-1}\Theta$ and the Peltier coefficient is $\Pi=TS$. In experiments, the heat current is often measured such that $\mathbf{J}=0$, in which case the proportionality constant between $\mathbf{J}^{Q}$ and $-\mathbf{\nabla}T$ is $\kappa=K-T\Theta\sigma^{-1}\Theta$ [14, 2].

V.3 Charge conductivity

In order to derive the coefficients of DC conductivity, we apply a constant electric field $\mathbf{E}$. The unperturbed distribution function is given by

[TABLE]

We need to solve the equation (48) in the following simplified form

[TABLE]

for $\lambda=+$ and $\lambda=-$ to obtain $h_{\lambda}(\mathbf{p})$. In equation (V.3), the right-hand side denotes the linear order in the perturbation of the collision integral derived in section VI. The left-hand side is derived in the Green’s function formalism as (180). The suggested ansatz in this calculation is

[TABLE]

and we solve for $\chi_{\lambda}(p)$ numerically. The second term in the ansatz becomes relevant when we have a magnetic field. The charge current is given by (50) and the DC conductivity can be directly read off. Due to the symmetry of the collision integral that will be discussed in section VII, we can show that $\sigma_{xx}=\sigma_{yy}$ because of rotational invariance and $\sigma_{xy}=\sigma_{yx}=0$ in the absence of magnetic field due to parity. The external magnetic field $\mathbf{B}$ breaks parity which gives us $\sigma_{xy}=-\sigma_{yx}\neq 0$.

V.4 Thermal conductivity

We consider a spatially dependent background temperature $T(\mathbf{x})=T+\delta T(\mathbf{x})$. The local equilibrium distribution function takes the form

[TABLE]

We now consider a constant gradient in temperature by introducing the space-time independent driving force $\mathbf{F}^{T}=-\frac{\nabla_{\mathbf{x}}\delta T}{T}$. We then need to solve equation (48) in the following simplified form

[TABLE]

for $\lambda=+$ and $\lambda=-$ to obtain $h_{\lambda}(\mathbf{p})$. The left-hand side is obtained from (185). The suggested ansatz is

[TABLE]

From the heat current along with equations (51) we can read off the thermal conductivity. For the thermopower, we consider the same ansatz as for the thermal conductivity, but calculate the charge current which is given by (50).

V.5 Viscosity

To calculate the DC shear viscosity, we consider a background velocity for the particles and holes. Therefore, the local equilibrium distribution function takes the form

[TABLE]

where $\mathbf{u}_{+(-)}(\mathbf{x})$ is the perturbed background velocity of electrons (holes). We apply a constant shear with the explicit form

[TABLE]

where $F_{\lambda}$ is a space-time independent perturbation and the definition of strain is

[TABLE]

We need to solve equation (48) in the following simplified form

[TABLE]

for $\lambda=+$ and $\lambda=-$ to obtain $h_{\lambda}(\mathbf{p})$. The left-hand side comes from (193). The suggested ansatz is

[TABLE]

The stress tensor $T^{12}_{\lambda}$ is given by (52) and the shear viscosity coefficients are given by the definition

[TABLE]

In the experiment [48], the authors found that quasi-particle collisions can importantly impact the transport in monolayer graphene. The results showed that the electrons behave as a highly viscous fluid due to the electron-electron interactions in the clean limit. Even though there has not been an analogous experiment for BLG yet, we expect that highly viscous behaviour of BLG will be found in the near future. The viscosity coefficients will play an important role for simulation of electronic transport in BLG and for comparison against experimental results.

VI Collision integral

We now focus on the right-hand side of the QBE—the collision integral. We discuss the contribution from quasi-particle interactions, scattering on disorder and scattering off the boundary separately in subsections VI.1, VI.2 and VI.3 respectively. In the quasi-particle scattering channel, we ignore Umklapp processes at low energies near charge neutrality. Since in our regime $k_{F}a\ll 1$, Umklapp scattering is negligible due to the lack of available phase space. Inter-valley scattering is also ignored due to the long range nature of the Coulomb interaction. Up to linear order in the perturbation, the generalized collision integral on the right-hand side of the kinetic equation (48) includes contributions from quasi-particle interactions, scattering on disorders and finite size effect and scattering on phonons respectively

[TABLE]

In the following subsection, we will discuss in detail each contribution of (69).

VI.1 Quasi-particles’ Coulomb interaction

The first contribution to the collision integral that we consider is that coming from the Coulomb interaction of the quasi-particles. We are interested in the experimental regime of sufficiently clean BLG [18] in which the transport properties are dominated by quasi-particle interactions. In this section, we formulate the quasi-particle interactions of the form (13) via the screened Coulomb potential $V_{C}(\mathbf{q})$ that was derived previously in section III. To derive the contribution $I^{(1)}_{\lambda,\textrm{int}}$, we generalize the Kadanoff-Baym equations [9] to BLG. We again only consider the diagonal component of the Green’s functions (47) and calculate the collision integral contribution due to the interaction (13). The technical details of the derivation will be left for Appendix C, in this subsection we only quote the main result. The collision integral due to interactions for each band index $\lambda$ is then given by

[TABLE]

where we follow [13] and define the form factor

[TABLE]

as well as the channel dependent scattering matrix

[TABLE]

The collision integral vanishes when we substitute the Fermi distribution (57). Linear order in perturbation of the collision integral is given by the perturbation (53). The linearized collision integral is then given by

[TABLE]

where we defined $h_{\lambda}(\mathbf{k})$ in equation (53). The collision integral (73) shares similarities with one of monolayer graphene in Ref [13]. However, due to the difference in the quasi-particle dispersion relation, which is quadratic for BLG and linear for monolayer graphene, their allowed scattering channels differ qualitatively. In the case of BLG, we have to consider the scattering channel where one quasi-particle decays to two quasi-particles and one hole. On the other hand, this scattering channel is kinematically forbidden in monolayer graphene because of momentum and energy conservation. This contribution was missed in a previous publication on the kinetic theory of BLG [32]. However, due to kinematic restrictions, the phase space for this scattering process is small and therefore this channel does not contribute significantly to the collision integral. We have checked this statement numerically.

VI.2 Contribution from disorder

Due to the Galilean invariance of our system in the absence of disorder, the collision integral is unchanged under a Galilean boost. However, under a uniform boost of all particles by $\mathbf{u}$, the current density transforms as $\mathbf{J}\to\mathbf{J}+en\mathbf{u}$. So as long as the charge density $n\neq 0$ (ie $\mu\neq 0$) we change the current density by boosting frames and therefore the conductivity is ill-defined in the absence of a momentum-relaxing mechanism. Including one or several momentum-relaxing scattering channels is therefore crucial for calculating the transport coefficients away from $\beta\mu=0$.

One such momentum relaxing process is the scattering of electrons off impurities in the sample. For this calculation, we put our system in a box of side length $L$ with periodic boundary conditions. We follow [13] and consider a disorder Hamiltonian

[TABLE]

where $V_{\textrm{dis}}$ is the interaction potential between an electron and the impurities, which we take to be charges $Ze$ located at random positions $\mathbf{x}_{i}$ and having number density $n_{\textrm{imp}}=N_{\textrm{imp}}/L^{2}$. We use the screened Coulomb interactions to obtain

[TABLE]

From the interaction (74), we can calculate the scattering rate of quasi-particles off the disorder. We then obtain the contribution to the collision integral from disorder up to linear order in the perturbation

[TABLE]

where we define a short hand notation for the impurity scattering rate

[TABLE]

The corresponding dimensionless parameter is $\alpha_{\textrm{imp}}\equiv\beta\tau_{\textrm{imp}}^{-1}=1/2(8\pi^{2}Z/N_{f}\epsilon_{r})^{2}\beta n_{\textrm{imp}}/m$. The detailed derivation of (76) is left for Appendix D.

VI.3 Effect of finite system size

In very clean samples of bilayer graphene, it is expected that the scattering length due to impurity scattering is longer than the system size $L$, which is currently limited in suspended graphene samples [28]. In this case, in order to have a well-defined conductivity, we need to include the effect of the finite size of the system. There will be scattering of the electrons off the boundary, which effectively acts as an additional scattering time. Assume the scattering time due to collisions with the boundary is $\tau(p)=\frac{L}{v}=\frac{Lm}{p}$ where $L$ is the size of the sample up to a factor depending on the geometry of the BLG sample. Here, we are making the simplifying assumption, that that the scattering does not depend on the direction of the momentum. We neglect the angular dependence of the boundary scattering which is likely to contribute a geometric factor to the scattering time. The collision integral is then

[TABLE]

which is just the Bhatnagar-Gross-Krook (BGK) collision operator [49] with $\tau$ given by $\tau(p)$. The corresponding dimensionless parameter is $\alpha_{L}=\frac{\sqrt[]{\beta}}{\sqrt[]{m}L}$.

VI.4 Phonon scattering

We should also consider the effect of the electrons scattering off phonons. The maximum energy of an acoustic phonon is $\varepsilon_{\textrm{max}}=2c\sqrt{2m\textrm{max}(k_{B}T,\mu)}$, where $c$ is the speed of sound in graphene. In the experimental setting, we are at high temperatures compared to the Bloch-Grüneisen temperature

[TABLE]

and additionally we have $T\gg 2\gamma_{1}/k_{B}(c/v)^{2}$. Thus we are in the high temperature regime $k_{B}T\gg\varepsilon_{\textrm{max}}$, where can treat the phonons as introducing another scattering time [50]888We ignore the scattering channel from conductance electrons to valence electrons due to the emission and absorption of phonons because of the suppression in the scattering matrix.[50].

[TABLE]

where $D$ is the deformation potential and $\rho$ is the mass density. Then the collision integral is

[TABLE]

The corresponding dimensionless parameter is $\alpha_{\textrm{ph}}=\beta\tau_{\textrm{phonon}}^{-1}$. It is crucial to note that whereas $\alpha_{\textrm{imp}}=\alpha_{\textrm{imp}}(T)$ and $\alpha_{\textrm{L}}=\alpha_{\textrm{L}}(T)$, $\alpha_{\textrm{ph}}$ does not depend on temperature.

VII Symmetries

VII.1 Spatial symmetries

The electrical conductivity is rotationally symmetric. The only rotationally symmetric tensors in 2d are $\delta_{ij}$ and $\epsilon_{ij}$ so any rotationally invariant tensor $\sigma_{ij}$ can be written as

[TABLE]

In the absence of a magnetic field, we have an additional symmetry, namely 2D parity $y\to-y$, which implies

[TABLE]

With a magnetic field, 2D parity implies

[TABLE]

The thermal conductivity and thermopower satisfy the same relations.

VII.2 Particle-hole symmetry

Under the particle-hole transformation, we have $\lambda\to-\lambda$, $\mu\to-\mu$ and $B\to-B$. First consider the electrical conductivity. Due to the particle-hole symmetry, we have

[TABLE]

Now consider the viscosity which we only calculate for $B=0$. Particle-hole symmetry implies that we have

[TABLE]

These symmetries follow directly from the form of the collision integral (73).

VIII Two-fluid model

We introduce the two-fluid model, which reproduces the salient features of our numerical results. Motivated by comparison with experiment [52] we choose to only include the phonon scattering as a momentum relaxing mechanism. We multiply the kinetic equation by $\lambda\mathbf{p}/m$ and integrate over momentum space in order to derive the evolution of the mean fluid velocities as

[TABLE]

[TABLE]

where we defined the electron and hole velocities as

[TABLE]

The coefficients $\Lambda^{e,h}$ account for the fact that the average entropy per particle is $\Lambda k_{B}$

[TABLE]

[TABLE]

The definitions (92) and (93) follow from the $\nabla T$ term in the QBE when (89) and (90) are derived. The Coulomb drag term can be derived explicitly from the collision integral

[TABLE]

This allows us to calculate $\tau_{eh}$ and $\tau_{he}$. We perform the calculation at charge neutrality and then use (108) to extrapolate. $\tau_{se}$ is the momentum-relaxing scattering time for electrons and $\tau_{sh}$ is the corresponding time for holes ($s$ stands for ”scattering”). They are given by

[TABLE]

[TABLE]

where where $\tau_{\textrm{phonon}}$ and $\tau_{\textrm{imp}}$ are given by (80) and (77) respectively and the scattering times off the boundary are

[TABLE]

[TABLE]

We consider the steady state $\partial_{t}\mathbf{u}^{e}=\partial_{t}\mathbf{u}^{h}=0$ and calculate the electric current and energy current

[TABLE]

[TABLE]

where the number densities calculated from the Fermi distribution are

[TABLE]

From this, we can derive the thermoelectric coefficients in the absence of a magnetic field

[TABLE]

[TABLE]

[TABLE]

where

[TABLE]

[TABLE]

For momentum conservation we require

[TABLE]

We verify explicitly that the Onsager relations for the thermoelectric coefficients are satisfied if equation (107) is satisfied. Thus we can choose

[TABLE]

This ansatz agrees with the full numerical result obtained from (94) to within 10% in the entire range of $\beta\mu$ and is therefore a satisfactory approximation. By evaluating the collision integral in (94) numerically, we find

[TABLE]

We define the dimensionless electrical conductivity as $\sigma_{ij}=\frac{N_{f}e^{2}}{2\hbar}\tilde{\sigma}_{ij}$. With a magnetic field and at CN, we calculate from the two fluid model that the Hall conductivity at small fields behaves like

[TABLE]

where $\alpha_{s}=\alpha_{\textrm{imp}}+\alpha_{\textrm{ph}}$ and we have neglected the boundary scattering. We plot this quantity in Fig. 3 and show that the result from the two-fluid model agrees perfectly with the numerical result.

IX Numerical results

Armed with the full formalism for the QBE, we are in a position to numerically calculate the transport properties. In our companion paper [52] we plot the three thermoelectric coefficients as a function of $\beta\mu$. In this section we therefore focus on the behaviour of the transport coefficients at CN and in a magnetic field and we discuss the viscosity.

IX.1 Transport at CN

In Fig. 2 we show how the electrical conductivity at charge neutrality depends on temperature. In order to obtain a non-trivial temperature dependence we need to go beyond the Coulomb and phonon scattering. We assume that collisions off impurities can be neglected, as claimed in the experimental work [18]. With Coulomb interactions and phonons alone, the conductivity at CN would be independent of temperature, since at CN, the conductivity would only depend on the dimensionless parameter $\alpha_{\textrm{ph}}$, which is temperature-independent. So the temperature-dependence is entirely due to the finite size scattering, which comes with the dimensionless number $\alpha_{L}$, which does depend on temperature. This figure shows qualitative agreement with the behaviour seen in Fig. 4 of [18]. The thermal conductivity at charge neutrality shows the same type of behaviour as the electrical conductivity and for the same reasons.

IX.2 Lorenz number

From the Lorenz number $\mathcal{L}=\kappa_{xx}/\sigma_{xx}T$ and the Hall Lorenz number $\mathcal{L}_{H}=\kappa_{xy}/\sigma_{xy}T$ we deduce a further signature of the two-fluid model. The Lorenz number is enhanced relative to the Wiedemann-Franz (WF) law which predicts $\mathcal{L}=\pi^{2}/3(k_{B}/e)^{2}$. The violation of the WF law has been reported in a recent theoretical work [53]. On the other hand, the violation of the WF law is much less severe for the Hall Lorenz number $\mathcal{L}_{H}$. Both these observations can be explained in the following simple picture. We find the Lorenz number at charge neutrality

[TABLE]

where $\Lambda^{e}=\Lambda^{h}\equiv\Lambda$ at CN and $\tau_{s}^{-1}=\tau_{\textrm{imp}}^{-1}+\tau_{\textrm{phonon}}^{-1}$. We neglect scattering off the boundary in this section. From Drude theory, $\kappa_{xx}\propto\tau_{\kappa}$, where $\tau_{\kappa}^{-1}$ is the scattering rate due to all collisions that relax the energy current. At CN for an applied thermal gradient, $\mathbf{u}^{e}=\mathbf{u}^{h}$ so there is no Coulomb drag between the electrons and holes and hence only the momentum relaxing scattering limits the thermal conductivity $\tau_{\kappa}^{-1}=\tau_{s}^{-1}.$ On the other hand the Coulomb drag is important for the electrical conductivity, hence $\sigma_{xx}\propto\tau_{\sigma}$ where $\tau_{\sigma}^{-1}=\tau_{eh}^{-1}+\tau_{he}^{-1}+\tau_{s}^{-1}$, where we have added the scattering rates according to Matthiessen’s rule 999Matthiessen’s rule is only valid at CN, where the electrons and holes have equal densities and scattering times.. This immediately yields equation (111).

The Hall Lorenz number on the other hand is

[TABLE]

We can again derive this result using simple arguments. For an applied thermal gradient, electrons and holes move in the same direction, thus the only scattering along the direction of the gradient, x, is $\tau_{s}^{-1}$. Electrons and holes will be deflected in opposite directions by the applied magnetic field, so the friction in the perpendicular y direction is $\tau_{s}^{-1}+\tau_{eh}^{-1}$. This will increase the friction between the two fluids and limit the value of $\kappa_{xy}$. For an applied electric field, the electrons and holes move in opposite direction, so the friction in the x direction is $\tau_{s}^{-1}+\tau_{eh}^{-1}$. The magnetic field will deflect them in the same direction and the two fluids will feel the reduced friction $\tau_{s}^{-1}$ between them in the y direction. From these considerations, (112) can be shown.

Thus,

[TABLE]

From the numerics we indeed find $\mathcal{L}\approx 25(k_{B}/e)^{2}$ and $\mathcal{L}_{H}\approx 6(k_{B}/e)^{2}$.

IX.3 Viscosity

We calculate the shear viscosity as the response of the stress tensor when a shear flow is applied to either of the particle species. The viscosity tensor $\eta_{\lambda\lambda^{\prime}}$ is then defined via $T_{\lambda}^{12}=-\eta_{\lambda\lambda^{\prime}}F_{\lambda^{\prime}}$ and $\eta_{\lambda,\lambda^{\prime}}=\frac{N_{f}m}{\beta}\tilde{\eta}_{\lambda,\lambda^{\prime}}$. $\eta_{\lambda,\lambda^{\prime}}$ is a measure of the friction between particle species $\lambda$ and $\lambda^{\prime}$. In the numerical data FIG. 4 we see that at large $\beta\mu$, $\eta_{++}$ dominates, since electron-electron collisions are the dominating ones. Conversely $\eta_{+-}$ decreases at large $\beta\mu$ since there are less holes present to exert a friction on the electrons. We also note that $\eta_{+-}\ll\eta_{++}$. This can be understood from the kinematics of collisions. Energy and momentum conservation constrain the available phase space more for electron-hole collisions than for electron-electron collisions. In addition, the matrix elements for electron-hole collisions favour large momentum exchange, which is suppressed by the potential. This justifies the two-fluid model since this shows that the intra-fluid collisions of the electron and hole fluids dominate over the inter-fluid collisions and therefore we can treat the two fluids as weakly interacting.

A possible probe of the viscosity is via the negative non-local resistance [4, 5]. In order to make this measurement quantitative, the relation between viscosity and negative nonlocal resistance must be determined which requires a full solution of the fluid equations in the relevant geometry. Another method for measuring the viscosity of graphene has been proposed using a Corbino-disc device [55].

We note that the famous KSS result [56] provides a lower bound for the ratio of the shear viscosity $\eta$ to the entropy density $s$ in a strongly interacting quantum fluid, $\eta/s\geq 1/(4\pi k_{B})$. In our case the entropy density is $s=n^{e}k_{B}\Lambda^{e}+n^{h}k_{B}\Lambda^{h}=(2\pi/3)mk_{B}/\beta$ and so

[TABLE]

since $\tilde{\eta}\gtrsim 0.5$ from Fig. 4. Since the bound is saturated for an infinitely strongly coupled conformal field theory, the fact that we are away from the bound is consistent with the previous arguments that we are in a weakly coupled regime and the semiclassical method is valid.

IX.4 Detailed benchmarking

In order to assess the usefulness of the two-fluid model, we now perform a detailed analysis of the agreement between the QBE and the two-fluid model for a large range of the parameter-space of the problem. In the top rows of Figs. 6 and 7 we check the agreement for phonon or impurity scattering, which is described by dimensionless strength $\alpha_{\textrm{ph}}$ or $\alpha_{\textrm{imp}}$. These two cases are identical in both the QBE and two-fluid model at fixed $\alpha$, they only differ in the temperature dependence of the dimensionless number $\alpha$. The bottom row shows the corresponding results for finite-size scattering, which is described by the dimensionless strength $\alpha_{\textrm{L}}$.

We see from Figs. 6 and 7 that in general the agreement is very good for weak momentum-relaxing scattering, ie. small values of $\alpha_{\textrm{ph}}$ or $\alpha_{L}$. In this limit, the Coulomb-mediated electron-electron collisions are dominant and the hydrodynamic description works well. For larger values the agreement gets worse, especially in the case of the thermal conductivity. This is due to the fact that our two-fluid model only includes the equation for the first moment of the QBE. To get the thermal conductivity accurately, one would have to include the second moment as well, however this renders the two-fluid model too complicated to solve analytically, defeating the purpose of introducing it in the first place. We note however, that there is no reason to trust either the QBE or the two-fluid model in the strongly coupled regime $\alpha\gtrsim 1$. We also note the at least for the conductivity, the agreement is significantly better for the phonon and impurity scattering compared to the boundary scattering. The reason is that the boundary scattering has a scattering time that depends on momentum and in the two-fluid model we parametrize this by an average scattering time.

We note that even when the agreement with the two-fluid model fails, our QBE solution still satisfies the symmetries listed in section VII, since these are exact symmetries of the QBE. We have checked our numerical solution to find that the symmetries are indeed obeyed to an accuracy of $10^{-7}$.

X Conclusion

This paper sets up the Quantum Boltzmann formalism for bilayer graphene. It will serve as a reference work for numerical studies of the QBE that can be compared to experimental results. The experimentally-relevant transport quantities that we focus on are the thermo-electric coefficients (electrical conductivity, thermopower, thermal conductivity) and the shear viscosity. So far, only the electrical conductivity has been measured in experiment.

The calculation of the transport coefficients requires two ingredients: Firstly, we need to calculate the conserved currents associated with the coefficients in terms of the distribution function. In the case of the viscosity for instance, one has to calculate the stress tensor. This requires working out the coupling of BLG to a curved background metric, a calculation that is performed in this paper for the first time. Secondly, we need to work out the change in the distribution function due to the applied external fields. We use the Kadanoff-Baym equations as a starting point. The most technical part of this derivation is the calculation of the collision integral, which is performed in detail in the appendices. Once we have these ingredients, we can plug the change of the distribution function into the expressions for the conserved currents to find the linear response to the applied external fields. This allows us to read off the transport coefficients.

The dominant term in the collision integral in the hydrodynamic regime of BLG will be the Coulomb interactions. However, in order to obtain a finite conductivity, we need to break the Galilean invariance of the system. There are three possible terms that can be added to the collision integral: the scattering of the electrons off (1) phonons, (2) impurities, (3) the boundary of the sample. Depending on the experimental parameters, one or the other may dominate and in this work we have calculated all three contributions.

In the case of monolayer graphene, the electrons obey a linear dispersion relation. Energy conservation together with momentum conservation then places tight constraints on the phase space of collisions and this allows analytical results for the collision integral to be obtained [13]. A similar simplification in the case of BLG is not possible due to the quadratic energy dispersion. Due to the analogy with monolayer graphene, some previous authors have neglected scattering terms that are forbidden for monolayer graphene but allowed for BLG. We explicitly included these terms in our work. The collision integral must be evaluated numerically.

We derived from the QBE the two-fluid model — a simple hydrodynamic model for the evolution of the mean fluid velocity of the electron and hole fluids. There is Coulomb drag between the two fluids and they are both subject to scattering off phonons. This model is simple enough to be able to obtain analytical formulae for the transport coefficients. We show that the two-fluid model provides quantitatively accurate results in the hydrodynamic regime where the electron-electron collisions are dominant and momentum-relaxing collisions are subdominant.

Our predictions for the temperature and magnetic field dependence of the conductivity can be verified experimentally. It should be possible to add a magnetic field to the experiment and perform the measurement of the electrical and thermal conductivities. This is another interesting probe of the hydrodynamic regime in BLG and can be used to check the agreement between the experimental behaviour and the theoretical predictions.

Our formalism can be adapted to consider BLG far from charge neutrality by modifying the screening calculation. Due to the quadratic dispersion, we expect that in the $\beta\mu\gg 1$ regime, one recovers the standard Fermi liquid results.

It is also possible to generalize the formalism to treat multilayer graphene. A further possible avenue of research is to extend the present formalism to finite frequencies. Besides adding an extra term to the collision integral corresponding to the time-derivative in the Boltzmann equation, this may require taking into account the off-diagonal components of the probability distribution matrix as well as considering the Zitterbewegung contributions.

Finally, another direction of research is the calculation of the Hall viscosity, which has been measured experimentally in monolayer graphene and BLG experiments [57].

XI Acknowledgements

We would like to acknowledge helpful discussions with Philipp Dumitrescu. We thank Christopher Herzog and Dam Thanh Son for their comments on the KSS bound. This work was supported by EP/N01930X/1 and EP/S020527/1. Statement of compliance with EPSRC policy framework on research data: This publication is theoretical work that does not require supporting research data.

Note added: During the completion of this work we became aware of related work [53] which looks at the thermoelectric properties for BLG. However, in our work, we present a detailed derivation of the quantum Boltzmann formalism and we have a different form of the collision integral. In addition, we deduce the effect of the experimentally relevant finite-size effect as well as phonon scattering. We aim to provide quantitative results that can be compared to current and future experimental data.

Appendix A Detailed calculation for Coulomb screening in a homogeneous metric

In this subsection, we will show how to modify the calculation for screening momentum in the homogeneous metric

[TABLE]

where $\delta g_{ij}$ is space-time independent. The effective Hamiltonian is given by

[TABLE]

where here we define

[TABLE]

Note that this definition only works for a homogeneous metric, for a general metric, we need to replace $|\mathbf{x}-\mathbf{y}|$ by the geodesic distance. We follow [40] and define the rescaled field as

[TABLE]

We can rewrite the Hamiltonian in the momentum space as

[TABLE]

where the Fourier transformation is given as follows

[TABLE]

and we define

[TABLE]

The screened Coulomb interaction is given by

[TABLE]

where $V(|\mathbf{q}|)$ can be read off from (119)

[TABLE]

From the Hamiltonian (119), we see that the propagator of $\tilde{\hat{\Psi}}$ is given by

[TABLE]

We can repeat the calculation in section III to obtain the susceptibility

[TABLE]

From the definition (121), we introduce the new variable

[TABLE]

and obtain

[TABLE]

By the means of transformation (126), we have

[TABLE]

Combining equations (127), (128), and (125) we arrive at

[TABLE]

Thus we have

[TABLE]

Note that the calculation in this subsection only valid for a homogeneous metric. The screening potential for general metric should be calculated in a completely different manner.

Appendix B Detailed derivation of energy current and stress tensor operators

In this Appendix, we will present the detailed derivation of equation (31) and (38). he exact low energy wave function can be calculated directly by diagonalizing the Hamiltonian (1), we obtain

[TABLE]

[TABLE]

where $k=k_{x}+ik_{y}$ and $\bar{k}=k_{x}-ik_{y}$. We see that when we take the approximation $v_{F}|\mathbf{k}|\ll\gamma_{1}$, we recover the low energy results up to a gauge.

B.1 Derivation of (31)

In order to calculate the commutator in (30), we need to write down the explicit form of the Hamiltonian in second quantization language

[TABLE]

with $\hat{H}_{\xi}(\mathbf{k})=\hat{\mathbf{H}}_{\xi}(\mathbf{k},0)$ and

[TABLE]

Combining equations (29) and (133) and using the anti-commutation relation

[TABLE]

we have

[TABLE]

We now will calculate the last contribution. We can write down explicitly the commutator as

[TABLE]

Up to linear order in perturbation, and due to the fact that $\hat{h}_{\xi aa}=0$ , the only nonzero contribution comes from $\mathbf{q}^{\prime}=0$. We can rewrite the above equation as

[TABLE]

where $N_{0}$ is the background charge. The contribution of the kinetic part to the energy current comes from the first term of equation (136). We can read off the spatially independent current density by taking the $\mathbf{q}\rightarrow 0$ limit of equation (30)

[TABLE]

where

[TABLE]

[TABLE]

The contribution from low energy bands can be calculated explicitly by substituting wave functions (131), (132) into the above equation and then making the approximation $v_{F}|\mathbf{k}|\ll\gamma_{1}$. We can obtain the kinetic part of the energy current by adding the contribution from $K$ and $K^{\prime}$ valleys to derive (31).

B.2 Derivation of (38)

In this subsection, we will present the detailed derivation of (38). We rewrite the definition (36) in terms of the vierbein instead of the metric as follows. First, we introduce the vierbein $e^{i}_{a}$ with following definition

[TABLE]

Using the definition (36) and the symmetry of the stress tensor operator $T^{ij}$, we have the new definition of $T^{ij}$

[TABLE]

The zero momentum component of the stress tensor $T^{ij}(\mathbf{q}=0)$ is given by the response of the system to a homogeneous perturbation of the local metric $\delta g_{ij}(\mathbf{x})=\delta g_{ij}$. From the definition (36), the only contribution to the stress tensor $T^{ij}(\mathbf{q}=0)$ comes from the Hamiltonian 101010We see that the first term of the effective action $i\tilde{\hat{\Psi}}_{\xi}^{\dagger}\partial_{t}\tilde{\hat{\Psi}}_{\xi}$ has no variation with respect to the homogeneous perturbation of the metric.. The coupling of the kinetic part of the Hamiltonian with the spatially uniform vierbein is given by 111111Since we only consider a spatially independent perturbation of the metric, the spin connection vanishes.

[TABLE]

with

[TABLE]

Combining (143) and (144), one can derive the kinetic contribution to the spatially independent part of the stress tensor operator from each valley

[TABLE]

We plug in the explicit form of wave functions (131) and (132) of the low energy band and use the approximation $v_{F}|\mathbf{k}|\ll\gamma_{1}$ in the above equation and obtain (38).

B.3 Derivation of (IV.3.2)

In this subsection, we will present the detailed derivation for the interaction contribution to the stress tensor. We rewrite the Coulomb interaction (42)

[TABLE]

where in $\tilde{\rho}^{\text{eff}}(\mathbf{p})$, we replace the field operator $c_{\lambda f}(\mathbf{k})$ by the rescaled one with homogeneous metric perturbation

[TABLE]

The metric dependence of $\hat{V}_{C}$ in the linear transport formalism is in the screened Coulomb potential

[TABLE]

We take the derivative of $\mathcal{V}_{C}(\mathbf{p})$ with respect to $g_{ij}$ and obtain

[TABLE]

We then plug (150) and (147) into the definition (36) with $\delta g_{ij}(\mathbf{x})=\delta g_{ij}$ to obtain (IV.3.2).

Appendix C Density matrix formalism for BLG

C.1 Generalizing the quantum Boltzmann equation to charge conductivity

In this section, we will derive the density matrix formalism for BLG. The derivation follows closely Ref [13] and the classic book [9]. The effective Hamiltonian for each flavor is given by

[TABLE]

We begin with the modification of equations (8.27) and (8.28) of Ref. [9]. We define the density matrix for each flavor as

[TABLE]

[TABLE]

where $a,b=1,2$ are the sub-lattice indices, $f=1,..,4$ is the flavor index. We define the perturbed expectation value as

[TABLE]

and we define the unitary transformation as

[TABLE]

where $U(\mathbf{x},t)$ is the applied scalar potential and $\rho(\mathbf{x},t)$ is the density operator. By construction, we notice that

[TABLE]

since the perturbation and interactions conserve the flavor.

The Green’s functions $G^{<,>}_{ab}(1,1^{\prime},U)$ satisfies the following, so-called Kadanoff-Baym equations of motion

[TABLE]

[TABLE]

in which we omit the indices, the equations (157) and (158) need to be considered as matrix equations. The self-energy matrix in Born collision approximation is given by

[TABLE]

Equation (159) is a matrix equation, both left and right sides are $8\times 8$ matrix.

At this point we are ready to derive the equation of motion for the density matrix. We subtract (158) from (157) to obtain

[TABLE]

We want to use the approximation of the Green’s function $G^{\lessgtr}(\mathbf{x}_{1},t_{1},\mathbf{x}_{1^{\prime}},t_{1^{\prime}};U)$ for slowly varying applied potential $U(\mathbf{R},T)$ as a function of

[TABLE]

We also want to consider $G^{\lessgtr}(\mathbf{x}_{1},t_{1},\mathbf{x}_{1^{\prime}},t_{1^{\prime}})$ to be sharply peaked about $\mathbf{r}=0$ and $t=0$, where

[TABLE]

We can rewrite the Green’s function as

[TABLE]

In the DC case, we only consider the static component of Green’s function in which satisfies

[TABLE]

We consider the left hand side of of (160), which we can rewrite as

[TABLE]

We Fourier transform the relative coordinate $\mathbf{r}$ and $t$ by multiplying by $e^{-i\mathbf{p}\cdot\mathbf{r}+i\omega t}$ and integrating over $\mathbf{r}$ and $t$ to obtain

[TABLE]

For an applied static electric field, we have

[TABLE]

Let’s look at the approximation that we applied to the Green’s function more closely. In the Weyl-Wigner formulation, we only consider the Green’s $G^{<}(\mathbf{r},t,\mathbf{R},T;U)$ that is slowly varying in $\mathbf{R}$. Physically, it means that we only consider the perturbation such that

[TABLE]

is slowly varying in $\mathbf{R}$. Fourier transforming (168), we obtain

[TABLE]

where $\mathbf{k}$ ($\mathbf{K}$) is the momentum conjugate to $\mathbf{r}$ ($\mathbf{R}$). Since we are interested in spatially homogeneous distributions, that means we set $\mathbf{K}=0$ and only consider the Green’s function of the form

[TABLE]

We now can convert the Green’s function with sub-lattice indices to the Green’s function with band indices using the relation

[TABLE]

where $\mathcal{U}^{f}(\mathbf{k})$ and $\mathcal{U}^{f\dagger}(\mathbf{k})$ can be read off from the explicit form of the band wave function in Section II

[TABLE]

We can transform the equation (166) for each flavor

[TABLE]

where we omit $\mathbf{R}$ and $T$ and replace $\partial_{\mathbf{R}}=0$ and $\partial_{T}=0$, we also denote the Green’s function in band indices as

[TABLE]

In the Weyl-Wigner formulation, the LHS of the equation for $G^{>}(\mathbf{p},\omega;U)$ can be written similarly as

[TABLE]

In our calculation, we are interested in the DC transport, in which we omit the off-diagonal part of the density matrix due to the condition (164) 121212Because of the time-dependence of $c^{(\dagger)}_{\lambda}(\mathbf{k})$, the off-diagonal component of the Green’s function $g^{\lessgtr}$ depends on the central time $T$.. We then linearize the kinetic equation up to linear order in perturbation. The equation (173) is rewritten as

[TABLE]

where $g_{0}^{<}(\mathbf{p},\omega)$ is given by

[TABLE]

where fermionic distribution functions is given by

[TABLE]

We integrate over $\omega$ to obtain the equal time Green funtion, the equation (177) becomes

[TABLE]

C.2 Generalizing the quantum Boltzmann equation to thermal conductivity

In order to derive the kinetic equation for thermal conductivity, we turn on a gradient of the temperature $T=T(\mathbf{R})$. The local equilibrium distribution function is given by

[TABLE]

The equation of motion for $G^{<}(1,1^{\prime};T(\mathbf{R}))$ is given by (160). We again use Weyl-Wigner formulation and rewrite the left hand side of equation (160) as

[TABLE]

If we ignore the off diagonal part of density matrix, equation (182) can be rewritten as

[TABLE]

Up to linear order in perturbation, we replace $g^{<}(\mathbf{p},\omega,\mathbf{R};T(\mathbf{R}))$ by the equilibrium one

[TABLE]

Again, we integrate over $\omega$ to obtain the equal time Green’s function. We consider a constant gradient in temperature by introducing the space-time independent driving force $\mathbf{F}^{T}=-\frac{\nabla_{\mathbf{R}}T}{T}$, the equation (183) becomes

[TABLE]

C.3 Generalizing the quantum Boltzmann equation to shear viscosity

In order to derive the kinetic equation for shear viscosity, we assume that the particles and holes have a spatially dependent local velocity. We consider the local equilibirium distribution function

[TABLE]

We can follow the last subsection to obtain the kinetic equation for $G^{<}(1,1^{\prime};\{\mathbf{u}_{\lambda}\})$ (160). We can use Weyl-Wigner coordinate and obtain the left hand side of equation (160) as

[TABLE]

If we ignore the off-diagonal part of the density matrix and consider the linearized version of equation (187), we obtain

[TABLE]

Up to linear order in perturbation, we replace $g^{<}_{\lambda\lambda}(\mathbf{p},\omega,\mathbf{R};\{\mathbf{u}_{\lambda}\})$ by the equilibrium one

[TABLE]

Again, we integrate over $\omega$ to obtain equal time Green’s function, the equation (188) becomes

[TABLE]

We turn on the shear which, by definition, is the divergence-free background flow $\partial_{i}u_{\lambda}^{i}=0$. We define the shear tensor

[TABLE]

We can turn on the space-time independent off-diagonal part of the shear tensor so that

[TABLE]

where $F_{\lambda}$ is space-time independent. The equation (188) becomes

[TABLE]

We need to solve the kinetic equations and calculate equal time Green’s function $g^{<}_{\lambda}(\mathbf{p};\{F_{\lambda}\})$, and derive the stress tensor $T_{\lambda}^{ij}$. The viscosity can be read off from the equation

[TABLE]

C.4 The collision integral induced by Coulomb interaction

In this section, we will derive the right-hand side of (160) in the Weyl-Wigner formulation. Let’s simplify the collisional part of the kinetic equation. With the assumption that, in Weyl-Wigner coordinates, $G^{\lessgtr}(\mathbf{r},t,\mathbf{R},T)$ and $\Sigma^{\lessgtr}(\mathbf{r},t,\mathbf{R},T)$ vary slowly in $\mathbf{R}$ and $T$, we can perform the Fourier transformation and rewrite the collision integral. In the rest of this section, we will omit $X$,$\mathbf{R}$ and $T$ to simplify the notation. We define the following notation

[TABLE]

The explicit formula of self energy is given by

[TABLE]

Since flavor index is conserved at each vertex, in the following formulae, we will omit the flavor index. The transformation of the Green’s function between sub-lattice index and band index is given by (174) and (175). We can rewrite the self energy as

[TABLE]

where $\alpha$,$\beta$,$\gamma$ and $\delta$ now are sub-lattice indices only. The factor of $N_{f}$ in the first term comes from the summation over the flavor index of the loop diagram. The simplification comes from the assumption that the perturbed Green’s functions for each flavor are the same. Now we consider only the diagonal part of the density matrix in band index

[TABLE]

The right hand side of (160) for each band index $\lambda$, after integrating over $\omega$ and multiplying by $\mathcal{U}^{\dagger}(\mathbf{p})$ on the left and $\mathcal{U}(\mathbf{p})$ on the right, is the collision integral for this band index and given by

[TABLE]

where

[TABLE]

and

[TABLE]

Combining equations (200), (201) and (202) gives us the contribution to collision integral from the Coulomb interaction between quasi-particles (70).

Appendix D Detailed derivation of (76)

Fourier transforming (74) and writing this in terms of the creation and annihilation operators, we find

[TABLE]

where

[TABLE]

From the interacting Hamiltonian (203), we work out the matrix element

[TABLE]

We square the matrix element and average over disorder realizations with $N_{\textrm{imp}}$ impurities

[TABLE]

From Fermi’s Golden Rule, the scattering rate is

[TABLE]

where the density of states in 2d is

[TABLE]

so

[TABLE]

The collision integral for the impurity scattering is

[TABLE]

[TABLE]

Now write

[TABLE]

and linearize the collision integral (in this case the exact collision integral is already linear in $h$).

[TABLE]

Now let us take the fully screened potential

[TABLE]

to obtain

[TABLE]

Using

[TABLE]

leads to

[TABLE]

where we have defined

[TABLE]

and

[TABLE]

[TABLE]

which gives us equation (76).

Appendix E Detailed numerical calculations

In this section, we explain in detail our numerical setup as well as the steps of calculations. In particular, we explain how to turn the QBE into a matrix equation. We follow [14]. The Boltzmann equation is

[TABLE]

The suggested ansatz in this calculation is

[TABLE]

Expand in terms of basis functions

[TABLE]

such that $a$ is dimensionless. Here the basis functions are taken to be

[TABLE]

where $K=\sqrt{\beta/m}k$ is the dimensionless momentum. For all powers $n>2$ we multiply by an exponential factor so the basis function is $K^{n}e^{-K/2}$. We expand in up to 16 basis functions. Increasing the number of basis function changes the results only marginally. Use the fact that this must be valid for all $\mathbf{E}$, sum over $\lambda$, multiply separately by $\mathbf{\hat{p}}g_{m}(\lambda,p)$ and $(\mathbf{\hat{p}}\times\mathbf{\hat{z}})g_{m}(\lambda,p)$ and integrate over $\mathbf{p}$. This yields two equations that can be summarized in matrix form as

[TABLE]

where we defined the dimensionless matrices

[TABLE]

and

[TABLE]

and the dimensionless vector

[TABLE]

(225) can be inverted to yield

[TABLE]

where

[TABLE]

The charge current is

[TABLE]

The DC conductivity is read off as

[TABLE]

[TABLE]

where we have exceptionally restored $\hbar$ and where the dimensionless vector

[TABLE]