Theory of Conductivity of Chiral Particles

Theory of Conductivity of Chiral Particles

Janik Kailasvuori International Institute of Physics, Universidade Federal do Rio Grande do Norte, 59078-400 Natal-RN, Brazil Max-Planck-Institut für Physik komplexer Systeme, Nöthnitzer Str. 38, 01187 Dresden, Germany    Břetislav Šopík Central European Institute of Technology, Masaryk University, Kamenice 735, 62500 Brno, Czech Republic    Maxim Trushin University of Konstanz, Fachbereich Physik M703, 78457 Konstanz, Germany
July 15, 2019
Abstract

In this methodology focused paper we scrutinize the application of the band-coherent Boltzmann equation approach to calculating the conductivity of chiral particles. As the ideal testing ground we use the two-band kinetic Hamiltonian with an -fold chiral twist that arise in a low-energy description of charge carriers in rhombohedrally stacked multilayer graphene. To understand the role of chirality in the conductivity of such particles we also consider the artificial model with the chiral winding number decoupled from the power of the dispersion. We first utilize the approximate but analytically solvable band-coherent Boltzmann approach including the ill-understood principal value terms that are a byproduct of several quantum-many body theory derivations of Boltzmann collision integrals. Further on, we employ the finite-size Kubo formula with the exact diagonalization of the total Hamiltonian perturbed by disorder. Finally, we compare several choices of Ansatz in the derivation of the Boltzmann equation according to the qualitative agreement between the Boltzmann and Kubo conductivities. We find that the best agreement can be reached in the approach where the principle value terms in the collision integral are absent.

I Introduction

The isolation of the one-atomic carbon layers Novoselov et al. (2004) has boosted an extremely intensive research in the field of quantum transport associated with Dirac electrons.Castro Neto et al. (2009) One of the primary reasons of such an overwhelmed interest is the bizarre chiral nature of the Dirac carriers which makes the conduction and valence states entangled with each other. This coupling leads to many nontrivial effects,Castro Neto et al. (2009) including the famous Klein tunneling Katsnelson et al. (2006).

The chiral carriers are usually characterized by means of the pseudospin which is an additional quantum number formally similar to the real spin in spin-orbit coupled systems. In a single layer graphene, the pseudospin keeps its radial orientation along the Fermi circle resulting in the winding number . The winding number is directly related to the Berry phase responsible for the anomalous quantum Hall effect observed at the very beginning of the graphene rush.Novoselov et al. (2005); Zhang et al. (2005) There are quite a few speculations on the important role of the pseudospin-coherent contribution in the electrical dc conductivity of Dirac fermions at low carrier concentrations. Auslender and Katsnelson (2007); Trushin and Schliemann (2007); Kailasvuori and Lüffe (2010); Culcer and Winkler (2008); Liu et al. (2008) There is however no agreement whether this contribution really matters and can be extracted from the total conductivity.

There are two aspects of this problem worthy of discussion. The first one is experimental: It is very difficult experimentally to detect the pseudospin-coherent contribution because the subtle quantum mechanical effects responsible for this contribution are easily smeared out by temperature, ripples, and charged inhomogeneities inherited in graphene. In spite of the impressive progressMayorov et al. (2012) made with suspended graphene, this contribution has not been extracted yet. The second problem is purely theoretical: It is still not clear how large the pseudospin contribution should be since the numerical methods can give only an estimation Nomura and MacDonald (2007) while the approximated Boltzmann-like approaches Auslender and Katsnelson (2007); Trushin and Schliemann (2007); Kailasvuori and Lüffe (2010); Culcer and Winkler (2008); Liu et al. (2008) are unreliable at low carrier concentrations and may lead to different answers depending on the approximations involved.

The main goal of the present work is to compare the exact numerical (finite-size Kubo) and approximated analytical (Boltzmann-like) pseudospin-coherent conductivities and in that way to figure out the correct approximation needed to derive the most reasonable analytical expression. In particular, we focus on the dependence of pseudospin-coherent conductivity on the winding number . Surprisingly, the semiclassical pseudospin-coherent conductivity acquires qualitatively different behavior on depending on the approximation done. The subsequent comparison with the numerical calculation helps us to rule out the “wrong” approximations and keep the “correct” ones.

Note, that the chiral particles with an arbitrary winding number can be realized in the chirally stacked multilayer graphene. The recent investigations of the band structure in a few layer graphene Latil and Henrard (2006); Guinea et al. (2006); Avetisyan et al. (2010) indicate its strong sensitivity to the stacking pattern. Graphene layers placed together in the most natural way do not lie exactly one on top of the other but are shifted in such a way that only half of the carbon atoms have a neighbor in another layer and the other half are projected right into the middle of the graphene’s honeycomb cell. If the third layer aligns with the first (and the layer with the -th) then we arrive at the most abundantly occurring arrangement of graphene layers known as Bernal (or ABAB…) stacking. An alternative stacking is obtained when the third layer aligns neither with the first nor with the second but is shifted with respect to both, and it is the fourth layer that again aligns with the first. This arrangement is known as rhombohedral (or ABCABC…) stacking. The simplest model used in the study of the latter system Min and MacDonald (2008) is the model that we shall use as our testing ground in the investigation on theoretical methods to calculate the conductivity.

