Electronic transport in bilayer graphene

We present theoretical studies on the transport properties and localization effects of bilayer graphene. We calculate the conductivity by using the effective mass model with the self-consistent Born approximation, in the presence and absence of an energy gap opened by the interlayer asymmetry. We find that, in the absence of the gap, the minimum conductivity approaches the universal value by increasing the disorder potential, and the value is robust in the strong disorder regime where mixing with high-energy states is considerable. The gap-opening suppresses the conductivity over a wide energy range, even in the region away from the gap.We also study the localization effects in the vicinity of zero energy in bilayer graphene. We find that the states are all localized in the absence of the gap, while the gap-opening causes a phase transition analogous to the quantum Hall transition, which is accompanied by electron delocalization.

In this paper, we present theoretical studies on the transport properties of bilayer graphene from two different viewpoints. First, we calculate the conductivity within the self-consistent Born approximation (SCBA). The SCBA successfully works in describing the disorder effects in various kinds of system and has been applied to monolayer graphene [36]- [39]. For bilayer graphene, we previously calculated the SCBA conductivity in the absence of an energy gap, using the effective two-band model valid in the vicinity of zero energy [14]. Here, we extend the analysis to the four-band model that is valid in a broader energy range, and study the effects of higher bands as well as the asymmetry gap. In the second part, we study the electron localization properties in bilayer graphene [40]. Generally, the localization effect gives a quantum correction to the conductivity which is not taken in the SCBA, and such effects are investigated for bilayer graphene [20,21] as well as monolayer [41,42]. Here, we numerically calculate the energy levels in disordered bilayer graphene of finite system sizes, to estimate the localization length from its scaling behavior. We formulate an effective mass model of bilayer graphene in section 2. We present the calculation of the SCBA conductivity in section 3, the study of the localization effects in section 4 and the conclusion in section 5.

