Electronic properties of bilayer and multilayer graphene

Electronic properties of bilayer and multilayer graphene

Johan Nilsson Department of Physics, Boston University, 590 Commonwealth Avenue, Boston, MA 02215, USA Instituut-Lorentz, Universiteit Leiden, P.O. Box 9506, 2300 RA Leiden, The Netherlands    A. H. Castro Neto Department of Physics, Boston University, 590 Commonwealth Avenue, Boston, MA 02215, USA    F. Guinea Instituto de Ciencia de Materiales de Madrid, CSIC, Cantoblanco E28049 Madrid, Spain    N. M. R. Peres Center of Physics and Departamento de Física, Universidade do Minho, P-4710-057, Braga, Portugal
January 18, 2007
Abstract

We study the effects of site dilution disorder on the electronic properties in graphene multilayers, in particular the bilayer and the infinite stack. The simplicity of the model allows for an easy implementation of the coherent potential approximation and some analytical results. Within the model we compute the self-energies, the density of states and the spectral functions. Moreover, we obtain the frequency and temperature dependence of the conductivity as well as the DC conductivity. The c-axis response is unconventional in the sense that impurities increase the response for low enough doping. We also study the problem of impurities in the biased graphene bilayer.

pacs:
81.05.Uw 73.21.Ac 71.23.-k

I Introduction

The isolation of single layer graphene by Novoselov et al. Novoselov et al. (2004) has generated enormous interest in the physics community. On the one hand, the electronic excitations of graphene can be described by the two-dimensional (2D) Dirac equation, creating connections with certain theories in particle physics.Castro Neto et al. (2006) Moreover, the “relativistic” nature of the quasiparticles, albeit with a speed of propagation, , 300 times smaller than the speed of light, leads to unusual spectroscopic, transport, and thermodynamic properties that are at odds with the standard Landau-Fermi liquid theory of metals.Castro Neto et al. (2007) On the other hand, graphene opens the doors for an all-carbon based micro-electronics.Geim and Novoselov (2007)

Due to the strong nature of the bonds in graphene, and strong mechanical stability of the graphene lattice, miniaturization can be obtained at sizes of order of a few nanometers, beyond what can obtained with the current silicon technology (the smallest size being of the order of the benzene molecule). Furthermore, the same stability allows for creation of entire devices (transistors, wires, and contacts) carved out of the same graphene sheet, reducing tremendously the energy loss, and hence heating, created by contacts between different materials.Berger et al. (2004) Early proposals for the control of the electronic properties in graphene, such as the opening of gaps, were based on controlling its geometry, either by reducing it to nanoribbons,Nakada et al. (1996) or producing graphene quantum dots.Silvestrov and Efetov (2007) Nevertheless, current lithographic techniques that can produce such nanostructures do not have enough accuracy to cut graphene to Ångstrom precision. As a result, graphene nanostructures unavoidably have rough edges which have strong effects in the transport properties of nanoribbons.Han et al. (2007) In addition, the small size of these structures can lead to strong enhancement of the Coulomb interaction between electrons which, allied to the disorder at the edge of the nanostructures, can lead to Coulomb blockade effects easily observable in transport and spectroscopy.Sols et al. (2007)

Hence, the control of electronic gaps by finite geometry is still very unreliable at this point in time and one should look for control in bulk systems which are insensitive to edge disorder. Fortunately, graphene is an extremely flexible material from the electronic point of view and electronic gaps can be controlled. This can be accomplished in a graphene bilayer with an electric field applied perpendicular to the plane. It was shown theoretically McCann and Fal’ko (2006); McCann (2006) and demonstrated experimentally Castro et al. (2007a); Oostinga et al. (2007) that a graphene bilayer is the only material with semiconducting properties that can be controlled by electric field effect. The size of the gap between conduction and valence bands is proportional to the voltage drop between the two graphene planes and can be as large as eV, allowing for novel terahertz devicesCastro et al. (2007a) and carbon-based quantum dotsJ. Milton Pereira et al. (2007) and transistors.Nilsson et al. (2007)