The beauty of the effective low-energy Hamiltonian describing the carriers in the chirally stacked multilayer graphene has attracted many theoreticians working in the field of quantum transport. Among the most recent works are the theory of chirality-dependent phonon scattering Min et al. (2011a), the semiclassical Boltzmann transport theory Min et al. (2011b), the investigations of Landau levelsKoshino and McCann (2011), Berry phase Koshino and McCann (2009), effects of screening Koshino (2010); Min et al. (2012), optical conductivityMin and MacDonald (2009), and THz absorption Trushin and Schliemann (2012). We would like to warn the reader, that although our ultimate interest is to understand conductivity in graphene type of systems, the aim of this paper is purely methodological, and we undertake no discussion about the connection to realistic systems and to experimental findings. For example, we will for simplicity and analytical solvability assume a fully homogenous background, whereas realistic graphene might have to be modeled with a variable potential to account for ripples and other sources of inhomogeneity.

The outline of the paper is the following. In section II we discuss some preliminaries that will be used in later sections. Section IV offers in more detail some of the background and motivation of the present paper. In section III we present an analytical heuristic derivation of the conductivity by evaluating the Kubo formula using the chiral plane-wave states of free electrons of multilayer graphene. The derivation, interesting in its own right, will serve among other as a easy introduction to band-coherent effects and to the distinction between intra-band and inter-band contributions to the conductivity. In section V we will discuss the conductivity derived with a band-coherent Boltzmann equation approach. In section VI we present the results of the numerical studies of the Kubo formula. In section VII we wrap up with a discussion.

Ii Preliminaries

Many properties of bilayer graphene can be well understood in terms of two independent, time-reversal related copies of the effective two-band HamiltonianMcCann and Fal’ko (2006)

 H0=12m∗(0(px−ipy)2(px+ipy)20). (1)

This hamiltonian exhibits a coupling between momentum and a spin type of variable, a pseudo-spin that in this case derives from the non-equivalent lattice sites, one from each layer. For -layer rhombohedrally stacked graphene one can with similar approximations derive an effective low-energy hamiltonian that also turns also out to be a two band model, this time given by Min and MacDonald (2008)

 H0=γ(0(px−ipy)N(px+ipy)N0), (2)

with (like above) being determined by material parameters of the system. As a realistic models for rhombohedral multilayer graphene the hamiltonian (2) is expected to become increasingly poor with an increasing number of layers. We wish to exploit the hamiltonian (2) for theoretical and methodological studies because of its combination of simplicity and yet very non-trivial chiral properties. For our studies it will turn out to be interesting to generalize the hamiltonian (2) to the more general although artificial hamiltonian

 H0=\boldmath{b}⋅\boldmath{σ% }=b(k)\boldmath{^b}(\boldmath{^k})⋅% \boldmath{σ}=γ(ℏk)Nd(0e−iNcθeiNcθ0), (3)

where is the absolute value of the wave number , is the angle of the direction of the wave number and is the (pseudo)spin-orbit field. In the hamiltonian (3) the power of in the dispersion has been decoupled from the winding number of the direction of the pseudospin-orbit field as a function of , which determines the chiral properties of the hamiltonian. In the special case this hamiltonian reduces to the hamiltonian (2). The benefit of introducing such an artificial Hamiltonian is that one can take advantage of some special features of a quadratic dispersion, to be explained later, but still keep the chiral properties general.

For both the numerical and analytical treatments of the conductivity we will need the velocity matrix, defined as . For the multilayer Hamiltonian (2) one has in the basis of the sublattice orbitals (that is, in the basis in which equation (2) is written) the -component

 vx=γℏN−1(0(kx−iky)N−1(kx+iky)N−10). (4)

For us it will sometimes be convenient to express matrix quantities in the eigenbasis of the Hamiltonian, which we also call the chiral basis. For the Hamiltonian (3) the eigenstates are given by

 ψs\boldmath{k}=ei% \boldmath{k}⋅\boldmath{x}L(1seiNcθ) (5)

where is the pseudospin index and where is the angle of the momentum vector in the momentum plane. The eigenenergies are . The velocity operator in the chiral basis reads

 (6)

where and . The velocity operator can formally be decomposed into two parts. The intraband part accounts for the ordinary velocity components for each band as where the bands not talking to each other. In the chiral basis it is diagonal. On the other hand, the interband part couples the two bands and is an off-diagonal matrix in the chiral basis. The interpretation of these matrix components is less intuitive. They reflect the Zitterbewegung property of Dirac electrons. Both intraband and interband components of the velocity can contribute to the current. Within the Kubo formula picture we will see this in next section. In the Boltzmann picture, we note that for a matrix-valued distribution function the current in the -direction is given by

 jx=∫d2k(2π)2Tr(vxf)=2∫\boldmath{k}vk(\boldmath{^k}xf\boldmath{^b}+\boldmath{^θ}xf\boldmath{^c}NcNd) (7)

in terms of the components of the matrix-valued distribution function

 f=f01+\boldmath{f}⋅% \boldmath{σ}=f01+(f\boldmath{^b}% \boldmath{^b}+f\boldmath{^c}\boldmath{^c% }+fz\boldmath{^z})⋅\boldmath{σ}=(f0+f\boldmath{^b}fz−if% \boldmath{^c}fz+if\boldmath{^c}f0−f\boldmath{^b})ch (8)