Effective mass model of bilayer graphene
Bilayer graphene is composed of a pair of hexagonal networks of carbon atoms, which include A 1 and B 1 atoms on layer 1 and A 2 and B 2 on layer 2. As shown in figure 1, the two layers are arranged in Bernal stacking, where A 2 atoms are located directly below B 1 atoms [10]. The lattice constant within a layer is given by a ≈ 0.246 nm and the layer spacing by d ≈ 0.334 nm. We adopt the Hamiltonian based on the effective mass model for bulk graphite [5,6]. We include parameters γ 0 , γ 1 and γ 3 , where γ 0 represents the intralayer coupling A j ↔ B j , γ 1 the interlayer coupling B 1 ↔ A 2 , and γ 3 another interlayer coupling between A 1 and B 2 . For typical values, we quote [43] γ 0 ≈ 3.16 eV, γ 1 ≈ 0.39 eV and γ 3 ≈ 0.315 eV, which will be used in the following numerical calculations throughout. We neglect the parameter γ 4 for the couplings between A 1 and A 2 , and B 1 and B 2 , which have a relatively minor role in the electronic structure [6,44].
The low-energy spectrum is given by the states around the K and K points in the Brillouin zone. Let |A j and |B j be the Bloch functions at the K point, corresponding to the A and B γ 0 Figure 1. Atomic structure of bilayer graphene. +δ and −δ represent the potential of the top and bottom layers, respectively. sublattices, respectively, of layer j. The Hamiltonian for the basis (|A 1 , |B 1 , |A 2 , |B 2 ) is given by [10] where p ± = p x ± i p y with p = −ih∇, δ describes a possible potential difference between layers 1 and 2, v = √ 3γ 0 a/2 is the band velocity of monolayer graphene and v 3 = √ 3γ 3 a/2 is another velocity associated with γ 3 . The effective Hamiltonian for K can be obtained by exchanging p + and p − , thereby giving the equivalent spectrum.
Since a unit cell has four atoms, there are four electronic bands, as shown in figure 2. When δ = 0 and γ 3 = 0, the eigenenergies of equation (1) are given by [10,11] ε µ,s (p) = s where µ = ±, s = ± and p = p 2 x + p 2 y . The index µ = − gives a pair of bands closer to zero energies, and µ = + another pair repelled away by approximately ±γ 1 . In each pair, s = + and − represent the electron and hole branches, respectively. The higher band states of µ = + correspond to the dimer states formed mainly from the sites B 1 and A 2 , which are coupled by γ 1 , while the lower band states of µ = − come from the remaining sites A 1 and B 2 . The electron and hole bands of µ = − touch at zero energy in δ = 0, where the dispersion is approximately quadratic with an effective mass [10] m * = γ 1 2v 2 .
(3) Figure 2 shows the energy band dispersion given by equation (1) with several δs. A finite δ opens a gap around zero energy, as it gives a difference of energy between A 1 and B 2 . In the absence of γ 3 , the dispersion has a 'Mexican hat' structure where the band bottom is isotropically located at the off-center momentum [11]. The parameter γ 3 gives a trigonal warping effect in which the isotropic band dispersion is distorted in a trigonal symmetry [10,14]. At δ = 0, it produces a fine structure in the vicinity of zero energy, where the equi-energetic line splits into four pockets [10,14], as shown in the contour plots in figure 3(a). The separation occurs in small energy |ε| < ε trig with When δ and γ 3 coexist, the Mexican hat structure is trigonally warped into three pockets [45]. Figure 3(b) shows the contour plot in δ = 0.2γ 1 , where the trigonal pockets are largely enhanced as compared to δ = 0 in figure 3(a). Note here that the axes are differently scaled between (a) and (b). Figure 2(b) shows the density of states (DOS) for several values of δ. At δ = 0, the DOS linearly increases with |ε| apart from the fine structure |ε| ε trig . When this structure is smeared 5 out, the DOS at zero energy is roughly given by g v g s ρ 0 , with ρ 0 = m * 2πh 2 (5) and g v = g s = 2 being the valley and spin degeneracies, respectively. The steps at ε = ±γ 1 correspond to the higher band edges. In the presence of δ, we have a singular enhancement of the DOS around the edge of the energy gap, corresponding to the flattened band bottom shown in figure 2(a). The Hamiltonian of equation (1) can be reduced for the low-energy bands with basis |A 1 , |B 2 as [10] This effective Hamiltonian describes the low-energy part of the dispersion quite well including the fine structure caused by the trigonal warping. The approximation fails roughly in |ε| γ 1 /4 ≈ 100 meV, where vp becomes of the order of γ 1 . When the trigonal warping is neglected, the DOS of this model becomes constant at g v g s ρ 0 independently of energy.

Disorder potential
For the disorder potential, we assume that the length scale is much longer than the atomic scale and neglect the intervalley scattering. We also assume that the disorder length scale is shorter than the typical wavelength of 2π/k with k being the wavenumber from K or K points, so as to be modeled as a short-ranged potential within the valley-decoupled Hamiltonian [14,36]. We suppose that there is no difference between the disorder potentials of layers 1 and 2. The disorder term is then expressed as We assume an equal amount of positive and negative scatterers u i = ±u and a total density per unit area n imp . The total Hamiltonian for the K valley, H = H 0 + U , belongs to the orthogonal class only when both δ and γ 3 are zero; otherwise it becomes unitary class [20]. To show this, we define the effective time-reversal operator T as where ψ is a wavefunction and S is a matrix, At δ = γ 3 = 0, we can immediately show that HT = T H and T 2 = 1, i.e. the Hamiltonian H is orthogonal. It should be noted that the operation T is an effective time-reversal symmetry within a single valley, and distinguished from the real time-reversal symmetry connecting K and K that always exists at any values of δ and γ 3 .

