The electronic properties of bilayer graphene
Abstract
We review the electronic properties of bilayer graphene, beginning with a description of the tightbinding model of bilayer graphene and the derivation of the effective Hamiltonian describing massive chiral quasiparticles in two parabolic bands at low energy. We take into account five tightbinding parameters of the SlonczewskiWeissMcClure model of bulk graphite plus intra and interlayer asymmetry between atomic sites which induce band gaps in the lowenergy spectrum. The Hartree model of screening and bandgap opening due to interlayer asymmetry in the presence of external gates is presented. The tightbinding model is used to describe optical and transport properties including the integer quantum Hall effect, and we also discuss orbital magnetism, phonons and the influence of strain on electronic properties. We conclude with an overview of electronic interaction effects.
Contents
 I Introduction

II Electronic band structure
 II.1 The crystal structure and the Brillouin zone
 II.2 The tightbinding model
 II.3 Effective twoband Hamiltonian at low energy
 II.4 Interlayer coupling : massive chiral electrons
 II.5 Interlayer coupling : trigonal warping and the Lifshitz transition
 II.6 Interlayer coupling and onsite parameter : electronhole asymmetry
 II.7 Asymmetry between onsite energies: band gaps
 II.8 Nextnearest neighbour hopping
 II.9 Spinorbit coupling
 II.10 The integer quantum Hall effect
 III Tuneable band gap
 IV Transport properties
 V Optical properties
 VI Orbital magnetism
 VII Phonons and strain
 VIII Electronic interactions
 IX Summary