Nevertheless, just as single layer graphene,Peres et al. (2006) bilayer graphene is also sensitive to the unavoidable disorder generated by the environment of the SiO substrate: adatoms, ionized impurities, etc. Disorder generates a scattering rate and hence a characteristic energy scale which is the order of the Fermi energy ( is the Fermi momentum and is the planar density of electrons) when the chemical potential is close to the Dirac point (). Thus, one expects disorder to have a strong effect in the physical properties of graphene. Indeed, theoretical studies of the effect of disorder in unbiased Nilsson et al. (2006a) and biased Nilsson and Castro Neto (2007) graphene bilayer (and multilayer) show that disorder leads to strong modifications of its transport and spectroscopic properties. The understanding of the effects of disorder in this new class of materials is fundamental for any future technological applications. In this context it is worth to mention the transport theories based on the Boltzmann equation,Katsnelson (2007); Adam and Sarma () a study of weak localization in bilayer graphene,Kechedzhi et al. (2007) and also corresponding further experimental characterization.Morozov et al. (); Gorbachev et al. (2007) DC transport in few-layer graphene systems have been studied in Ref. Nakamura and Hirasawa, , both without and in the presence of a magnetic field.

In this paper, we study the effects of site dilution (or unitary scattering) on the electronic properties of graphene multilayers within the well-known coherent potential approximation (CPA). While the CPA does not take into account electron localization,not (); Ziegler (2006) it does provide quantitative and qualitative information on the effect of disorder in the electronic excitations. Furthermore, this approximation allows for analytical results of electronic self-energies, allowing us to compute physical quantities such as spectral functions (measurable by angle resolved photoemission, ARPES Zhou et al. (2006a, b); Bostwick et al. (2007); Ohta et al. (2006); Zhou et al. (2007)) and density of states (measurable by scanning tunneling microscopy, STMStolyarova et al. (2007); Mallet et al. (2007); Li and Andrei (2007); Brar et al. (2007)), besides standard transport properties such as the DC and AC conductivities.Nilsson et al. (2006a) Furthermore, in the case of the semi-infinite stack of graphene planes we can compute the c-axis response of the system which is rather unusual since it increases with disorder at low electronic densities, in agreement with early transport measurements in graphite.Brandt et al. (1988)

The paper is organized as follows. In Sec. II we discuss the band model of the graphene bilayer within the tight-binding approximation. We also connect our notation with the one established for graphite, namely the Slonsczewki-Weiss-McClure (SWM) parameterization. In Sec. III we introduce several simplified band models and compare the electronic bands in different approximations. The Green’s functions that we use later on in the paper are given in Sec. IV.

We employ a simplified model for the disordered graphene bilayer in Sec. V and work out the consequences on the single particle properties encoded in the self-energies, the density of states (DOS) and the spectral function. Sec. VI contains results for the graphene multilayer. In Sec. VII we introduce the linear response formulas that we use to calculate the electronic and optical response. The results for the conductivities in the bilayer are presented in Sec. VIII, while those for the multilayer can be found in Sec. IX.

The rest of the paper concerns the problem of impurities in the biased graphene bilayer. The model of the system and some of its basic properties are discussed in Sec. X. In Sec. XI we solve the problem of a Dirac delta impurity exactly within the effective mass approximation. A simple estimate of when the interactions among impurities becomes important is presented in Sec. XII. We treat more general impurity potentials with variational methods in Sec. XIII, and the special case of a potential well with finite range is studied in Sec. XIV. In Sec. XV we study the problem of a finite density of impurities in the coherent potential approximation (CPA). The effects of trigonal distortions on our results for the biased graphene bilayer are discussed briefly in Sec. XVI. Finally, the conclusions of the paper are to be found in Sec. XVII. We have also included four appendices with technical details of the calculations of the minimal conductivity in bilayer graphene (App. A), the DOS in multilayer graphene (App. B), the conductivity kernels (App. C), and the Green’s function in the biased graphene bilayer (App. D).

Ii Electronic bands of the graphene bilayer

Many of the special properties of the graphene bilayer have their origin in its lattice structure that leads to the peculiar band structure that we discuss in detail in this section. A simple way of arriving at the band structure of the graphene bilayer is to use a tight-binding approximation. The positions of the different atoms in the graphene bilayer are shown in Fig. 1 together with our labeling convention.

The advantage of this notation is that one can discuss collectively about the A (B) atoms that are equivalent in their physical properties such as the weight of the wave functions and the distribution of the density of states etc. This notation was used in early work on graphite.McClure (1957); Slonczewski and Weiss (1958) Many authors use instead a notation similar to , , , and . In this notation the relative orientation within the planes of the A () and B () atoms are the same; but for the other physical properties the equivalent atoms are instead A (B) and (). Because the other physical properties are often more relevant for the physics than the relative orientation of the atoms within the planes we choose to use the, in our view, most “natural” labeling convention.

ii.1 Monolayer graphene

