Time-reversal odd transport in bilayer graphene: Hall conductivity and Hall viscosity

We compute the time-reversal odd effective action for the low-energy model of bilayer graphene in the presence of out-of-plane magnetic field. A generating functional for the effective action that captures the electromagnetic response to all orders in momentum and frequency is presented and evaluated to the third order in space-time gradient $\mathcal O(\partial^3)$. In addition, we calculate the Hall viscosity and derive an explicit relationship with the $q^2$ coefficient of the Hall conductivity. It is reminiscent of the Hoyos-Son relation in the Galilean invariant systems, which can be recovered in the limit of large filling factor $N$ .


I. INTRODUCTION
The family of graphite multilayer is one of the most intriguing in the realm of modern condensed matter. The simplest model in the group, graphene, has drawn enormous attention from both theoretical and experimental communities 1 . Its linear dispersion at low energy makes it a low dimension example of particle-hole symmetric ultra-relativistic fermion. Its unconventional electronic property distinguishes itself from other semiconductors made up of ordinary non-relativistic quasi-particles. When placed in a magnetic field, it possesses the anomalous Hall conductivity 1 2 e 2 h at the filling factor ν = 0 2 , which deeply relates to the concept of parity anomaly for a single Dirac cone defined in (2+1) dimensions 3 .
Even more fruitful physics emerges as one stacks two layers of graphene on top of each other. In the AB-stacked configuration as shown in Fig.2, the low-energy projection of the model forms another family without any relativistic analogue: a particle-hole symmetric two-band semiconductor with parabolic dispersion 4 . It also serves as a model possessing a Fermi surface with a Berry phase of 2π as well as a candidate for the dual description of ν = 1 fractional quantum Hall state in the context of fermion-vortex duality 5 .
As a background magnetic field is turned on, Landau levels form (Fig.3) in bilayer graphene as well. The spectrum, nevertheless, possesses two zero-energy bands 6 . This feature reshapes our understanding of low energy physics in quantum Hall systems 7 . In particular, in large magnetic field, the conventional wisdom and machinery of lowest Landau level projection has to be modified since the low-energy sector contains more than holomorphic wave functions in symmetric gauge and the system can host exotic quantum phases in the lowest Landau level 8 . The dielectric property and the low energy excitations in the zero-energy bands under this circumstance were investigated in Ref. 9 and Ref. 10.
This work concentrates on the time-reversal odd responses of bilayer graphene to dynamical and inhomogeneous external perturbations in unique in (2+1) dimensions, which is partially addressed recently in an independent work Ref. 11. Here we adopt a different approach that systematically generates the effective action as a functional of external gauge field A µ to all orders in momentum and frequency, using which we compute the effective action as in gradient expansion to explore the large scale dynamics. As a result, we will show for the low-energy model for bilayer graphene in a background mag-netic field, the Hall conductivity in an inhomogeneous and dynamical electric field E i (ω, q) at filling factor N is where and ω c are the magnetic length and cyclotron frequency. We shall also establish a relationship between Hall viscosity and the coefficient of (q ) 2 term. In particular, with a high energy Landau level cutoff N c , the static ω = 0 result reduces to It is reminiscent of the Hoyos-Son relation 12 for Galilean invariant systems. An explicit formula is derived to compute higher order corrections in N −1 and N −1 c . This paper is organized as follows. In Sec.II we review the calculation of the one-loop effective action from a microscopic model and introduce the low-energy two-band model for bilayer graphene. The Feynman rules are specified and a generating functional for polarization tensor is derived. The time-reversal odd effective action is computed to cubic order in space-time gradients. We then focus of the coefficient of Hall conductivity at the order of (q ) 2 and investigate its relationship with the Hall viscosity. In Sec.III, we construct the stress tensor and compute Hall viscosity and orbital magnetic susceptibility for our model. This provides numerical support for the observation established in Sec.II. Finally we revisit the conductivity tensor using the Kubo formula in Sec.IV and derive an exact algebraic relation that connects the Hall conductivity and Hall viscosity in the absence of space-time symmetry. We then conclude the paper. Details of computation are given in the appendix for comprehensiveness.