Single-valley Hall conductivity
In bilayer graphene with a nonzero δ, interestingly, the Hamiltonian at a single valley has a non-zero Hall conductivity even in a zero magnetic field. The Hall conductivity is estimated by the Kubo formula, where S is the area of the system, v x and v y are the velocity operators, η is the positive infinitesimal, f (ε) is the Fermi distribution function, and |α and ε α describe the eigenstate and the eigenenergy of the system, respectively. We estimate the Hall conductivity in the Hamiltonian H (eff) 0 of equation (23) without disorder potential and at zero temperature. When γ 3 is neglected, it reads with sgn(x) = x/|x|. In the absence of δ, we always have σ K x y = 0 in all energies. We note that σ K x y in the energy gap (|ε F | < |δ|) is quantized into different values in the negative and positive δs. We can show that the quantized value at the energy gap is never changed in the presence of γ 3 . The result for the K point is given by σ K x y = −σ K x y , so that the net Hall conductivity is exactly zero as it should be. The single-valley Hall conductivity can never be directly observed, but we will see in section 4 that it strongly influences the electron localization property.

Self-consistent Born approximation
We estimate the conductivity with the SCBA. We define Green's function as Green's function averaged over the impurity configurations satisfies Dyson's equation where is the average over configurations of the disorder potential, α is an eigenstate of H 0 , and G (0) α = (ε − ε α ) −1 , with ε α being the eigenenergy of the state α in H 0 . In the SCBA, the self-energy is given by [36] α,α (ε) = 7 Equations (13) and (14) need to be solved self-consistently. The conductivity is calculated by the Kubo formula as Re where G R = G(ε + i0) and G A = G(ε − i0) are the retarded and the advanced Green's functions, respectively and v x = ∂H 0 /∂ p x is the velocity operator. This is transformed as Re withṽ RA In the SCBA,ṽ x should be calculated in the ladder approximation. The SCBA conductivity was previously calculated in the reduced 2 × 2 model of equation (6) in the zero-gap case [14]. The energy broadening due to the disorder potential at zero energy is characterized by When ε trig , i.e. the disorder is strong enough to smear out the fine band structure of the trigonal warping, we can approximately solve the SCBA equation in an analytic form. The selfenergy equation (14) becomes diagonal with respect to the index α and becomes a constant, The conductivity at the Fermi energy ε is written as The third term in the square bracket arises from the vertex correction of equation (17) due to the trigonal warping. In high energies |ε| , σ approximates as which increases linearly with energy. The value at zero energy becomes In the strong disorder limit ε trig / → 0, the correction arising from the trigonal warping vanishes so that the conductivity approaches the universal value g v g s e 2 /(π 2h ). As the 2 × 2 model works well in the vicinity of energy, we naively expect that the conductivity estimated above would be no longer valid in the strong disorder regime where mixing to the higher energy part is considerable. To check this, we go back to the original 4 × 4 Hamiltonian of equation (1) and solve the self-consistent equations (13) and (14) with a numerical iteration. In practice, we first replace G with G (0) in equation (14), and obtain the selfenergy (1) . Then we substitute (1) for in equation (13) in order to obtain G (1) as a solution of G . We iterate the process until G (i+1) and G (i) are equal within the numerical precision. We introduce the momentum-space cut-off at vp = 10γ 1 , of the order of the graphene band width, to avoid the divergence of the self-energy [36]. Figure 4(a) shows the DOS calculated for several disorder strengths at δ = 0. We see that the states of the larger energy part are pushed toward zero energy by the disorder potential, leading to an enhancement of the DOS around zero energy. This is in contrast to the 2 × 2 model, where the DOS is almost constant at g v g s ρ 0 regardless of the disorder strength.
We show the conductivity σ (ε) and its expanded plot in figures 4(b) and (c), respectively. For comparison, we append in (c) the results for the 2 × 2 model expressed in equation (20) as dashed curves. Significantly, the zero-energy conductivity barely shifts from that of the 2 × 2 model (equation (20)), and actually approaches the universal value g v g s e 2 /(π 2h ) as increases. It is surprising that this stands in the stronger disorders such as /γ 1 = 0.15, even though the DOS at zero energy becomes nearly twice as large as in the 2 × 2 model. In the strong disorder limit, γ 1 such that the interlayer coupling γ 1 is negligible, it is expected that the minimum conductivity is just twice as large as that of the monolayer graphene [36], and is again 2 × g v g s e 2 /(2π 2h ). In higher energy |ε| > , the conductivity increases linearly with |ε| in qualitative agreement with equation (20) of the 2 × 2 model, while the gradient is generally steeper. The conductivity has a small dip at the higher band edge around |ε| ∼ γ 1 , and this is because the frequency of electron scattering is strongly enhanced by introducing the higher band states. Figures 5(a) and (b) show DOS and conductivity, respectively, calculated for several values of δ with a fixed disorder strength /γ 1 = 0.03. In increasing δ, the conductivity decreases around the gap region as is naturally expected. In large δs, we observe a considerable suppression of the conductivity in the overall energy region even when away from the gap. There the band velocity is strongly suppressed due to a relatively flat band bottom shown in figure 2(a).

