The electronic properties of bilayer graphene
We review the electronic properties of bilayer graphene, beginning with a description of the tight-binding 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 tight-binding parameters of the Slonczewski-Weiss-McClure model of bulk graphite plus intra- and interlayer asymmetry between atomic sites which induce band gaps in the low-energy spectrum. The Hartree model of screening and band-gap opening due to interlayer asymmetry in the presence of external gates is presented. The tight-binding 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.
- I Introduction
II Electronic band structure
- II.1 The crystal structure and the Brillouin zone
- II.2 The tight-binding model
- II.3 Effective two-band 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 on-site parameter : electron-hole asymmetry
- II.7 Asymmetry between on-site energies: band gaps
- II.8 Next-nearest neighbour hopping
- II.9 Spin-orbit 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
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 Dirac-like 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 Dirac-like Hamiltonian of monolayer graphene and the second (after the monolayer) in a family of chiral Hamiltonians that appear at low energy in ABC-stacked (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 (); zhang-wang11 ()); 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 (); high-frequency transistors xia10 (); thermoelectric devices wang11 (); photonic devices including plasmonic devices yan12 () and photodetectors yan-kim12 (); 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 low-energy 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 low-energy 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 tight-binding Hamiltonian and resulting band structure describing the low-energy 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 density-dependence of the band gap by taking into account screening by electrons on the bilayer device. The tight-binding 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 electronic-interaction effects. Note that this review considers Bernal-stacked (also known as AB-stacked) bilayer graphene; we do not consider other stacking types such as AA-stacked 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
where is the lattice constant, the distance between adjacent unit cells, Åsaito (). Note that the lattice constant is distinct from the carbon-carbon 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 ‘non-dimer’ 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 (); kosh-mcc10 () 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
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 tight-binding model
ii.2.1 An arbitrary crystal structure
In the following, we will describe the tight-binding 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
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
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
where is a column vector, . The transfer integral matrix and overlap integral matrix are matrices with matrix elements defined by
The band energies may be determined from the generalised eigenvalue equation (5) by solving the secular equation
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 tight-binding model to graphene, and refer the reader to tutorial-style 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
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 off-diagonal 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 nearest-neighbour sites, labelled by index :
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 nearest-neighbour sites is identical for every site. A hopping parameter may be defined as
which is positive. Then, the matrix element may be written as
The other off-diagonal matrix element is given by . The function describing nearest-neighbor hopping, equation (11), is given by
where is the in-plane wave vector. The calculation of the off-diagonal elements of the overlap integral matrix is similar to those of . A parameter is introduced to describe the possibility of non-zero overlap between orbitals on adjacent sites, giving .
Gathering the matrix elements, the transfer and overlap integral matrices of monolayer graphene may be written as
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 non-equivalent (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 Dirac-like 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 tight-binding 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 tight-binding parameters are defined as
Here, we use the notation of the Slonczewski-Weiss-McClure (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.0111This parameter was not determined by the given experiment, the value quoted was taken from previous literature.||-||3.16(3)||3.1111This 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.39111This parameter was not determined by the given experiment, the value quoted was taken from previous literature.|
|0.315(15)||0.10||0.3111This parameter was not determined by the given experiment, the value quoted was taken from previous literature.||-||0.38(6)||0.315111This parameter was not determined by the given experiment, the value quoted was taken from previous literature.|
The upper-left and lower-right square blocks of describe intra-layer terms and are simple generalisations of the monolayer, equation (13). For bilayer graphene, however, we include parameters describing the on-site 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 ():
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 non-dimer 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 upper-right and lower-left square blocks of describe inter-layer 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 in-plane hopping. Parameter describes interlayer coupling between non-dimer orbitals and , and describes interlayer coupling between dimer and non-dimer orbitals and or and . Both and couplings are ‘skew’: they are not strictly vertical, but involve a component of in-plane hopping, and each atom on one layer (e.g., for ) has three equidistant nearest-neighbours (e.g., for ) on the other layer. In fact, the in-plane component of this skew hopping is analogous to nearest-neighbour hopping within a single layer, as parameterised by . Hence, the skew interlayer hopping (e.g., ) contains the factor describing in-plane hopping.
It is possible to introduce an overlap integral matrix for bilayer graphene mucha10 ()
with a form that mirrors . Here, we only include two parameters: describing non-orthogonality of intra-layer nearest-neighbours and describing non-orthogonality 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 non-orthogonality 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 tight-binding 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 Bernal-stacked trilayer graphene by observation of Landau level crossings tay11 (). Note that there are seven parameters in the Slonczewski-Weiss-McClure (SWM) model of graphite sw58 (); mcclure57 (); mcclure60 (); dressel02 () because the next-nearest 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 non-dimer 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 anti-bonding pair arising from the strong coupling (by interlayer coupling ) of the orbitals on the dimer and sites, whereas the ‘low-energy’ bands arise from hopping between the non-dimer and sites. In pristine bilayer graphene, the Fermi level lies at the point where the two low-energy 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 tight-binding 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 tight-binding parameters in this way, generally in line with the experimental ones listed in Table 1. The tight-binding 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 tight-binding 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 non-orthogonality of adjacent orbitals has been neglected here, but it contributes electron-hole asymmetry which is particularly prevalent near the point at the centre of the Brillouin zone saito (); mcc12 ().
ii.2.4 Effective four-band 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
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
where we introduced the effective velocities, and .
At zero magnetic field Hamiltonian (30) yields four valley-degenerate bands . A simple analytic solution may be obtained by neglecting the terms , proportional to , and by considering only interlayer asymmetry in the on-site energies: and . Then, there is electron-hole symmetry, i.e., energies may be written , , mcc06a () with
where is the polar angle of momentum . Energy describes the higher-energy bands split from zero energy by the interlayer coupling between the orbitals on the dimer sites , .
Low-energy bands are related to orbitals on the non-dimer sites , . In an intermediate energy range it is possible to neglect the interlayer asymmetry and terms proportional to (i.e., set ), and the low-energy bands may be approximated mcc06a () as
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 two-component Hamiltonian describing the orbitals on the non-dimer 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 two-band Hamiltonian at low energy
In this section we focus on the low-energy 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 low-energy 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 low-energy and dimer components:
The second row of (33) allows the dimer components to be expressed in terms of the low-energy ones:
Substituting this into the first row of (33) gives an effective eigenvalue equation written solely for the low-energy components:
where . The second equation is accurate up to linear terms in . Finally, we perform a transformation :
This transformation ensures that normalisation of is consistent with that of the original states:
ii.3.2 Bilayer graphene
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
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 electron-hole 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 time-inversion symmetry at zero magnetic field. The intrinsic terms , , , and satisfy spatial-inversion symmetry because the bilayer crystal structure is spatial-inversion symmetric, but terms and , with on the diagonal, are imposed by external fields and they violate spatial-inversion 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 Dirac-like Hamiltonian of monolayer graphene, but with a quadratic-in-momentum term on the off-diagonal rather than linear. For example, the term accounts for an effective hopping between the non-dimer 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
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 Dirac-like 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 ABC-stacked (rhombohedral) multilayer graphene mcc06a (); guinea06 (); koshino07 (); manes07 (); nak08 (); min08a (); min08b (); koshinoABC (); zhang10 ():
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 iso-energetic 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 iso-energetic 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 non-dimer and sites. The influence of on the band structure is described by equation (31). In the two-band 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 iso-energetic 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
in agreement with equation (31) for ,