A. Methodology
The information of electromagnetic response can be compactly summarized in the effective action as a functional of U(1) gauge potential. Starting with a microscopic model for matter field ψ defined by the action S[ψ], one gauges the charge symmetry by coupling the charge and current densities j µ = (ρ, j i ) with the gauge field A µ . The effective action arXiv:2003.07870v1 [cond-mat.mes-hall] 17 Mar 2020 is defined as follows.
If the external electric and magnetic fields E and B take constant value, Eq. (3) can be computed and serves as an example of Euler-Heisenberg effective action 13 . As far as the linear response is concerned, S eff is usually expanded as a multinomial in A µ . In the d-dimensional Fourier space, the effective action under Gaussian approximation assumes the form wherej denotes the average charge in the ground state. The polarization tensor Π µν (q) encodes the response functions. In (2+1) dimensions, gauge symmetry alone fixes the form of effective Lagrangian L eff . Organized in number of derivatives, the Lagrangian reads Computing the functional determinant for a specific model yields the parameters k, , µ, α, and β. For a fermionic system, computing Eq. (3) amounts to evaluating the functional determinant of the fermion action. Formally one can decompose the microscopic Lagrangian for the fermion ψ into the free part iψ † D −1 ψ and the potential part −v(x)ψ † ψ. Such decomposition is straightforward for a free fermion system in an external potential. If a two-particle interaction is present, one can first introduce an auxiliary field by the Hubbard-Stratonovich transformation to decompose the two-particle term. This way the fermion part of the theory can be integrated in the path integral, resulting in a functional of v. Expanding the functional to the second order of v, it reads This formula has an intuitive diagrammatic representation shown in Fig.1. The inverse of action D corresponds to the Feynman propagator of the free theory and the potential v is the vertex. One can systematically compute Eq. (6) using the standard perturbation methods in quantum field theories.

B. The model for bilayer graphene
Let us now apply the above machinery to the low energy model of AB-stacked bilayer graphene 6 . We depict the lattice structure in Fig.2. The minimal low energy model for each valley in the Brillouin zone contains two copies of Dirac fermions and a hopping amplitude 2m bridging sites B and A. In the low energy regime ω m , the four-band model can be projected to the dominant two bands. The model for valley K is where π = π x − iπ y and π i = p i + A i is the kinematic momentum. ψ K is the spinor storing the dominant two bands. For this model, k and the dielectric constant were discovered in Ref. 9. The focus of this paper is to compute α and β. We first derive the propagator D for model (7) in a magnetic field. Turning on a finite background magnetic field in symmetric gauge (Ā x ,Ā y ) = B 2 (−y, x), the spectrum of the system consists of non-trivial Landau levels. It also introduces the length scale of magnetic length = B −1/2 . To solve the spectrum, it is convenient to introduce the ladder operators which satisfies [a, a † ] = 1. The Hilbert space of ( r, p) operators is then organized partly using the eigenstates of the operator n = a † a, {|n |n ∈ Z + }. The Hamiltonian (7) then can be shown to have the eigenstates 15 associated with the spectrum with the cyclotron frequency ω c = B/m . The degeneracy of each Landau level, as in 2DEG in symmetric gauge, is encoded the other set of ladder operators [b, b † ] = 1 which generates the angular momentum z = 1 2 ( r × p − p ×r) z . (See Ref. 16.) Note that in this system we no longer treat z as the canonical angular momentum. The algebraic structure of the second pair of commuting ladder operators still holds. Denoting the Hilbert space with 1 . (11) Next, we turn on the perturbation on top ofĀ µ →Ā µ + A µ , leading to the variation of Hamiltonian with the momentum operators In order to perform gradient expansion, the Fourier transform needs to be introduced properly since the coordinates x = (x, y) are treated as operators on Hilbert space. Given a cvalued vector k, we recognize that where the dimensionless complex momenta are k = √ 2 (k x − ik y ) andk = k * . Each operator in v(t, x) is Fourier-expanded as an integral of these exponential operators weighted by the Fourier coefficients. For concrete instances, we Fourier transform terms linear in external field.
One could notice that because of the ±i involved in the definitions of (a, a † ), perturbation in x direction now resembles to σ y and vice versa. This twist will lead to an extra minus in the effective action. We have equipped the machinery for computing the traces in Eq. (6). As far as the response functions are concerned, the second diagram in Fig.1 is the only non-trivial one. The first diagram gives rise to the ground-state charge density, while the third one, the contact term or diamagnetic current, vanishes exactly for this model. For vertices v in Eqs. (15a), (15b), (15c), the trace assumes the general form The denominator does not depend on the angular momenta m and m . Thus the trace over {|m , |m } space can be computed separately. To this end, we decompose the vertices (15a), (15b), and (15c) into This way we trace out m|e ib †k e ibk |m m |e ib †k e ibk |m , obtaining a delta function 2π 2 e |k| 2 δ (2) (k+k ). The time and frequency integrals for (t, t ) and (Ω, Ω ) can be performed with the help of δ function and residue calculus. The non-trivial summands left incorporate only the virtual processes between filled and empty Landau levels (sgn(ξ n )sgn(ξ n ) < 0). Suppose the fermi energy is pinned in between the N th and the (N + 1)th Landau levels as shown in Fig.3. We then have with γ µ nn denoting the matrix element (n|γ µ |n ). Eq. (18) is an exact result of the trace to all orders in frequency and wave numbers. The matrix elements γ µ nn (k) can be written as linear combinations of associated Laguerre polynomials L α n (|k| 2 ). The long-wavelength and low-frequency effective theory is derived upon expanding the Π µν (k) in small ω and |k|. By comparing Eq. (18) and (3), the polarization tensor is identified as