Let us briefly review the tight-binding model of monolayer graphene.Wallace (1947) The band structure can be described in terms of a triangular lattice with two atoms per unit cell. The real-space lattice vectors can be taken to be and . Here () denotes the nearest neighbor carbon distance. Three vectors that connect atoms that are nearest neighbors are , , and ; we take these to connect the atoms to the atoms. In terms of the operators that creates (annihilates) an electron on the lattice site at position and lattice site [ denotes the atom sublattice and () denotes the plane ]: (), the tight-binding Hamiltonian reads:

(1)

Here () is the energy associated with the hopping of electrons between neighboring orbitals. We define the Fourier-transformed operators,

(2)

where is the number of unit cells in the system. Throughout this paper we use units such that unless specified otherwise.

Because of the sublattice structure it is often convenient to describe the system in terms of a spinor: , in which case the Hamiltonian can be written as:

(3)

where

(4)

The reciprocal lattice vectors can be taken to be and as is readily verified. The center of the Brillouin zone (BZ) is denoted by , but for the low-energy properties one can expand close to the K point of the BZ, which has coordinates . One then finds , where and . Note that along the line of the BZ and that it increases anti-clockwise. With these approximation one finds that the spectrum of Eq. (3) is that of massless 2D Dirac fermions: .

Figure 1: [color online] Lattice structure of the graphene bilayer. The A (B) sublattices are indicated by the darker (lighter) spheres and the planes are labeled by 1 and 2.

ii.2 Bilayer graphene

Since the system is 2D only the relative position of the atoms projected on to the --plane enters into the model. The projected position of the different atoms are shown in Fig. 2 together with the BZ.

Figure 2: [color online] The real space lattice structure of the graphene bilayer projected onto the - plane showing the relative positions of the different sublattices. The upper right corner shows the BZ of the graphene bilayer including the labeling of the high symmetry points.

Since the A atoms are sitting right on top of each other in the lattice, the hopping term between the and atoms are local in real space and hence a constant that we denote by in momentum space. Referring back to Section II.1 we note that the hopping [] gives rise to the factor [], with defined in Eq. (4). Since the geometrical role of the A and B atoms are interchanged between plane 1 and plane 2 we immediately find that in Fourier space the hopping [] gives rise to the factor []. Furthermore, the direction in the hopping (projected on to the - plane) is opposite to that of hopping . Thus we associate a factor to the hopping , where the factor is needed because the hopping energy is instead of . Similarly, the direction of hopping (projected on to the - plane) is the same as and therefore the term goes with the hopping . The minus sign in front of follows from the conventional definition of in graphite, as are discussed below. Continuing to fill in all the entries of the matrix the full tight-binding Hamiltonian in the graphene bilayer becomes:

(5)

where the spinor is . Here we have also introduced the conventional (from graphite) that parametrizes the difference in energy between A and B atoms. In addition we included the parameter which gives different values of the potential energy in the two planes, such a term is generally allowed by symmetry and is generated by an electric field that is perpendicular to the two layers. The system with is called the biased graphene bilayer and has a gap in the spectrum, in contrast the spectrum is gapless if .

It is also possible to include further hoppings into the tight-binding picture, this was done for graphite by Johnson and Dresselhaus.Johnson and Dresselhaus (1973) The inclusion of such terms is necessary if one wants an accurate description of the bands throughout the whole BZ. If we expand the expression in Eq. (5) close to the K point in the BZ we obtain the matrix:

(6)

where was introduced after Eq. (4).

The typical behavior of the bands obtained from Eq. (6) is shown in Fig. 3. Two of the bands are moved away from the Dirac point by an energy that is approximately given by the interplane hopping term for . In the figure we have taken ; but for there is no gap for the two bands closest to zero energy (i.e. the Dirac point).

Figure 3: [color online] Band dispersions near the K points in the bilayer along the direction , with and . Left: bands obtained from the full model in Eq. (6) with , and ; Right: bands obtained from the simplified model in Eq. (10).

ii.3 The Slonsczewki-Weiss-McClure (SWM) model

First we make the observation that the graphene bilayer in the A-B stacking is just the unit cell of graphite that we depict in Fig. 1. Therefore, if the two planes are equivalent much of the symmetry analysis of graphite is also valid for the graphene bilayer. Thus we could alternatively use the SWM for graphite with the proper identification of the parameters. The SWM model for graphite,McClure (1957); Slonczewski and Weiss (1958) is usually written as

(7)

where

(8a)
(8b)
(8c)
(8d)
(8e)
(8f)

Here , and , with being the interplane distance. Typical values of the parameters from the graphite literature are shown in Table 1.