with the vector part or spin part of decomposed along the -dependent basis vectors with . The densities (in phase space) in the upper and lower band, respectively, we find in the chiral basis on the diagonal. Note however, that for the type of Hamiltonians we consider, only the spin component but not the total charge density contributes to the current.

The conductivities (the driving electric field is assumed to be in the -direction) will in the present paper be given as functions of the dimensionless variable

 ℓkF=τtrvFkF=NdτtrϵFℏ. (9)

The momentum relaxation time is derived from integral

 I(k)=∫d2k′(2π)2δ(ϵk−ϵk′)W\boldmath{k}% \boldmath{k′}cos2NcΔθ2(1−cosΔθ) (10)

where gives the collision probability rate according to Fermi’s golden rule for spinless particles scattered by an impurity potential and the unusual trigonometric factor comes from the spin overlap between the chiral eigenstates at a momentum angle difference . In particular, the Drude conductivity is given by

 σD=e2h⋅ℓkF2, (11)

independently of the dispersion and winding number.

The Kubo formula for the conductivity reads

 σKubo=−iℏe2L2∑n,n′fFD(ϵn)−fFD(ϵn′)ϵn−ϵn′⟨n|vx|n′⟩⟨n′|vx|n⟩ϵn−ϵn′+iη, (12)

where are the eigenstates of the total Hamiltonian including impurities but excluding the driving electric potential and is the Fermi-Dirac distribution. The parameter describes the broadening of the eigenstates. For a pure sample it would be infinitessimal. In our case it is finite, as we describe a system with impurities. We will find it useful to study the quantities

 σintra := σKubo(vx→vintrax) σinter := σKubo(vx→vinterx). (13)

Numerical checks will show that there is at the most a negligible contribution from the remaining possible terms containing matrix elements or   .

It is rather clear that the Drude conductivity as a pure intraband contribution resides completely within . However, the above definitions of and are in general not a clean decomposition of intraband and interband contributions to the conductivity, wherefore might on top of contain some interband contributions. For the problem we study in the present paper the Boltzmann approach gives support for this decomposition faithfully separating pure intraband contributions form interband contributions, but in other cases, for example for the much more complex case of monolayer graphene (), the Boltzmann approach (see ref. Kailasvuori and Lüffe, 2010) comes with the nonintuitive message that there can be important spin-coherence originated contributions also in the part of the distribution function that contributes to the current through , although at a first glance seems to be a pure intraband contribution as it describes the polarization along the quantization axis of the Hamiltonian.

Iii Heuristic evaluation of the Kubo formula using the plane wave basis

In the analytical derivation based on the Kubo formula (12) presented in this section we will assume that contains all intraband contributions, in this case only the Drude conductivity , whereas contains all interband contributions. (As mentioned, support for this assumption will come from the Boltzmann approach.) We therefore choose such that we obtain and then derive with the given value of , hoping that we obtain a sensible result. This is one of the assumptions make the derivation a heuristic one.

What furthermore allows an analytical treatment but also contributes to the heuristic nature of the derivation is that we assume that we—rather than using the eigenstates of the full problem including impurities—can use the the eigenstates (5) of the free Hamiltonian, that is the plane wave states, which makes an analytic evaluation simple. At high densities, that is, at high Fermi momenta, the use of plane wave states should be a good approximation as the weak impurities deform the eigenstates only slightly from the plane wave states. However, at small Fermi momenta, close to the Dirac point, one would expect the eigenstates to be substantially deformed by the impurities and we know of no solid argument for why the use of plane wave states would give sensible results also close to the Dirac point. Nonetheless, the outcome appears to be sensible. In any case, the derivation gives a nice introduction into the non-intuitive interband contributions to the conductivity.

The calculation is done in the chiral basis. For the evaluation we use that

 (s,s′)⟨% \boldmath{k}s|vx|\boldmath{k′}s′⟩⟨\boldmath{k′}s′|vx|% \boldmath{k}s⟩(+1,+1)v2kδ\boldmath{k}\boldmath{k′}δ\boldmath{k′}\boldmath{k}(cosθ)2intraband(−1,−1)v2kδ\boldmath{k}\boldmath{k′}% δ\boldmath{k′}\boldmath{k}(−cosθ)2intraband(+1,−1)v2kδ\boldmath{k}\boldmath{k′}% δ\boldmath{k′}\boldmath{k}(−iNcNdsinθ)(iNcNdsinθ)interband(−1,+1)v2kδ\boldmath{k}\boldmath{k′}% δ\boldmath{k′}\boldmath{k}(iNcNdsinθ)(−iNcNdsinθ)interband. (14)

Thus, we can compactly summarize this as

 σ=−iℏe2L2∑ss′%\boldmath$k$v2fFD(ϵsk)−fFD(ϵs′k)ϵsk−ϵs′k(NcNd)1−ss′(1+ss′cos2θ)/2ϵsk−ϵs′k+iη. (15)

For the intraband part we find the zero in the denominator. However, also the numerator contains a zero in . Since the arguments of the two Fermi-Dirac functions are ’close’ to each other, we can Taylor expand the difference as

 fFD(ϵsk)−fFD(ϵs′k)=0+(ϵsk−ϵs′k)∂fFD(ϵsk)∂ϵsk+… (16)