C. Results
To avoid the ambiguity of the doubly degenerate bands, we fill them simultaneously by taking N > 2 and compute Π 0i and Π ij to the any specified limit or desired order in momentum. For instance, with a homogeneous external electric field of form E = E ω e −iωt , Π 00 = Π 0i = 0 and Π ij can be computed exactly to all orders in frequency: .
The denominator implies a resonance in transport at two distinct frequencies In large N , one of them reduces to the cyclotron frequency ω 2 = ω 2 c and the other one diverges linearly in N , i.e. ω ∼ ω c N . The former pole indicates the system respects Kohn's theorem 17 asymptotically in large filling factor. At large scale, the effective action is organized by powers of ( ∂ i ) and ∂ t /ω c . Evaluating the sum to order O(∂ 3 ), the parity-odd part of the Lagrangian is At the first glance, one may be concerned with the precise meaning of the level of Chern-Simons term. For nonrelativistic fermions, the level corresponds precisely to the number of filled Landau level. The level here N is the index of the largest filled Landau level. Owing to the filled negative energy bands, the number of filled Landau levels and thus the interpretation of N could be ambiguous. It is resolved by understanding the net contribution from the Fermi sea. By redoing the computation for N = −1, it is straightforward to confirm negative energy bands contributes to 1 4π A dA. Therefore, N + 1 is understood as the number of filled bands with non-negative energies.
We refer the Hall conductivity σ H to the coefficient of the 18 . To order (q ) 2 and ω 2 , the Hall conductivity is with the DC value σ H (0) = N/(2π). To make more sense out of this result, we can look at the static large N limit, in which the hole band becomes insignificant and we should recover the physics of non-relativistic 2DEG. Taking N → ∞, Eq. (22) can be organized in the form Formally, it produces the Hoyos-Son relation 12,19 for 2DEG: if we identify the number density as N/(2π 2 ).
Let us remind ourselves the Hoyos-Son result establishes a relation between the Hall viscosity η H (ω) at long wavelength, orbital magnetic susceptibility − ∂ 2 ∂B 2 (B) and the coefficient of Hall conductivity at order (q ) 2 based on the Galilean symmetry of the microscopic physics. In addition, the 2DEG result is physically understood as the ratio of η H to charge density n. Starting from model Eq. (7), we would not expect this relation to persist at finite N because it lacks the Galilean symmetry and the charge density depends on the regularization of the bottom of Fermi sea. Nonetheless, motivated by the observation at large N , we would like to explore if a similar or approximate relation exists without any particular spacetime symmetry. In particular, we wish to clarify the roles of Hall viscosity and orbital magnetic susceptibility played in between.