3.16 0.39 -0.02 0.315 0.044 0.038 0.008 -.024
3.12 0.377 -0.020 0.29 0.120 0.0125 0.004 -.0206
Table 1: Values of the SWM parameters for the band structure of graphite. Upper row from Ref. [Brandt et al., 1988] and lower row from Ref. [Chung, 2002].

It is straightforward to show that by identifying and taking , and the matrices in Eq. (6) and Eq. (7) are equivalent up to a unitary transformation. Hence they give rise to identical eigenvalues and band structures. This completes the correspondence between the tight-binding model and the SWM model (see also Refs. [Johnson and Dresselhaus, 1973; Partoens and Peeters, 2006] for a discussion on the connection between the tight-binding parameters and those of SWM.)

The accepted parameters from the graphite literature results in electrons near the K point [] and holes near the H point [] in the BZ as sketched in Fig. 4. These electron and hole pockets are mainly generated by the coupling that in the tight-binding model corresponds to a hopping between the B-atoms of next-nearest planes. Note that this process involves a hopping of a distance as large as .

Figure 4: [color online] Left: graphite lattice; Right: three-dimensional BZ with the symmetry points K and H indicated. The accepted parameters for graphite results in electron pockets near the K points and hole pockets near the H points as sketched in the figure.

Finally, it is interesting to note that at the H-point in the BZ, , and therefore the two planes “decouple” at this point. Furthermore, if one neglects the spectrum is that of massless Dirac fermions just like in the case of graphene. Note that in graphite A and B atoms are different however, and that the term parametrized by , that breaks sublattice symmetry in each plane, opens a gap in the spectrum leading to massive Dirac fermions at the H-point. Since the value of in the literature is quite small, the almost linear massless behavior should be observed by experimental probes that are not sensitive to this small energy scale.

The values of the parameters used in the graphite literature are consistent with a large number of experiments. The most accurate ones are various magneto-transport measurements discussed in Ref. [Brandt et al., 1988]. More recently, angle-resolved photo-emission spectroscopy (ARPES) was used to directly visualize the dispersion of massless Dirac quasi-particles near the H point and massive quasi-particles near the K point in the BZ.Zhou et al. (2006a, b); Ohta et al. (2006)

The band structure of graphite has been calculated and recalculated many times over the years, a recent reference is Ref. [Charlier et al., 1991]. It is also worth to mention that because of the (relatively) large contribution of the non-local van der Waals interaction to the interaction between the layers in graphite, the usual local density approximation or semilocal density approximation schemes are off by an order of magnitude when the binding energy of the planes are calculated and compared with experiments. For a discussion of this topic and a possible remedy, see Ref. [Rydberg et al., 2003].

Iii Simplified electronic band models

In this section we introduce three simplified models that we employ for most of the calculations in this paper. We also show how to obtain an effective two-band model that is sometimes useful.

iii.1 Unbiased bilayer

For the unbiased bilayer, a minimal model includes only the nearest neighbor hopping energies within the planes and the interplane hopping term between A atoms; this leads to a Hamiltonian matrix of the form:

(9)

near the K point in the BZ. Here we write , where and is the appropriate angle. Note that the absolute value of the angle can be changed by a gauge transformation into a phase of the wave functions on the B sublattices. This reflects the rotational symmetry of the model. If one includes the “trigonal distortion” term parametrized by the rotational symmetry is broken and it is necessary to keep track of the absolute value of the angle. From now on in this paper, we most often use units such that for simplicity.

This Hamiltonian has the advantage that it allows for relatively simple calculations. Some of the fine details of the physics might not be accurate but it works as a minimal model and capture most of the important physics. It is important to know the qualitative nature of the terms that are neglected in this approximation, this is discussed later in this section. It is also an interesting toy model as it allows for (approximately) “chiral” particles with mass (i.e., a parabolic spectrum) at low energies.McCann and Fal’ko (2006)

iii.2 Biased graphene bilayer

For the biased bilayer, a minimal model employs Eq. (9) augmented with the bias potential :

(10)

This model was introduced in Refs.  [Guinea et al., 2006; McCann, 2006]. It correctly captures the formation of an electronic gap of size at the K point and the overall features of the bands as can be seen in Fig. 3. Nevertheless, the fine details of the bands close to the band edge are not captured correctly in this simple model; this fact is illustrated in Figs. 5-6. In particular the simple model is cylindrically symmetric; whereas the ”trigonal distortion” breaks this symmetry. Thus the inclusion of results in a ”trihorn” structure for small values of and a weaker modulation of the band edge for larger values of as illustrated in Fig. 5.