I Introduction
The production by mechanical exfoliation of isolated flakes of graphene with excellent conducting properties novo04 () was soon followed by the observation of an unusual sequence of plateaus in the integer quantum Hall effect in monolayer graphene novo05 (); zhang05 (). This confirmed the fact that charge carriers in monolayer graphene are massless chiral quasiparticles with a linear dispersion, as described by a Diraclike effective Hamiltonian d+m84 (); semenoff84 (); haldane88 (), and it prompted an explosion of interest in the field cnreview ().
The integer quantum Hall effect in bilayer graphene novo06 () is arguably even more unusual than in monolayer because it indicates the presence of massive chiral quasiparticles mcc06a () with a parabolic dispersion at low energy. The effective Hamiltonian of bilayer graphene may be viewed as a generalisation of the Diraclike Hamiltonian of monolayer graphene and the second (after the monolayer) in a family of chiral Hamiltonians that appear at low energy in ABCstacked (rhombohedral) multilayer graphene mcc06a (); guinea06 (); koshino07 (); manes07 (); nak08 (); min08a (); min08b (). In addition to interesting underlying physics, bilayer graphene holds potential for electronics applications, not least because of the possibility to control both carrier density and energy band gap through doping or gating mcc06a (); guinea06 (); ohta06 (); mcc06b (); min07 (); oostinga (); castro ().
Not surprisingly, many of the properties of bilayer graphene are similar to those in monolayer cnreview (); novo12 (). These include excellent electrical conductivity with room temperature mobility of up to cmVs in air dean10 (); the possibility to tune electrical properties by changing the carrier density through gating or doping novo04 (); novo06 (); ohta06 (); high thermal conductivity with room temperature thermal conductivity of about WmK ghosh10 (); balandin11 (); mechanical stiffness, strength and flexibility (Young’s modulus is estimated to be about TPa neekamal10 (); zhangwang11 ()); transparency with transmittance of white light of about nair08 (); impermeability to gases bunch08 (); and the ability to be chemically functionalised elias09 (). Thus, as with monolayer graphene, bilayer graphene has potential for future applications in many areas novo12 () including transparent, flexible electrodes for touch screen displays bae10 (); highfrequency transistors xia10 (); thermoelectric devices wang11 (); photonic devices including plasmonic devices yan12 () and photodetectors yankim12 (); energy applications including batteries sugawara11 (); kanetani12 (); and composite materials gong12 (); young12 ().
It should be stressed, however, that bilayer graphene has features that make it distinct from monolayer. The lowenergy band structure, described in detail in Section II, is different. Like monolayer, intrinsic bilayer has no band gap between its conduction and valence bands, but the lowenergy dispersion is quadratic (rather than linear as in monolayer) with massive chiral quasiparticles novo06 (); mcc06a () rather than massless ones. As there are two layers, bilayer graphene represents the thinnest possible limit of an intercalated material sugawara11 (); kanetani12 (). It is possible to address each layer separately leading to entirely new functionalities in bilayer graphene including the possibility to control an energy band gap of up to about meV through doping or gating mcc06a (); guinea06 (); ohta06 (); mcc06b (); min07 (); oostinga (); castro (). Recently, this band gap has been used to create devices  constrictions and dots  by electrostatic confinement with gates goossens12 (). Bilayer or multilayer graphene devices may also be preferable to monolayer ones when there is a need to use more material for increased electrical or thermal conduction, strength gong12 (); young12 (), or optical signature yan12 ().
In the following we review the electronic properties of bilayer graphene. Section II is an overview of the electronic tightbinding Hamiltonian and resulting band structure describing the lowenergy chiral Hamiltonian and taking into account different parameters that couple atomic orbitals as well as external factors that may change the electron bands by, for example, opening a band gap. We include the Landau level spectrum in the presence of a perpendicular magnetic field and the corresponding integer quantum Hall effect. In section III we consider the opening of a band gap due to doping or gating and present a simple analytical model that describes the densitydependence of the band gap by taking into account screening by electrons on the bilayer device. The tightbinding model is used to describe transport properties, section IV, and optical properties, section V. We also discuss orbital magnetism in section VI, phonons and the influence of strain in section VII. Section VIII concludes with an overview of electronicinteraction effects. Note that this review considers Bernalstacked (also known as ABstacked) bilayer graphene; we do not consider other stacking types such as AAstacked graphene liu09 (), twisted graphene lopes07 (); berger06 (); hass08 (); mele10 (); li10 (); luican11 () or two graphene sheets separated by a dielectric with, possibly, electronic interactions between them schmidt08 (); ni08 (); min08c (); khar08 (); schmidt10 (); pon11 ().
Ii Electronic band structure
ii.1 The crystal structure and the Brillouin zone
Bilayer graphene consists of two coupled monolayers of carbon atoms, each with a honeycomb crystal structure. Figure 1 shows the crystal structure of monolayer graphene, figure 2 shows bilayer graphene. In both cases, primitive lattice vectors and may be defined as
(1) 
where is the lattice constant, the distance between adjacent unit cells, Åsaito (). Note that the lattice constant is distinct from the carboncarbon bond length Å, which is the distance between adjacent carbon atoms.
In monolayer graphene, each unit cell contains two carbon atoms, labelled and , figure 1(a). The positions of and atoms are not equivalent because it is not possible to connect them with a lattice vector of the form , where and are integers. Bilayer graphene consists of two coupled monolayers, with four atoms in the unit cell, labelled , on the lower layer and , on the upper layer. The layers are arranged so that one of the atoms from the lower layer is directly below an atom, , from the upper layer. We refer to these two atomic sites as ‘dimer’ sites because the electronic orbitals on them are coupled together by a relatively strong interlayer coupling. The other two atoms, and , don’t have a counterpart on the other layer that is directly above or below them, and are referred to as ‘nondimer’ sites. Note that some authors guinea06 (); nil08 (); li09 (); zhang08 () employ different definitions of and sites as used here. The point group of the bilayer crystal structure is latil06 (); manes07 (); koshmcc10 () consisting of elements (), and it may be regarded as a direct product of group () with the inversion group (). Thus, the lattice is symmetric with respect to spatial inversion symmetry .
Primitive reciprocal lattice vectors and of monolayer and bilayer graphene, where and , are given by
(2) 
As shown in figure 1(b), the reciprocal lattice is an hexagonal Bravais lattice, and the first Brillouin zone is an hexagon.
ii.2 The tightbinding model
ii.2.1 An arbitrary crystal structure
In the following, we will describe the tightbinding model ash+mer (); saito (); mcc12 () and its application to bilayer graphene. We begin by considering an arbitrary crystal with translational invariance and atomic orbitals per unit cell, labelled by index . Bloch states for a given position vector and wave vector may be written as
(3) 
where is the number of unit cells, labels the unit cell, and is the position vector of the th orbital in the th unit cell.
The electronic wave function may be expressed as a linear superposition of Bloch states
(4) 
where are expansion coefficients. There are different energy bands, and the energy of the th band is given by where is the Hamiltonian. Minimising the energy with respect to the expansion coefficients saito (); mcc12 () leads to
(5) 
where is a column vector, . The transfer integral matrix and overlap integral matrix are matrices with matrix elements defined by
(6) 
The band energies may be determined from the generalised eigenvalue equation (5) by solving the secular equation
(7) 
where ‘’ stands for the determinant of the matrix.
In order to model a given system in terms of the generalised eigenvalue problem (5), it is necessary to determine the matrices and . We will proceed by considering the relatively simple case of monolayer graphene, before generalising the approach to bilayers. In the following sections, we will omit the subscript on and in equation (5), remembering that the number of solutions is , the number of orbitals per unit cell.
ii.2.2 Monolayer graphene
Here, we will outline how to apply the tightbinding model to graphene, and refer the reader to tutorialstyle reviews saito (); mcc12 () for further details. We take into account one orbital per atomic site and, as there are two atoms in the unit cell of monolayer graphene, figure 1(a), we include two orbitals per unit cell labelled as and (the atoms and the atoms are each arranged on an hexagonal Bravais lattice).
We begin by considering the diagonal element of the transfer integral matrix , equation (6), for the site orbital. It may be determined by substituting the Bloch function (3) for into the matrix element (6), which results in a double sum over the positions of the unit cells in the crystal. Assuming that the dominant contribution arises from those terms involving a given orbital interacting with itself (i.e., in the same unit cell), the matrix element may be written as
(8) 
This may be regarded as a summation over all unit cells of a parameter that takes the same value in every unit cell. Thus, the matrix element may be simply expressed as . Similarly, the diagonal element for the site orbital can be written as , while for intrinsic graphene is equal to as the two sublattices are identical. The calculation of the diagonal elements of the overlap integral matrix , equation (6), proceeds in the same way as that of , with the overlap of an orbital with itself equal to unity, . Thus, .
The offdiagonal element of the transfer integral matrix describes the possibility of hopping between orbitals on and sites. Substituting the Bloch function (3) into the matrix element (6) results in a sum over all sites and a sum over all sites. We assume that the dominant contribution arises from hopping between adjacent sites. If we consider a given site, say, then we take into account the possibility of hopping to its three nearestneighbour sites, labelled by index :
(9)  
where are the positions of three nearest atoms relative to a given atom, which may be written as , , .
The sum with respect to the three nearestneighbour sites is identical for every site. A hopping parameter may be defined as
(10) 
which is positive. Then, the matrix element may be written as
(11) 
The other offdiagonal matrix element is given by . The function describing nearestneighbor hopping, equation (11), is given by
(12) 
where is the inplane wave vector. The calculation of the offdiagonal elements of the overlap integral matrix is similar to those of . A parameter is introduced to describe the possibility of nonzero overlap between orbitals on adjacent sites, giving .
Gathering the matrix elements, the transfer and overlap integral matrices of monolayer graphene may be written as
(13)  
(14) 
The corresponding energy may be determined saito () by solving the secular equation (7). For intrinsic graphene, i.e., , we have
(15) 
The parameter values are listed by Saito et al saito () as eV and .
The function , equation (12) is zero at the corners of the Brillouin zone, two of which are nonequivalent (i.e., they are not connected by a reciprocal lattice vector). For example, corners and with wave vectors are labelled in Figure 1(b). Such positions are called points or valleys, and we will use a valley index to distinguish points . At these positions, the solutions (15) are degenerate, marking a crossing point and zero band gap between the conduction and valence bands. The transfer matrix is approximately equal to a Diraclike Hamiltonian in the vicinity of the point, describing massless chiral quasiparticles with a linear dispersion relation. These points are particularly important because the Fermi level is located near them in pristine graphene.
ii.2.3 Bilayer graphene
In the tightbinding description of bilayer graphene, we take into account orbitals on the four atomic sites in the unit cell, labelled as . Then, the transfer integral matrix of bilayer graphene mcc06a (); guinea06 (); part06 (); mucha08 (); nil08 (); mucha10 () is a matrix given by
where the tightbinding parameters are defined as
(17)  
(18)  
(19)  
(20) 
Here, we use the notation of the SlonczewskiWeissMcClure (SWM) model sw58 (); mcclure57 (); mcclure60 (); dressel02 () that describes bulk graphite. Note that definitions of the parameters used by authors can differ, particularly with respect to signs.
Parameter  Graphite dressel02 ()  Bilayer mal07 ()  Bilayer zhang08 ()  Bilayer li09 ()  Bilayer kuz09b ()  Trilayer tay11 () 