III. HALL VISCOSITY AND ORBITAL MAGNETIC SUSCEPTIBILITY OF BILAYER GRAPHENE
To proceed with the last observation, we wish to compute the Hall viscosity η H and orbital magnetic susceptibility of the model Eq. (7). Within the framework of linear response, η H is given by the correlation function of stress tensors τ xx τ xy , which can be expressed using the Kubo formula as 19,20 where a, b label bands above and below Fermi energy. Defining the stress tensor nevertheless requires some caution for the following reason. The stress tensor can be defined as the response of Hamiltonian with respect to the variation of metric g ij . On the other hand, there is no obvious covariant way of coupling (7) to a general curved manifold. The way we circumvent this conceptual obstacle goes as follows. Instead of the projected low-energy Hamiltonian (7), we revisit the original model consisting of two Dirac spinors, where one can consistently define the action over a curved manifold 21 , and hence the stress tensor τ ij . We derive the stress tensor in this manner and project the components again to the low-energy bands 22 . The components of τ ij are found to be In Appendix B, we present another distinct way of deriving the above results. A natural conjecture for the model on a general curved manifold is postulated and τ ij is obtained by varying the Hamiltonian with respect to the vielbein field. The recent work Ref.11 derives the stress tensor by constructing strain generators and computing their commutators with the model Hamiltonian 23 . This approach does not evade from the difficulty elaborated above. Taking the pseudospin degree of freedom for instance, it is not transparent the pseudospin matrices generate physical rotations and should be part of the strain generators. Plugging τ xx and τ xy into Eq. (24) yields the Hall viscosity to all orders in frequency. Particularly in the static limit ω → 0, we have Caution is required for computing orbital magnetic susceptibility − (B), the negative of the second derivative of the energy density as a function of background field. Suppose we naively evaluate the energy density per Landau level (B).
After canceling summands from n = −N to n = N , the remaining sum is proportional to n(n − 1), which is severely divergent. To extract a finite result, we regularize the sum by introducing a natural cutoff N c obeying ω c N c (N c − 1) = m , which constrains the validity of the low energy model. The regularized the sum now reads The sum can be evaluated approximately with the Euler-Maclaurin formula. In a double expansion of large N and small ω c /m , the leading contribution is Since N is bounded from above by N c ∼ m /ω c , at large N the logarithmic term is sub-leading. Consequently, both Eq. (26) and Eq. (28) coincide with 2DEG results in the limit N → ∞ 24 . In the rest of the article, we exploit these observations and establish a concrete algebraic identity.

IV. AN ALGEBRAIC RELATION FROM KUBO FORMULA
To better understand the relationship, we reexamine the Hall conductivity from the more conventional perspective of the Kubo formula. It is algebraically equivalent to the computation of the two-point functions Π µν . Nonetheless, it can make our speculation more transparent. The current operator in the first quantized form is given by the symmetrized velocity operator Applying the Kubo formula for Hall conductivity to the current (29), we obtain This formula can be quickly justified by looking at 1 iω Π ij , which reduces to conductivity tensor in temporal gauge A 0 = 0. The above formula can be straightforwardly expanded in small q = (q ). To proceed, let us exploit rotational invariance and take q = (0, q).
To organize the products, we first observe the current operators map the n-th Landau level to (n ± 1)-th, whereas the stress tensor operators only generate transitions between n → n, n ± 2. Their products thus has no finite summand and as a result, there is no term linear in q. Another non-trivial observation is at the level of matrix elements (τ xx ) ba (τ xy ) ab = (τ xx ) ba (iσ z τ yy ) ab even though τ xy = (iσ z τ yy ) ab in general. Consequently, expanding the numerator of (30) to the second order in (q ) yields the following identity: which unambiguously identifies the role of Hall viscosity as part of the coefficient of (q ) 2 . The sum in the second line of Eq. (31) is not expressed directly with a physical observable at finite ω. In the static limit, it sums up to − N 2 2π = − (B) − 1 16 ln N Nc , and it reduces to Eq. (2) in the static limit ω → 0. At finite frequency, despite its lack of lucid physical interpretation, Eq. (31) provides with a decomposition that generates corrections to Hoyos-Son relation in powers of ω 2 and N −1 .