Figure 5: [color online] Band dispersions near the band edge (note the energy scale) in the biased graphene bilayer. Left: , Right: . The solid line is the simplified model in Eq. (10) that is cylindrically symmetric. The other lines are along different directions in the BZ for the full model in Eq. (6): (dotted), (dash-dotted), and (dashed). The parameters are the same as in Fig. 3
Figure 6: [color online] Contour plots of the band dispersions near the band edge in the biased graphene bilayer for . Left: full model in Eq. (6), Right: simplified model in Eq. (10). The parameters are the same as in Fig. 3

iii.3 Multilayer graphene

In the graphene multilayer, a minimal model for the bands is again given by Eq. (9) with the understanding that the momentum label also includes the perpendicular direction: . The only change is that we must make the substitution everywhere. Note that this is exactly the -factor appearing in the SWM model discussed in Section II.3. In the following we often use units such that the interplane distance is set to , then – since the unit cell holds two layers – the allowed values of lies in the interval . We note that this band model was used already in the seminal paper by Wallace as a simple model for graphite.Wallace (1947) More recent works on the band structure of few-layer graphene systems include Refs. Latil and Henrard, 2006; Partoens and Peeters, 2006; Guinea et al., 2006; Koshino and Ando, 2007; Nakamura and Hirasawa, .

iii.4 Approximate effective two-band models

There are two main reasons for constructing approximate two-band models. Firstly, on physical grounds the high-energy bands (far away from the Dirac point) should not be very important for the low-energy properties of the system. Secondly, it is sometimes easier to work with instead of matrices. Nonetheless, it is not always a simplification to use the 2-band model when one is studying inhomogeneous systems as it generally leads to two coupled second order differential equations whereas the 4-band model involves four coupled linear differential equations. The matching of the wave functions in the 2-band case then involves both the continuity of the wave function and its derivative whereas in the 4-band model only continuity of the wave function is necessary. We note that the two-band description of the problem is only valid as long the electronic density is low enough, that is, when the Fermi energy is much smaller than . At intermediate to high densities a 4-band model is required in order to obtain the correct physical properties.Kusminskiy et al. (2007)

In this section, we derive the low-energy effective model by doing degenerate second order perturbation theory. The quality of the expansion is good as long as . We first present the general expression for the second-order effective Hamiltonian, thereafter various simplified forms are introduced. Analyses similar to the one here were previously described in Refs. [McCann and Fal’ko, 2006; Nilsson et al., 2006b].

First we introduce the projector matrices () that projects onto (out of) the low-energy subspace of the B atoms. Then we split the Hamiltonian in Eq. (6) according to: , with

(11)
(12)
(13)

Introducing the vectors that to zeroth order only have components in the low energy subspace (i.e., ) and following the standard procedure (see e.g., Ref. [Sakurai, 1994]) for degenerate perturbation theory we arrive at:

(14)

where we explicitly assume that is of the same order as . This expression is correct to second order in . Note that we are doing second order perturbation theory for all of the components of the -matrix in the low-energy subspace. Working to first order in (and then setting ) one obtains for this matrix (taking for brevity):

(15)

Taking leads to an even simpler expression, in particular the effective spectrum becomes:

(16)

That this approximation to the bands is excellent near the band edge for small values of the bias is illustrated in Fig. 7. For larger values of the bias the agreement is less accurate because the assumption of smallness of certain terms in the perturbation expansion is no longer valid.

Figure 7: [color online] Band dispersions near the band edge in the biased graphene bilayer along two directions in the BZ. Left: , Right: . The solid (dash-dotted) line is the effective model in Eq. (16) along (). The dotted (dashed) line is the full model in Eq. (6) along (). The parameters are the same as in Fig. 3. For the different curves are almost not discernible.

Iv Green’s function in the graphene bilayer

As discussed in Section III we use the minimal model Hamiltonian in Eq. (9). We note that the phases can be gauged away by an application of a unitary transformation defined by the matrix

(17)

It is also easy to compute the energy eigenvalues that are given by and . Before we solve for the Green’s functions it is convenient to allow for a local frequency dependent self-energy in the problem. In the general case the self-energies on all of the inequivalent sites in the problem are allowed to be different, and we explicitly introduce the matrix

(18)

to describe this. The Green’s function matrix is then given by the equation

(19)

In the case of the unbiased bilayer the A (B) sites in both of the layers are equivalent and we only need two self-energies: and which are local but we allow for a frequency () dependence. In this case the matrix inversion is simple since it factorizes into two matrices. An explicit form is given by