We see that the numerator contains the same zero as the denominator and we can cancel the two occurrences of to obtain the intralayer contribution

 σintra = −iℏe2L2∑s\boldmath{k}v2∂fFD(ϵsk)∂ϵskcos2θϵsk−ϵsk+iη=ℏe2ηL2∑s\boldmath{k}v2cos2θδ(ϵF−ϵsk)=ℏe2ηL2∑\boldmath{k}v2cos2θδ(ϵF−ϵ)⟶ (17) ⟶ ℏe2η∫d2k(2π)2v2cos2θδ(ϵF−ϵ)=ℏe2η∫dθ2πcos2θ∫dϵD(ϵ)v2δ(ϵF−ϵ)=ℏe22ηDFv2F=ℏe22hℏηkFvF,

where we used and where we have assumed zero temperature. We set in order for to reproduce the Drude conductivity .

We now turn to the interband contribution. We find

 σinter = −iℏe2L2Nc2Nd2∑s\boldmath{k}v2fFD(ϵsk)−fFD(ϵ−sk)ϵsk−ϵ−sksin2θϵsk−ϵ−sk+iη= (18) = −iℏe2L2Nc2Nd2∑\boldmath{k}v2sin2θ[Θ(ϵF−ϵ)−Θ(ϵF+ϵ)(2ϵ)(2ϵ+iη)+Θ(ϵF+ϵ)−Θ(ϵF−ϵ)(−2ϵ)(−2ϵ+iη)]= = −iℏe2Nc2Nd2∫dθ2πsin2θ∫dϵD(ϵ)v2[Θ(ϵF−ϵ)−Θ(ϵF+ϵ)][1(2ϵ)(2ϵ+iη)−1(−2ϵ)(−2ϵ+iη)]= = iℏe22Nc2Nd2∫∞ϵFdϵD(ϵ)v22ϵNd/4πℏ2−i2η(2ϵ)2+η2=e22hNc2Nd∫∞ϵFdϵη(2ϵ)2+η2=e24hNc2Nd[π2−arctan(2ϵFη)].

With the previously assumed we obtain

 σinter=e24hNc2Nd[π2−arctan(2ϵFτtrℏ)]. (19)

The total conductivity is therefore

 σ=σintra+σinter=e22hℓkF+e24hNc2Nd[π2−arctan(2ℓkFNd)] . (20)

The rewriting

 σinter=e24hNc2Ndarctan(Nd2ℓkF) (21)

with the arcus tangens function taylor expanded in powers shows how an infinite series of terms of pseudospin originated quantum corrections diverging at can sum up to something finite. In the Boltzmann regime the interband part vanishes and we are left with the Drude contribution. In the opposite limit the Drude conductivity vanishes and the interband contribution saturates at the finite value

 σ(ℓkF=0)=e2hπ8Nc2Nd. (22)

Note that the increasing the chiral winding number pushes up the conductivity at the Dirac point, whereas increasing the power of the dispersion pushes it down. For the multilayer Hamiltonian (2) () the total conductivity has a minimum at the Dirac point , and the value of this minimum is proportional to the number of layers . In particular, for the bilayer case this agrees with the result in ref. Trushin et al., 2010 derived analytically with a spin-coherent Boltzmann type approach. In the case the qualitative behavior is different. A minimum occurs at some non-zero whereas close enough to the conductivity shoots up to saturate at a local maximum, with the value given by (22). Examples are shown in Figs. 1 and 2.

Iv Detailed introduction to the background and motivation of the paper

This section complements the introduction with a more technical and detailed account of the background and motivation of the present study as for the methodological issues concerning the use of band-coherent generalizations of Boltzmann equations to study graphene-related systems.

The Boltzmann equation is a kinetic equation that is a widely used starting point for deriving the transport coefficients and other response functions, also for quantum systems, where it is used as a semiclassical approximation. The appeal of a Boltzmann type of approach is that it makes contact with a classically intuitive picture, in contrast with the kind of numerical evaluation discussed later in the paper. Furthermore, many problems can for a ”first feel” be suitably approximated to be amenable to an analytic solution. Both of these features offer great help for intuition. At the same time, the Boltzmann equation can be derived from a full many-body quantum theory, for example with the Kadanoff-Baym equations as the starting point. This means that it in principle should be possible to keep track of all the approximation and their range of validity, and it also means that there is a way to work out more refined and extended equations, when needed. This stands in contrast to a derivation such as the one in section III. However, whereas the Boltzmann approach is a very simple and instructive tool for simple problems such as deriving the Drude conductivity in an ordinary metal where spin is involved in a trivial way, it gets very complicated as soon as one considers more complicated systems such as graphene, where the pseudospin has a very non-trivial coupling with the momentum. Still, there is the possibility, as discussed in this paper, that a Boltzmann approach could make sense even for graphene-like systems and in some cases allow a fully analytical solution. This should be one clear motivation for investigating generalizations of Boltzmann equations in systems with non-trivial spin, in particular with some sort of spin-orbit coupling.