3.16(5)  2.9  3.0^{1}^{1}1This parameter was not determined by the given experiment, the value quoted was taken from previous literature.    3.16(3)  3.1^{1}^{1}1This parameter was not determined by the given experiment, the value quoted was taken from previous literature.  
0.39(1)  0.30  0.40(1)  0.404(10)  0.381(3)  0.39^{1}^{1}1This parameter was not determined by the given experiment, the value quoted was taken from previous literature.  
0.020(2)          0.028(4)  
0.315(15)  0.10  0.3^{1}^{1}1This parameter was not determined by the given experiment, the value quoted was taken from previous literature.    0.38(6)  0.315^{1}^{1}1This parameter was not determined by the given experiment, the value quoted was taken from previous literature.  
0.044(24)  0.12  0.15(4)    0.14(3)  0.041(10)  
0.038(5)          0.05(2)  
0.008(2)    0.018(3)  0.018(2)  0.022(3)  0.03(2)  
0.050(6)    0.018(3)  0.018(2)  0.022(3)  0.046(10)  
The upperleft and lowerright square blocks of describe intralayer terms and are simple generalisations of the monolayer, equation (13). For bilayer graphene, however, we include parameters describing the onsite energies , , , on the four atomic sites, that are not equal in the most general case. As there are four sites, differences between them are described by three parameters mucha10 ():
(21)  
(22)  
(23)  
(24) 
where
(25)  
(26)  
(27) 
The three independent parameters are to describe interlayer asymmetry between the two layers mcc06a (); ohta06 (); guinea06 (); mcc06b (); min07 (); castro (); oostinga (); aoki07 (); mcc07 (); guinea07 (); mcc07b (); guinea07b (); gava (); bouk08 (), for an energy difference between dimer and nondimer sites dressel02 (); nil08 (); zhang08 (); li09 (), and for an energy difference between and sites on each layer mucha09d (); mucha10 (). These parameters are described in detail in sections II.6 and II.7.
The upperright and lowerleft square blocks of describe interlayer coupling. Parameter describes coupling between pairs of orbitals on the dimer sites and : since this is a vertical coupling, the corresponding terms in (i.e., ) do not contain which describes inplane hopping. Parameter describes interlayer coupling between nondimer orbitals and , and describes interlayer coupling between dimer and nondimer orbitals and or and . Both and couplings are ‘skew’: they are not strictly vertical, but involve a component of inplane hopping, and each atom on one layer (e.g., for ) has three equidistant nearestneighbours (e.g., for ) on the other layer. In fact, the inplane component of this skew hopping is analogous to nearestneighbour hopping within a single layer, as parameterised by . Hence, the skew interlayer hopping (e.g., ) contains the factor describing inplane hopping.
It is possible to introduce an overlap integral matrix for bilayer graphene mucha10 ()
(28) 
with a form that mirrors . Here, we only include two parameters: describing nonorthogonality of intralayer nearestneighbours and describing nonorthogonality of orbitals on dimer sites and . In principle, it is possible to introduce additional parameters analogous to , , etc., but generally they will be small and irrelevant. In fact, it is common practice to neglect the overlap integral matrix entirely, i.e., replace with a unit matrix, because the influence of parameters and describing nonorthogonality of adjacent orbitals is small at low energy . Then, the generalised eigenvalue equation (5) reduces to an eigenvalue equation with Hamiltonian , equation (LABEL:Hbfull).
The energy differences and are usually attributed to extrinsic factors such as gates, substrates or doping. Thus, there are five independent parameters in the Hamiltonian (LABEL:Hbfull) of intrinsic bilayer graphene, namely , , , and . The band structure predicted by the tightbinding model has been compared to observations from photoemission ohta06 (), Raman mal07 () and infrared spectroscopy hen08 (); zhang08 (); li09 (); kuz09a (); kuz09b (); mak10 (). Parameter values determined by fitting to experiments are listed in Table 1 for bulk graphite dressel02 (), for bilayer graphene by Raman mal07 (); das () and infrared zhang08 (); li09 (); kuz09b () spectroscopy, and for Bernalstacked trilayer graphene by observation of Landau level crossings tay11 (). Note that there are seven parameters in the SlonczewskiWeissMcClure (SWM) model of graphite sw58 (); mcclure57 (); mcclure60 (); dressel02 () because the nextnearest layer couplings and , absent in bilayer, are present in graphite (and trilayer graphene, too). Parameter in the SWM model is related by to the parameter describing the energy difference between dimer and nondimer sites in bilayer graphene.
The energy bands are plotted in figure 3 along the axis in reciprocal space intersecting the corners , and the centre of the Brillouin zone [see figure 1(b)]. Plots were made using Hamiltonian , equation (LABEL:Hbfull), with parameter values determined by infrared spectroscopy eV, eV, eV, eV, eV, and kuz09b (). There are four bands because the model takes into account one orbital on each of the four atomic sites in the unit cell; a pair of conduction bands and a pair of valence bands. Over most of the Brillouin zone, each pair is split by an energy of the order of the interlayer spacing eV trickey (). Near the points, inset of figure 3, one conduction band and one valence band are split away from zero energy by an energy of the order of the interlayer coupling , whereas two bands touch at zero energy mcc06a (). The ‘split’ bands are a bonding and antibonding pair arising from the strong coupling (by interlayer coupling ) of the orbitals on the dimer and sites, whereas the ‘lowenergy’ bands arise from hopping between the nondimer and sites. In pristine bilayer graphene, the Fermi level lies at the point where the two lowenergy bands touch (shown as zero energy in figure 3) and, thus, this region is relevant for the study of electronic properties. It will be the focus of the following sections.
At low energy, the shape of the band structure predicted by the tightbinding model (see inset in figure 3) is in good agreement with that calculated by density functional theory latil06 (); min07 (); aoki07 () and it is possible obtain values for the tightbinding parameters in this way, generally in line with the experimental ones listed in Table 1. The tightbinding model Hamiltonian , equation (LABEL:Hbfull), used in conjuction with the parameters listed in Table 1, is not accurate over the whole Brillouin zone because the fitting of tightbinding parameters is generally done in the vicinity of the corners of the Brillouin zone and (as the Fermi level lies near zero energy). For example, parameter in equation (28) describing nonorthogonality of adjacent orbitals has been neglected here, but it contributes electronhole asymmetry which is particularly prevalent near the point at the centre of the Brillouin zone saito (); mcc12 ().
ii.2.4 Effective fourband Hamiltonian near the points
To describe the properties of electrons in the vicinity of the points, a momentum is introduced which is measured from the centre of the point. Expanding in powers of , the function , equation (12), is approximately given by which is valid close to the point, i.e., for , where . In monolayer graphene, the Hamiltonian matrix (13) is simplified by keeping only linear terms in momentum as
(29) 
where , , and is the band velocity. In the intrinsic case, , the eigen energy becomes , which approximates Eq. (15). In bilayer graphene, similarly, Eq. (LABEL:Hbfull) is reduced to
(30) 
where we introduced the effective velocities, and .
At zero magnetic field Hamiltonian (30) yields four valleydegenerate bands . A simple analytic solution may be obtained by neglecting the terms , proportional to , and by considering only interlayer asymmetry in the onsite energies: and . Then, there is electronhole symmetry, i.e., energies may be written , , mcc06a () with
(31)  
where is the polar angle of momentum . Energy describes the higherenergy bands split from zero energy by the interlayer coupling between the orbitals on the dimer sites , .
Lowenergy bands are related to orbitals on the nondimer sites , . In an intermediate energy range it is possible to neglect the interlayer asymmetry and terms proportional to (i.e., set ), and the lowenergy bands may be approximated mcc06a () as
(32) 
which interpolates between an approximately linear dispersion at large momentum to a quadratic one at small momentum, where the mass is (see inset in figure 3). This crossover occurs at . A convenient way to describe the bilayer at low energy and momentum is to eliminate the components in the Hamiltonian (30) related to the orbitals on dimer sites , , resulting in an effective twocomponent Hamiltonian describing the orbitals on the nondimer sites , , and, thus, the two bands that approach each other at zero energy. This is described in the next section, and the solutions of this Hamiltonian are shown to be massive chiral quasiparticles novo06 (); mcc06a (), as opposed to massless chiral ones in monolayer graphene.
ii.3 Effective twoband Hamiltonian at low energy
In this section we focus on the lowenergy electronic band structure in the vicinity of the points and at the corners of the first Brillouin zone, relevant for energies near the Fermi level. A simple model may be obtained by eliminating orbitals related to the dimer sites, resulting in an effective Hamiltonian for the lowenergy orbitals. First, we outline the procedure in general terms, because it may be applied to systems other than bilayer graphene such as stacked (rhombohedral) graphene multilayers koshinoABC (); zhang10 (), before applying it specifically to bilayer graphene.
ii.3.1 General procedure
We consider the energy eigenvalue equation, and consider separate blocks in the Hamiltonian corresponding to lowenergy and dimer components:
(33) 
The second row of (33) allows the dimer components to be expressed in terms of the lowenergy ones:
(34) 
Substituting this into the first row of (33) gives an effective eigenvalue equation written solely for the lowenergy components:
where . The second equation is accurate up to linear terms in . Finally, we perform a transformation :
(35) 
This transformation ensures that normalisation of is consistent with that of the original states:
where we used equation (34) for small : . Thus, the effective Hamiltonian for lowenergy components is given by equation (35):
(36)  
(37) 
ii.3.2 Bilayer graphene
The Hamiltonian (30) is written in basis . If, instead, it is written in the basis of lowenergy and dimer components , equation (33), then
Using the procedure described in the previous section, equations (36,37), it is possible to obtain an effective Hamiltonian for components . An expansion is performed by assuming that the intralayer hopping and the interlayer coupling are larger than other energies: . Then, keeping only terms that are linear in the small parameters and quadratic in momentum, the effective Hamiltonian mcc06a (); mucha10 () is
(38)  
where , . In the following sections, we discuss the terms in . The first term describes massive chiral electrons, section II.4. It generally dominates at low energy , so that the other terms in may be considered as perturbations of it. The second term , section II.5, introduces a triangular distortion of the Fermi circle around each point known as ‘trigonal warping’. Terms and , with on the diagonal, produce a band gap between the conduction and valence bands, section II.7, whereas and introduce electronhole asymmetry into the band structure, section II.6.
The Hamiltonian (38) is written in the vicinity of a valley with index distinguishing between and . In order to briefly discuss the effect of symmetry operations on it, we introduce Pauli spin matrices , , in the / sublattice space and , , in the valley space. Then, the first term in the Hamiltonian may be written as . The operation of spatial inversion is represented by because it swaps both valleys and lattice sites, time inversion is given by complex conjugation and , as it swaps valleys, too. Hamiltonian (38) satisfies timeinversion symmetry at zero magnetic field. The intrinsic terms , , , and satisfy spatialinversion symmetry because the bilayer crystal structure is spatialinversion symmetric, but terms and , with on the diagonal, are imposed by external fields and they violate spatialinversion symmetry, producing a band gap between the conduction and valence bands.
ii.4 Interlayer coupling : massive chiral electrons
The Hamiltonian in equation (38) resembles the Diraclike Hamiltonian of monolayer graphene, but with a quadraticinmomentum term on the offdiagonal rather than linear. For example, the term accounts for an effective hopping between the nondimer sites , via the dimer sites , consisting of a hop from to (contributing a factor ), followed by a transition between , dimer sites (giving a ‘mass’ ), and a hop from to (a second factor of ). The solutions are massive chiral electrons novo06 (); mcc06a (), with parabolic dispersion , . The density of states is per spin and per valley, and the Fermi velocity is momentum dependent, unlike the Fermi velocity of monolayer graphene.
The corresponding wave function is given by
(39) 
The wave function components describe the electronic amplitudes on the and sites, and it can be useful to introduce the concept of a pseudospin degree of freedom novo06 (); mcc06a () that is related to these amplitudes. If all the electronic density were located on the sites, then the pseudospin part of the wave function could be viewed as a pseudospin ‘up’ state, pointing upwards out of the graphene plane. Likewise, a state with density solely on the sites could be viewed as a pseudospin ‘down’ state. However, density is usually shared equally between the two layers, so that the pseudospin is a linear combination of up and down, , equation (39), and it lies in the graphene plane.
The Hamiltonian may also be written as where the pseudospin vector is , and is a unit vector. This illustrates the chiral nature of the electrons novo06 (); mcc06a (): the chiral operator projects the pseudospin onto the direction of quantization , which is fixed to lie in the graphene plane, but turns twice as quickly as the momentum . For these chiral quasiparticles, adiabatic propagation along a closed orbit produces a Berry’s phase berry84 () change of novo06 (); mcc06a () of the wave function, in contrast to Berry phase in monolayer graphene.
Note that the chiral Hamiltonian may be viewed as a generalisation of the Diraclike Hamiltonian of monolayer graphene and the second (after the monolayer) in a family of chiral Hamiltonians , , corresponding to Berry’s phase which appear at low energy in ABCstacked (rhombohedral) multilayer graphene mcc06a (); guinea06 (); koshino07 (); manes07 (); nak08 (); min08a (); min08b (); koshinoABC (); zhang10 ():
(40) 
where for monolayer, for bilayer, and for trilayer graphene. Since the pseudospin is related to the wavefunction amplitude on sites that are located on different layers, pseudospin may be viewed as a ‘which layer’ degree of freedom min08a (); zhang11 ().
ii.5 Interlayer coupling : trigonal warping and the Lifshitz transition
The Hamiltonian in equation (38) yields a quadratic, isotropic dispersion relation with circular isoenergetic lines, i.e., there is a circular Fermi line around each point. This is valid near the point, , whereas, at high energy, and momentum far from the point, there is a triangular perturbation of the circular isoenergetic lines known as trigonal warping, as in monolayer graphene and graphite. It occurs because the band structure follows the symmetry of the crystal lattice as described by the full momentum dependence of the function , equation (12) ando98 ().
In bilayer graphene mcc06a (), as in bulk graphite dre74 (); nak76 (); ino62 (); gup72 (), a second source of trigonal warping arises from the skew interlayer coupling between nondimer and sites. The influence of on the band structure is described by equation (31). In the twoband Hamiltonian, it is described by in equation (38), the second term of which arises from a quadratic term in the expansion of . This second term has the same momentum dependence as , and, thus, it actually only gives a small additional contribution to the mass . The first term in causes trigonal warping of the isoenergetic lines in directions , where at , at .
To analyse the influence of at low energy, we consider just and the first term in , and the resulting energy is given by
(41) 
in agreement with equation (31) for ,