Localization effects
We study in this section the localization property of bilayer graphene [40], neglected in the SCBA approach. Here we concentrate on the lower electron and hole bands (µ = −) and use the 2 × 2 Hamiltonian matrix shown in equation (23). Here we neglect γ 3 and assume that δ and vp are much smaller than γ 1 . The Hamiltonian of equation (6) is even simplified by neglecting the terms more than the third order in δ/γ 1 and vp/γ 1 as which is valid around zero energy in the small gap regime. For the disorder potential we use We numerically calculate the electronic states in the disorder potential by exactly diagonalizing the total Hamiltonian H (eff) 0 + U (eff) . We consider a finite square system with L × L and impose a boundary condition with phase factors exp(iφ x ) and exp(iφ y ) for xand y-directions, respectively. We introduce the unit wavenumber k 0 ash 2 k 2 0 /(2m * ) = /2 and the unit length λ 0 = 2π/k 0 . To make the matrix finite, we introduce the k-space cut-off k c = 6k 0 , which corresponds to ε c = 18 . From the obtained eigenenergies, we calculate the Thouless number g, which is the ratio of the shift δ E of each energy level due to the change of the boundary condition, to the level spacing [L 2 ρ] −1 with ρ being the DOS per unit area [46]. We estimate the energy shift by δ E = π |∂ 2 E(φ)/∂φ 2 | , where E(φ) represents the eigenenergy as a function of the boundary phase φ = φ x with fixed φ y , and represents the averaging over different levels in a small energy region around the energy ε in question. The localization length L loc is estimated by fitting the results to g(L) ∝ exp(−L/L loc ). When g 1, g is approximately related to the longitudinal conductivity σ x x by σ x x ∼ (e 2 /h)g [46]. We also calculate the Hall conductivity in disordered systems by substituting the eigenstates of the disordered system into the Kubo formula, equation (10). For every quantity, we take an average over a number of samples with different configurations of the disorder potential and boundary phase factors φ x and φ y . In figure 6(a), we plot the Thouless number g at zero energy and the inverse localization length 1/L loc , as a function of δ. At δ = 0, g drops exponentially as the system size expands, and thus the localization length becomes finite. This is in striking contrast to monolayer graphene where localization never takes place in the single-valley Hamiltonian [41,49]. As δ starts from zero, the localization becomes weaker (i.e. L loc larger), and at δ ∼ 0.8 , significantly, g becomes independent of the system size and thus the localization length diverges. At an even larger δ, the localization length goes back to a finite value and rapidly shrinks with increasing δ.
The corresponding plots for the single-valley Hall conductivity σ K x y at zero energy are shown in figure 6(b). σ K x y becomes completely 0 at δ = 0 as well as at the ideal limit. As δ increases, it generally goes up and appears to approach e 2 / h, the mid-gap value in the ideal limit in equation (11). At each single δ, we notice that σ K x y moves in a specific direction as the system size L expands. In δ 0.8 , σ K x y decreases as L increases, while in the opposite side δ 0.8 , in contrast, σ K x y increases as L increases. To relate the results of the Hall conductivity and the localization length above, we use an analogy of the conventional integer quantum Hall effect [47]. In a two-dimensional electron system in a magnetic field, the electronic states are quantized into Landau levels, which are generally broadened in energy by the disorder potential. In each broadened Landau level, the localization length diverges only at a single energy E c around the level center, while the states are localized otherwise. The Hall conductivity is quantized in different integers (in units of e 2 / h) between the regions E F < E c and E F > E c , and it jumps from one to another when E F crosses E c . The Hall conductivity never changes when E F stays in the localized region, and when it changes, it is always accompanied by the divergence of the localization length.
In the present system, we expect a similar quantum Hall transition (i.e. a jump in the quantized Hall conductivity) somewhere in the process of the gap opening, since σ K x y at zero energy actually transfers from zero to e 2 / h. In finite-size calculations, the Hall conductivity averaged over disorder configurations generally takes a non-integer value, while it gradually approaches an integer with increasing system size [47,48]. In figure 6(b), we naturally expect that the region of decreasing σ K x y becomes the Hall plateau with σ K x y = 0, and that of increasing σ K x y becomes another plateau with e 2 / h, in the infinite size limit. The point where the σ K x y switches its moving direction, δ ∼ 0.8 , is regarded as the critical point that separates the different Hall phases σ K x y = 0 and 1. Indeed this point coincides well with the divergence of the localization length L loc in figure 6(a), as expected from the analogy to the integer quantum Hall effect. This is strong evidence that the electron delocalization is associated with the change of the single-valley Hall conductivity.
To see the behaviors at other energies, we show in figure 7(a) plots of the Hall conductivities as a function of energy, at several δs. In δ 0.8, we have two critical energies off ε = 0, which separate different scaling behaviors of σ K x y (i.e. increasing or decreasing as L increases). Figure 7(b) is the phase diagram based on figure 7(a). Here each phase corresponds to the Hall plateau with the quantized σ K x y . We determine the phase boundary by taking points where σ K x y changes its moving direction. The phase diagram is symmetric with respect to δ = 0, while the sign of the Hall conductivity is opposite between negative and positive δs.
While the actual Hall conductivity is never directly observed, it would be possible to observe the divergence of the localization length in a temperature (T ) dependence of the conductivity σ x x , as in the conventional integer quantum Hall effect [50]. The trigonal warping effect due to γ 3 , which was neglected in this section, generally suppresses the electron localization at δ = 0, since it changes the symmetry class from orthogonal to unitary [20]. Nevertheless, since the jump of σ K x y in opening δ takes place regardless of γ 3 , we expect that the plateau transition still exists, while the localization length around δ = 0 would be longer than the present results. We leave the estimation of the localization length in the trigonal warping for future studies.

Conclusion
We studied the transport properties and the localization effect in bilayer graphene. We calculated the conductivity using the SCBA with a 4 × 4 effective mass model, and found that the minimum conductivity approximates the universal value g v g s e 2 /(π 2h ) apart from the correction due to trigonal warping. Compared to the 2 × 2 model [14], the conductivity at zero energy barely changes even when mixing with the higher band is considerable, while the value at higher energy generally becomes larger. We also studied the effect of the energy gap, and found that flattening of the band bottom leads to conductivity suppression over a wide energy range. In the latter part, we studied the localization properties. We found that the gap opening causes electron delocalization which can be interpreted as an analog of the quantum Hall transition at each single valley. This offers an example of quantum Hall physics that is purely controlled by the external electric field without any use of magnetic fields, although the Hall current cannot be directly measured.