(20)

where . Here D (ND) stands for diagonal (non-diagonal) in the layer index. The components of the g-matrices are given by

(21a)
(21b)
(21c)

where

(22)

Note that we often suppress the momentum and frequency dependence in the following when no confusion arises. We will come back to the biased case in Section XV.

V Impurity scattering: t-matrix and coherent potential approximation

We are interested in the influence of disorder in the bilayer. To model the impurities we use the standard t-matrix approach and the Coherent Potential Approximation (CPA). The effect of repeated scattering from a single impurity can be encoded in a self-energy which can be computed from:Jones and March (1985)

(23)

Here is a matrix that encodes both the strength and the lattice site of the impurity in question. For example, an impurity on an A1 lattice site of strength at the origin is encoded in Fourier space by the matrix

(24)

implying that the potential is located only on a single site. We have also introduced the quantity:

(25)

which is the local propagator at the impurity site; and in the second step the -sum is to be taken over the whole BZ. The last line is an approximate expression that is obtained by expanding the propagator close to the K points and taking the continuum limit with the introducing of the cutoff . We estimate the cutoff by a Debye approximation that conserves the number of states in the BZ. Then and in units of the cutoff we have (taking ). Due to the special form of the propagator and the impurity potential the self-energy we get from this is diagonal. The result for site dilution (or vacancies) is obtained by taking the limit , so that the electrons are not allowed to enter the site in question. We also introduce a finite density of impurities in the system. To leading order in the impurity concentration the equations for the self-energies then become:

(26a)
(26b)

The explicit form of the propagators in Eq. (21) makes it easy to compute the ’s. The t-matrix result for the self-energies is obtained by using the bare propagators on the right hand side of Eq. (26). In the CPA one assumes that the electrons are moving in an effective medium with recovered translational invariance which in this case is characterized by the local self-energies. To determine what the medium is, one must solve the equations self-consistently with the full propagators on the right hand side of Eq. (26). Because of the simple form of the propagators this is a simple numerical task in the model we are using. To simplify the equations further we assume that . This is a physical assumption since when the self-energies becomes of the order of the cutoff the effective theory breaks down. The self-consistent equations then reads:

(27a)
(27b)

This includes inter-valley scattering in the intermediate states. It is easy to obtain the non-disordered density of states from these equations by taking and (here is a positive infinitesimal) resulting in:

(28a)
(28b)

Observe in particular that the density of states on the A sublattice goes to zero in the limit of zero frequency, this fact is responsible for much of the unconventional physics in the graphene bilayer. In contrast the density of states on the B sublattice is finite at . We discuss how this result is changed with disorder and the solution of Eq. (27) in the following.

v.1 Zero frequency limit

One interesting feature of the CPA equations in Eq. (27) is that it is easy to see that they do not allow for a finite in the limit of . Since by setting the last term in Eq. (27b) must vanish, and this is not possible for finite values of . Then one also must have that there. This implies that the density of states on sublattice A is still zero even within the CPA in the limit . More explicitly, by defining one can show (assuming and ) that and are given asymptotically by:

(29a)
(29b)

and satisfies

(30)

Notice that is exactly the energy scale that is generated by disorder of the same kind in the single layer case.Peres et al. (2006) We have checked that the expressions in Eq. (29) seem to agree with the numerical calculations in the small frequency limit, and the frequency range in which they holds grows with increasing .

v.2 Uncompensated impurity densities

The divergence of the self-energy on the A sublattice in the CPA in the above is due to the fact that there is a perfect compensation between the number of impurities on the two sublattices: . For the more general case where the divergence is not present so that may become finite at . To make comparison with other work on the graphene bilayer it is fruitful to consider another extreme limit where only the B sites are affected by the disorder. Explicitly this means that we take and . The generalization of the CPA equations in Eq. (27) for this case then immediately imply that . In the limit of , is finite, purely imaginary, and given by:

(31)

v.3 Born scattering

Another often studied limit is the one of weak impurities, in particular Koshino and Ando have studied electron transport in the graphene bilayer in this approximation.Koshino and Ando (2006) This is the Born limit and it can be studied using perturbation theory in the strength of the impurities . The leading non-trivial contribution to the self-energies is given by the contribution to second order:

(32a)
(32b)

If one substitutes the bare propagators on the right hand side one finds and at the Dirac point. Thus, to leading order Born scattering is formally equivalent to the previous case with vacancies on only sublattice B exactly at . The frequency range for which the result is valid is different however.

v.4 Self-energy comparisons and the density of states