V. DISCUSSION AND CONCLUSION
The model Eq. (7) is often considered a hybridization of Dirac and non-relativistic fermions, capturing the particlehole structure of the former and the massive parabolic dispersion of the latter. The results derived in the main text entail in many manners that the features of non-relativistic fermion, or the restoration of Galilean symmetry, manifest asymptotically in the limit of large filling factor. Explicit examples are the forms of Hall conductivity, Hall viscosity, orbital magnetic susceptibility, poles of transport coefficients and the Hoyos-Son relation, although, in term of hydrodynamic relation, it is not yet clear in what sense the charge current and momentum density approach each other in the same limit. Established exact formulae can be utilized for interpolating between the asymmetric model and the Galilean symmetric paradigm.
To move forward, this work permits various open directions. Separate works under progress include the real-time effective theory at finite temperature using the Schwinger-Keldysh formalism 25 , generalization of the linear response theory to graphite multilayers as in Ref. 26, and interaction generated transport properties owing to either an instantaneous Coulomb interaction or a mixed-dimensional Maxwell term 27,28 . Equally interesting are the inclusion of lattice effect into the generalized Hoyos-Son relation 29 and clarifying the distinction between dual descriptions of ν = 1 fractional quantum Hall state 5,30 .
To conclude, we determined the time-reversal odd electromagnetic response for the low-energy model for bilayer graphene to quadratic order in momentum and frequency and established a conductivity-viscosity relationship in the absence of obvious space-time symmetry. We investigate the limit in which the symmetry is restored and provided supports from concrete computation. In addition to the conclusions in the above, the effective action derived in the main text can be further applied to exploring unknown faces of this model.

ACKNOWLEDGMENTS
The author thanks M. Lapa and F Setiawan for comments on the manuscript. This work is supported, in part, by U.S. DOE grant No. DE-FG02-13ER41958 and a Simons Investigator Grant from the Simons Foundation. Additional support was provided by the Chicago MRSEC, which is funded by NSF through grant DMR-1420709. Here we document the matrix elements for the vertices used in the main text. The fundamental ingredient is n|e ia † k e iak |n . Let us evaluate The sum can be expressed in terms of generalized Laguerre polynomials L m n (x). If n < n , we have to evaluate s at s = s + n − n, On the other hand if n ≥ n , we have to evaluate at s = s + n − n n|e ia † k e iak |n = (ik) n−n n ! n! L n−n n (|k| 2 ).
Together with eigenstates (9a), and (9b), it is then straightforward to compute γ µ nn (k) and γ µ n n (k) = (γ µ nn (−k)) * . In the following, we present some results. For definitely we consider the case n > 1 and n < n.
a. γ 0 nn (k) For n = 0 and 1, For |n | > 1, b. γ i nn (k) It is convenient to define in terms of which, For n > 1, Using these matrix elements, the generating functional (18) can be computed in a straightforward manner.

Appendix B: Another derivation of stress tensor
We derive here the stress tensor with a conjectured curved space generalization of model (7) 31 . Let us consider a curved space endowed with metric g ij (x). A matrix field vielbein e a j (x) can be introduced via the relation g ij = δ ab e a i e b j , where a, b = 1, 2 are local SO(2) rotation indices 32 . The inverse field E i a fulfills E i a e b i = δ b a and g ij = δ ab E i a E j b . With the vielbeins, the spatial gradient operator ∂ a can be promoted to curved manifold following ∂ a → E i a ∂ i . Let us further denote E i ± = E i 1 ± iE i 2 . The natural generalization of Hamiltonian (7) in the curved space reads where ψ and χ are the Grassmannian density of φ A and φB and h.c. refers to Hermitian conjugate. To derive the stress tensor relevant for viscosity computation, we turn on a slightly curved manifold with the metric g ij = δ ij + δg ij and assume δg ij to be homogeneous. Under this circumstance, τ ij assumes the following form (B2) Applying this formula to (B1), we arrive at the identical results, Eqs. (25a), (25b), and (25c).