In this paper it is shown that the simple but due to their chirality yet highly non-trivial model Hamiltonians (2) with the generalization (3) provide an ideal testing ground for scrutinizing and screening the kind of band-coherent generalizations of the Boltzmann equation that have been used in refs. Auslender and Katsnelson, 2007; Trushin and Schliemann, 2007; Culcer and Winkler, 2008; Liu et al., 2008; Culcer and Winkler, 2009; Kailasvuori and Lüffe, 2010; Trushin et al., 2010 to discuss the conductivity in regime close to the Dirac point, that is, at the low charge carrier filling where the chemical potential is close to the zero of energy given by the touching of the two bands, the one with pseudo-spin up and the other with pseudo-spin down along the momentum. In this regime it might be relevant to not only consider the energy eigenstates of the free Hamiltonian, which describe the population of respective pseudo-spin band, but also consider electrons in quantum coherent superpositions of these eigenstates. (We will interchangeably use the terms spin-coherent and band-coherent to make it clear that we do not only consider populations of the eigenstates but also quantum superpositions of these.) The mentioned band-coherent Boltzmann equation approaches treat the spatial degrees of freedom semi-classically as in the standard Boltzmann treatment, but keep the treatment of the the band index—that is, of the (pseudo)spin—fully quantum mechanical. Because of the band-coherent terms, which are elusive to classical intuition, the collision integral of these equation cannot sensibly be given by an educated guess, but have to be derived from full quantum many-body theory.

We are interested in to what extent these kinds of generalized Boltzmann approaches are in principal applicable, considering that the usual textbook criterion for a Boltzmann treatment breaks down close to the Dirac point. We are also interested in if this approach is in practice possible, considering that for Dirac systems the collision integrals in these approaches usually become horrendously complicated, in particular if one takes the often neglected principal value terms seriously. The latter turn the static and uniform Boltzmann equation—usually a differential equation problem—into a integro-differential problem. The mathematical challenges that come when solving these can be clearly seen in ref. Auslender and Katsnelson, 2007).

We show that the band-coherent Boltzmann equations considered in this paper are analytically solvable for the Hamiltonians (2) and (3) with the impurity potential given by point-like impurities, and most remarkable the equations remain so with moderate additional complexity even when principal value terms are kept. The analytical form of these results makes the lack of equivalence of the different Boltzmann approaches conspicuous and facilitate a neat matching against other approaches that are considered more reliable close to the Dirac point, in our case a numerical evaluation of the finite-size Kubo formula. In the special case of bilayer graphene with point-like impurities, most of these approaches coincide. In this case the prediction of a band-coherent Boltzmann approach agrees well with the numerical evaluation, as shown in ref Trushin et al. (2010).

The source of this multiplicity among Boltzmann approaches comes from there being a wide variety of derivations of Boltzmann equations from quantum many-body theory, with both different starting points and different approximations on the way. Even if the different approaches are sometimes motivated by different purposes, there are usually overlapping regimes where a comparison should be possible, like our case of weak and dilute impurities (implying that it should be enough to consider the lowest Born approximation and only non-crossed diagrams). In the context of graphene and spin physics the issue of their equivalence, or rather lack of such, was addressed by one of us Kailasvuori and Lüffe (2010). It was shown that although there turns out to be agreement in some standard electronic systems, there is no general equivalence. In particular, for electronic systems where coupling between momentum and some spin type of variable is strong—as is the case for the here studies Hamiltonians with the pseudo-spin to momentum coupled term being the dominant term—at least six feasible and yet different band-coherent generalizations of the Boltzmann equations can be derived, and, as discussed in ref. Kailasvuori and Lüffe, 2010, they all result in different spin-originated quantum corrections to the conductivity. On the other hand, in ordinary textbook systems with all degrees of freedom treated semi-classically, in particular for essentially spinless electrons with quadratic dispersion like in the ordinary metal, these derivations all agree, which is probably one of the reasons for the issue of equivalence of Boltzmann approaches not having received much attention.

This abundance of candidates for the Boltzmann equation comes above all from the non-equilibrium Green’s functions derivations. Here the Kadanoff-Baym (or Keldysh) equation is the starting point. In the derivation of a Boltzmann equation one has to choose whether the real part of the self-energies should be related to relaxation (through the collision integral), as is the fate of the imaginary part, or if the real part instead goes into the renormalisation of the free drift of the semiclassical quasiparticle. The most important source of discrepancy, however, comes from the choice of Ansatz. The reason that one needs the discussed Ansatz is that one needs some approximation for expressing a two-time non-equilibrium Green’s function (other variables are here left implicit) in terms of of the one-time distribution function in order to get a closed equation for latter. The necessity of the Ansatz is illustrated in appendix C in schematic explanation that has been abstracted and simplified from the kind of elements that occur in actual quantum many-body derivations of Boltzmann equations, see for example ref. Kailasvuori and Lüffe, 2010 and references therein.