We compare the self-energies obtained from the t-matrix and the CPA. Within the t-matrix the self-energies are given by

(33a)
(33b)

where the ’s are given in Eq. (28) and

(34a)
(34b)

The results for the self-energies in the two different approximations are shown in Figs. 8 and 9. Note that at least on the scale of the figures the diverges as in the CPA, as discussed above. The solution to the self-consistent equations also does not converge very well when they are pushed close to the limit of . The total DOS on the A-sublattice and B-sublattice is pictured in Fig. 10. Note in particular that the case of closely resemble the non-interacting case except for the new low-energy feature. A possible interpretation of the enhancement of the DOS on the B sublattice near is in terms of the “half-localized” states (meaning they do not decay fast enough to be normalizable at infinity) that have been studied for monolayer graphene.Pereira et al. (2006) Because these states have weight on only one sublattice (the opposite one of the vacancy) the construction in Ref. [Pereira et al., 2006] is valid also in the graphene bilayer when there is a vacancy on one of the A sublattices. For a discussion of the related problem of edge states in bilayer graphene see Ref. [Castro et al., ].

Figure 8: [color online] Self-energies within the t-matrix approximation in the bilayer as a function of the frequency . Left: sublattice A; Right: sublattice B.
Figure 9: [color online] Self-energies within the CPA in the bilayer as a function of the frequency . Left: sublattice A; Right: sublattice B.
Figure 10: [color online] Local density of states on the different sublattices in the bilayer. Left: t-matrix; Right: CPA. Top: sublattice A; Bottom: sublattice B.

v.5 Spectral function

The electron spectral function , which is observable in ARPES experiments, is defined by , so that in our case:

(35)

The spectral function in the plane, calculated within the CPA, is pictured in Fig. 11. As can be seen in the figures the low-energy branch becomes significantly blurred, especially for the higher impurity concentrations. Note also that the gap to the high-energy branch becomes slightly larger as the disorder value increases due to the fact that is not negligible there.

Figure 11: Intensity plot of the spectral function in the plane (normalized by the cutoff) in the bilayer for different impurity concentrations in the CPA approximation. From left to right: , , .

Examples of the momentum distribution curves (MDC’s) and the energy distribution curves (EDC’s) in the disordered graphene bilayer are shown in Figs. 12 and 13. The evolution of the peaks from delta functions to broader peaks with increasing disorder is clear in the figure.

Figure 12: MDC’s in the graphene bilayer. The three panels are for different impurity concentrations and are calculated in the CPA. From the top the energy cuts are at the energies: .0005, .0055, .0105, .0155, .0205, .0255, and .0305 in units of the cutoff . The curves are uniformly displaced for clarity.
Figure 13: EDC’s in the graphene bilayer. The three panels are for different impurity concentrations and are calculated in the CPA. From the top the cuts are at fixed values of the in-plane momentum : .0001, .01, .02, .03, .04, .05, and .06 in units of the cutoff . The curves are uniformly displaced for clarity.

Vi Green’s function and one-particle properties in multilayer graphene

We will use the extension of the bilayer model to the multilayer that we introduced in Section III.3. As discussed there we can immediately use the Hamiltonian in Eq. (9) with the understanding that the momentum label also includes the perpendicular direction: and by substituting everywhere. In particular the Green’s function including the self-energies are again given by the expressions in Eq. (21) and Eq. (22) with the substitution .

vi.1 Self-energies and the density of states

To get the CPA equations we must evaluate the local propagator that is given by the straightforward generalization of Eq. (25) to include an extra sum over :

(36)

Here is the number of unit cells in the perpendicular direction. This extra sum can be transformed into an integral using the relation . It is possible to perform the integrals analytically as we explain in App. B. Using the integrals defined there ( and ) we obtain for

(37a)
(37b)

From these equations one can easily obtain the non-interacting density of states by considering the clean retarded case and take and , from which we get:

(38a)
(38b)

The equivalent expression were previously obtained in Ref. [Guinea et al., 2006] using a different method. The self-energy within the t-matrix is again given by the expression in Eq. (33), with the ’s given by the non-interacting density of states in Eq. (38). The F’s are obtained from the real parts of the non-interacting local propagators:

(39a)
(39b)

The self-energies obtained within the t-matrix are shown in Fig. 14 while those obtained from the CPA are shown in Fig. 15. A comparison between the density of states in the different approximations is shown in Fig. 16. In general the curves are similar to the ones in the bilayer but somewhat smoother.