The Green’s function derivation can be contrasted with the more common and much more accessible quantum Liouville equation -based derivations of master equations (in our case the generalized Boltzmann equations) for the density matrix . Appendix B shows the origin of the Boltzmann collision integral in such a derivation. Here one is not confronted with the issues mentioned above. In the graphene related literature, almost all treatments have been of the Liouville type, see refs. Auslender and Katsnelson, 2007; Trushin and Schliemann, 2007; Culcer and Winkler, 2008, 2009; Trushin et al., 2010, the to us known exceptions being refs. Liu et al., 2008; Kailasvuori and Lüffe, 2010. In ref. Kailasvuori and Lüffe, 2010 we showed that none of the existing Ansatzes in the Green’s function literature gave a collision integral in the band-coherent Boltzmann equation that matched the collision integral we derived by a Liouville equation approach, and we invented yet another Ansatz to provide the missing link. With this we yet increased the number of Boltzmann equations that can be derived with a Green’s function approach.

One might think that a Liouville equation -based derivation is then the correct derivation as one does not seem to be confronted with these elements of choice. Our study does not support this view. Rather, we find that if any of the approaches can be ruled out, then it is the Liouville equation -based one, at least in its standard form.

A central technical issue in the present paper is the role of principal value terms in the collision integral. They arise as byproduct in all derivations of Boltzmann equations mentioned above, also in the Liouville type of derivation. Principal value terms have been discussed in other contexts (under various names), for example the Boltzmann equations with spinless electrons subject to strong interactions, but have in the context of Dirac electrons and spin-orbit coupled electrons simply been left out, with the to us known exceptions being refs. Auslender and Katsnelson, 2007; Kailasvuori and Lüffe, 2010. The usual Boltzmann collision integrals have as typical ingredients a Dirac delta function expressing the classically intuitive conservation of energy in a collision process. The quantum many-body derivations that lead to these expected collision terms will usually also result in similar terms, which instead of the delta functions involve principal values . (Appendix B shows explicitly where they come about). The latter are rather elusive to a semiclassical interpretation. They might be connected with the Zitterbewegung property of Dirac electrons and spin-orbit coupled electrons. We believe that is motivated to address their presence and role seriously, in particular when the distribution describes quantum coherent superpositions of eigenstates.

As mentioned above, the principal value terms make in most cases the Boltzmann equation horrendously complicated to solve. Maybe one reasons—on top of them being unintuitive—for them usually being neglected, most often without discussion. The special attractive feature of the models (2) and (3) together with point-like impurities is that there are non-vanishing principal value terms due to the band-coherence but an analytical solution of the different Boltzmann equations is nonetheless still possible to infinite order in in spin-originated quantum corrections and sums up to a neat result that remains finite even at the Dirac point . Remarkable here is that the presence of principal value terms only moderately but not qualitatively increases the complexity of the solution. Furthermore, it turns out that the principal value terms modify the conductivity in a very non-trivial way. The studied setting is therefore a very interesting testing ground for screening of the different approaches and an addressing of the role and influence of principal value terms as well as the applicability of band-coherent generalizations of the Boltzmann equation.

V The band-coherent Boltzmann treatment

In this section we are going to derive within the Boltzmann approach the pseudospin originated quantum corrections in the considered graphene type systems and in particular with principal value terms included. However, we have refrained from making the paper self-contained in this part, as this would make the paper far too ominous. For a technical understanding of the results familiarity is assumed with the derivation of Boltzmann equations in general as well as with ref. Kailasvuori and Lüffe, 2010, where technical details for the band-coherent case have been worked out in great detail and generality. This section functions as a guide for those readers who want to see how that work can be applied to the present context. Other reads might want to jump to the end of the section were the main results and conclusions have been gathered.

In the Boltzmann approach the particles are described in terms of a distribution function living in classical phase space. In the presence of spin this function obtains a spin index, usually without any semiclassical approximation, making a density matrix in the spin indices. The Boltzmann equation then reads

 i[H,f]+∂tf+12{vi,∂xif}+eEi∂kif=J[f], (23)

with the bracket standing for anti-commutation. The presence of spin leads among other to the presence of the spin-precession term . However, spin-precession has also a classical analogue in the precession of angular momentum. The truly quantum mechanical ingredients are found in the interband components of and in the collision integrals , for example in Pauli blocking factors and in the collision rates.

The Boltzmann equation can be derived from quantum many-body theory, with the main issue being the derivation of the collision integral. One common derivation is to use the quantum Liouville equation as a starting point. A derivation of a collision integral can be found in appendix B. Another starting point is the Kadanoff-Baym equations, from which one can derive

 [i∂t−H,G<] = ΣRG<−G<ΣA+Σ

After Wigner transformation followed by gradient expansion and integration over the energy the left hand side turns into the left hand side of (23), where by definition

 f(\boldmath{k},\boldmath{x},t)=∫dω2πG<(ω,\boldmath{k},% \boldmath{x},t), (25)

An equivalent form of (25) in a non-Wigner transformed coordinates is

 f(x1,x2,t)=G<(x1,t1,x2,t2)|t1=t2=t, (26)

that is, the distribution function is related to the time-diagonal of the lesser Green’s function. The right hand side of eq. (24)—involving self-energy terms which are functionals of the non-equilibrium Green’s function —turns into something that we can interpret as a collision integral. To derive a collision integral in the lowest Born approximation, the dressed Green’s functions and are replaced by their non-interacting analogues and and in the self-energies the series in orders of the interaction is only kept to second order in the interaction and first order in the impurity density. The retarder non-interacting Green’s function is in our case given by

 G0R\boldmath{k}(ω) = ∑s=±S\boldmath{^b}sω+i0+−ϵs\boldmath{k}=∑s=±S\boldmath{^b}s(Pω−ϵs\boldmath{k}−iπδ(ω−ϵs\boldmath{k})) (27)

with the projector

 S\boldmath{^b}s := 12(1+\boldmath{σ}⋅s\boldmath{^b}%\boldmath$^k$). (28)

The advanced Green’s function is obtained by hermitian conjugation.

As mentioned above, the collision integral is the ingredient of the Boltzmann equation that is the most demanding to derive. One reason is that it is non-linear in the distribution functions. It is not surprisingly in the collision integral that we will find differences between the different derivations. We will let the collision integral be divided up as , with containing the terms that contain delta functions and containing the terms that contain principal value terms. In the case of monolayer graphene (Hamiltonian (2) with ) both parts are sensitive to the approach, even for point-like impurities. (See sec. 6 in ref. Kailasvuori and Lüffe, 2010, in particular eqs. (65,66).) For the multilayer Hamiltonians with generalizations studied in the present paper, only the the principal value part will differ between the different approaches.

Ref. Kailasvuori and Lüffe, 2010 uses a nomenclature for the different approaches composed of a prefix and a suffix. The prefix G1 or G2 indicates whether the collision integral derives from all self-energy terms as suggested in equation (24) or only from the imaginary part of the self-energies, namely from the right hand side of (the with (24) equivalent expression)

 [i∂t−H−ReΣR,G<]−[Σ<,ReGR] = i{ImΣR,G<}−i{Σ<,ImGR}. (29)

In the latter case the real part of the self-energies to the left are instead assumed to renormalize the free drift. (As ref. Kailasvuori and Lüffe, 2010 was not able to carry out this renormalization in the spinful case, the statements and results for the approach have to be taken with more caution than those of .) The suffix GKBA (Generalized Kadanoff-Baym Ansatz Lipavský et al. (1986)), AA (”Alternative Ansatz” or ”Anti-ordered Ansatz” Kailasvuori and Lüffe (2010)) and SKBA (Symmetrized Kadanoff-Baym Ansatz Kailasvuori and Lüffe (2010)) refers to which Ansatz has been used. In the Green’s function derivations considered in this paper the issue of Ansatz is the central one.

The Ansatz gives an explicit relation between the non-equilibrium Green’s function and the distribution function . By definition one has , or in the Wigner-transformed language equation (25). This implicit relation is enough for the left hand side of (24) to turn all Green’s functions into distribution functions and this thanks to the left hand side being linear in Green’s functions. However, the right hand side is non-linear in Green’s functions, wherefore the implicit relation such as (25) is not enough, as explained in appendix C. Here one needs an explicit relation between and (and one needs explicit expressions for and , here given by (27)). Such a relation between a semiclassical object and a fully quantum mechanical object has of course to be an approximation. The Ansatz is the approximation that provides the explicit relation between and , allowing one to translate also in the other direction, that is, finding if one knows . This is necessary for distilling a closed equation for —as is the Boltzmann equation—out of equation (24). (See again appendix C for a more detailed explanation.)

There are several possible candidates for the required approximation involved in the Ansatz. The considered Ansatzes can for in the lowest Born approximation be written as

 GKBAG<(t1,t2)=iG0R(t1,t2)f(t2)−if(t1)G0A(t1,t2)AAG<(t1,t2)=if(t1)G0R(t1,t2)−iG0A(t1,t2)f(t2)SKBAG<(t1,t2)=12A0(t1,t2)f(t2)+12f(t1)A0(t1,t2) (30)

where is the spectral function. All Ansatzes are consistent with , but go beyond the definition by providing information also about the time-off-diagonal part of . Note that the Ansatz SKBA is the symmetrization of GKBA and AA. In equation (30) the spacial two-point indices have been suppressed for convenience, and so has on the right hand side also the convolution product in these. (In the time indices there is no convolution product.) For the background to these Ansatzes, as well as for the approximations that follow upon the Wigner transformation and with the Born approximation, see ref. Kailasvuori and Lüffe, 2010. Here we just remention that the combination G1AA is what provides the missing link between the Green’s function derived collision integral and the Liouville equation derived one (as derived in appendix B), at least in the lowest Born approximation, which is what we consider. Let us also stress that the list of approaches need not be exhaustive, but provides only a set of natural candidates.

The principal value part of the collision integral can be written as where is the part that does not change sign when one compares the approach G1GKBA to the approaches G1AA, whereas is the part the changes sign. The principal value terms for the other approaches can be composed from these two parts and is given by the table

 XYG1GKBA++\parG1AA+−\parG1SKBA+0\parG2GKBA0+\parG2AA0−\parG2SKBA00. (31)