The behavior of the self-energies at the Dirac point in the multilayer are similar to the case of the bilayer treated in Section V.1V.2 and V.3. We have more to say about this when we discuss the perpendicular transport in the multilayer in Section IX.1.

Figure 14: [color online] Self-energies within the t-matrix in the multilayer as a function of the frequency . Left: sublattice A; Right: sublattice B.
Figure 15: [color online] Self-energies within the CPA in the multilayer as a function of the frequency . Left: sublattice A; Right: sublattice B.
Figure 16: [color online] Local density of states on the different sublattices as a function of the frequency in the multilayer. Top: sublattice A; Bottom: sublattice B.

vi.2 Spectral function

The spectral function for the graphene multilayer is given by a generalization of Eq. (35):

(40)

Given the form of the Green’s function and the CPA self-energies it is straightforward to obtain this quantity. The spectral function is depicted in Fig. 17 for three values of the perpendicular momentum, since the model we use is electron-hole symmetric we only show the results for negative frequencies. We would like to stress that for a large part of the BZ the spectra are reminiscent of the bilayer spectra. At the edges of the BZ however, where , the spectrum is that of massless Dirac fermions.

Figure 17: Intensity plot of the electron spectral function in the plane for the multilayer for different disorder values and different values of in the CPA. From left to right: , , . From top to bottom: , , .

Examples of the momentum distribution curves (MDC’s) and the energy distribution curves (EDC’s) in the disordered graphene multilayer are shown in Figs. 18 and 19 for two different values of . The evolution of the peaks from delta functions to broader peaks with increasing disorder is clearly seen. One can also note that the influence of the impurities is more severe close to the H point of the BZ since the overlap of the peaks is larger there. The reason for this is that the particles there are dispersing linearly leading to peaks that are closer together than for particles with a parabolic dispersion.

Figure 18: [color online] MDC’s in the graphene multilayer for two values of the perpendicular momentum: (i.e. at the K point, parabolic dispersion, dashed line) and (i.e. at the H point, linear dispersion, solid line). The three panels are for different values of the density of impurities in the CPA. From the top the energy cuts are at the energies: .0005, .0055, .0105, .0155, .0205, .0255, and .0305 in units of the cutoff . The curves are uniformly displaced for clarity.
Figure 19: [color online] EDC’s in the graphene multilayer for two values of the perpendicular momentum: (i.e. between the K and H point, parabolic dispersion, dashed line) and (i.e. at the H point, linear dispersion, solid line). The three panels are for different values of the density of impurities in the CPA. From the top the cuts are at fixed values of the in-plane momentum : .0001, .01, .02, .03, .04, .05, and .06 in units of the cutoff . The curves are uniformly displaced for clarity.

Vii Electron transport in bilayer and multilayer graphene

Having worked out the self-energies in the previous sections we now turn to linear response (Kubo formula) to study electron transport. We saw in Sections V and VI that the low-energy states are mainly residing on the B sublattice. Nevertheless, electron transport coming from nearest neighbor hopping must go over the A sublattice. This is particularly important for the case of perpendicular transport, since in the simple model that we are using, hopping comes exclusively from states with weight on the A sublattice.

This feature implies that although the total density of states at half-filling is finite, because the density of states on the A sublattice goes to zero as the Dirac point is approached, the in-plane and out-of-plane transport properties are unconventional. The main purpose of this section and the two following sections is to show how this comes about in detail through concrete calculations. We calculate conductivities (or optical response) parallel and perpendicular to the planes in both bilayer graphene and multilayer graphene. The resulting conductivities are very anisotropic and we find a universal nonzero minimum value for the in-plane DC conductivity as a function of the chemical potential.

vii.1 Conductivity formulas

To calculate the conductivity we use the Kubo formula.Mahan (2000) We only consider the homogeneous () response, but we investigate both the temperature dependence and the frequency dependence of the various conductivities.

The conductivity is computed from the Kubo formula:

(41)

Here is the area of the system, is the current operator of interest, and the appropriate current-current correlation function. A contribution to from a term of the form where and can be shown by the usual methods to give a contribution to the real part of the conductivity of the form:

(42)

Here the imaginary parts only involve the frequency part and not the angular (spatial) parts of the propagators. In terms of the expression in Eq. (20) this imply that the imaginary parts involves and but not the spatial information encoded in and . With the inclusion of the two spin projections and two valleys we get (putting back and extracting the electric charge from the current operators):

(43)

Here is the Fermi distribution function. We have also introduced the kernel that for the case of the operators above becomes:

(44)

Thus the contribution to the in-plane DC conductivity at zero temperature is