Using the matrix decomposition in eq. (8) the Boltzmann equation including the collision integral will in a linear approximation in the electric field decompose in the Fourier modes of the angular variables of the momentum . For the terms , with being searched for linear-in-electric-field correction to the known equilibrium part one has according to ref. Kailasvuori and Lüffe, 2010 (for notational convenience the superscript is suppressed; the prime superscript indicates that the component of is a function of and not of )111 In ref. Kailasvuori and Lüffe, 2010 as well as in the three first arXiv versionsof that paper the following formulae contained important mistakes. These mistakes turn out to have no impact on the results presented for the first quantum correction in monolayer graphene discussed paper. However, they have an important impact on the multilayer case studied in the present paper. The following equations are corrected equations, also to appear in corrected versions of the arXiv version of ref. Kailasvuori and Lüffe, 2010 .

 JPX% \boldmath{^b}1=0JPY\boldmath{^b}1=+i∫% \boldmath{k′}W\boldmath{k}\boldmath{k′}2πP−sinNcΔθsinΔθf′z1,JPX\boldmath{^c}1=−∫% \boldmath{k′}W\boldmath{k}\boldmath{k′}2π(−P−cosNcΔθfz1+P+cosΔθf′z1)JPY\boldmath{^c}1=−∫% \boldmath{k′}W\boldmath{k}\boldmath{k′}2π(−P+fz1+P−cosNcΔθcosΔθf′z1)JPXz1=−∫\boldmath{k′}W\boldmath{k}\boldmath{k′}2π{P−cosNcΔθf\boldmath{^c}1−P+(isinNcΔθsinΔθf′% \boldmath{^b}1+cosNcΔθcosΔθf′\boldmath{^c}1)}JPYz1=−∫\boldmath{k′}W\boldmath{k}\boldmath{k′}2π(P+f\boldmath{^c}1−P−cosθΔf′% \boldmath{^c}1) (32)

with

 P±=Pb+b′±Pb−b′. (33)

In the case of monolayer graphene, both and are typically non-zero, and the six candidates in (31) lead to six results—all different—for the leading quantum correction. In contrast, for the case of (including the multilayer Hamiltonians (2) with ) combined with point-like impurities () most of the terms in (32) vanish trivially, namely all the terms that contain trigonometric functions. This leaves as only nonzero terms

 JPY% \boldmath{^c}1=+fz1∫\boldmath{k′}W\boldmath{k}\boldmath{k′}2πP+=:+fz1ΛkJPYz1=−f\boldmath{^c}1∫% \boldmath{k′}W\boldmath{k}\boldmath{k′}2πP+=:−f\boldmath{^c}1Λk.

The vanishing of all terms implies that we cannot distinguish between the approaches and . We henceforth assume that we deal with and suppress this prefix. More interestingly, the fact that some terms survive means that we can distinguish between the different Ansatzes GKBA, AA and SKBA. Furthermore, the integral in can be worked out analytically as (see appendix A)

 Λk = (34)

Note that the derivation requires .

Since the surviving principal value terms only involve but not we can plug out the distribution function out of the integrals, which allows us to turn the integro-differential equations into differential equations, as is typically the case for static and uniform Boltzmann equations. Thanks to this we essentially have the same type of Boltzmann equation as when the principal value terms were neglected. Section 6 in ref. Kailasvuori and Lüffe, 2010 treated the latter problem in a general setting and the results can be easily translated to the present case. Equation (69) in ref. Kailasvuori and Lüffe, 2010 is in the present case replaced by the following equation

 Ex−iEy2⎛⎜ ⎜ ⎜ ⎜ ⎜ ⎜ ⎜⎝∂kfeq0∂kfeq\boldmath{^b}iNckfeq\boldmath{^b}0⎞⎟ ⎟ ⎟ ⎟ ⎟ ⎟ ⎟⎠=−⎛⎜ ⎜ ⎜⎝I0000I0000I2b−νΛ00νΛ−2bI⎞⎟ ⎟ ⎟⎠⎛⎜ ⎜ ⎜ ⎜ ⎜ ⎜ ⎜⎝f(E)01f(E)\boldmath{^b}1f(E)\boldmath{^c}1f(E)z1⎞⎟ ⎟ ⎟ ⎟ ⎟ ⎟ ⎟⎠. (35)

Most importantly one has the simple shift in the strength of the precession term due to the principal value terms, with the sign depending on Ansatz according to

 GKBAν=+1AAν=−1SKBAν=0.

Furthermore, in the right hand side. The chiral origin of that can be found above equation (65) in ref. Kailasvuori and Lüffe, 2010. Other changes are that , and coalesce to , in the present paper denoted , and that in the present setup. The result in equation (75) in ref. Kailasvuori and Lüffe, 2010 becomes in the present case

 σI=σDσII=14πNc2∫∞kFdkbkI(2b−νΛ)2+I2. (36)

The result that there are no quantum corrections to gives an explanation from the perspective of the Boltzmann approach to why one of the central assumptions in section III made sense. There it was assumed that contains all the quantum corrections due to interband effects, whereas is identical to the Drude conductivity , so contains no corrections. The contribution relies on the term . However, according to definition of current in terms of the Boltzmann distribution fucntion in equation (7), the velocity component is what the distribution function component combines with. In the present case this component authors singel-handedly the correction . To see this, look at eq. (72) and eq. (75) in ref. Kailasvuori and Lüffe, 2010 and remember that in our case. Note that for the monolayer problem there is no such clean separation. There also , which contributes to the current through , contains quantum corrections, wherefore we expect that for the monolayer case the crucial assumption in section III not to be valid.

Even in the presence of principal value terms the correction term can be rendered on a neat analytical form. With and reintroduced, we arrive at the main result (see appendix A)

 σI = σD=e2hℓkF2=e2hNdbFτtr2ℏ