Gaps and tails in graphene and graphane

We study the density of states in monolayer and bilayer graphene in the presence of a random potential that breaks sublattice symmetries. While a uniform symmetry-breaking potential opens a uniform gap, a random symmetry-breaking potential also creates tails in the density of states. The latter can close the gap again, preventing the system to become an insulator. However, for a sufficiently large gap the tails contain localized states with nonzero density of states. These localized states allow the system to conduct at nonzero temperature via variable-range hopping. This result is in agreement with recent experimental observations in graphane by Elias {\it et al.}.


Introduction
Graphene is a single sheet of carbon atoms that is forming a honeycomb lattice. A graphene monolayer as well as a stack of two graphene sheets (i.e. a graphene bilayer) are semimetals with remarkably good conducting properties [1,2,3]. These materials have been experimentally realized with external gates, which allow for a continuous change in the charge-carrier density. There exists a non-zero minimal conductivity at the charge neutrality point. Its value is very robust and almost unaffected by disorder or thermal fluctuations [3,4,5,6].
Many potential applications of graphene require an electronic gap to switch between conducting and insulating states. A successful step in this direction has been achieved by recent experiments with hydrogenated graphene (graphane) [7] and with gated bilayer graphene [8,9,10]. These experiments take advantage of the fact that the breaking of a discrete symmetry of the lattice system opens a gap in the electronic spectrum at the Fermi energy. In the case of monolayer graphene (MLG), a staggered potential that depends on the sublattice of the honeycomb lattice plays the role of such symmetrybreaking potential (SBP). For bilayer graphene (BLG) a gate potential that distinguishes between the two graphene layers plays a similar role.
With these opportunities one enters a new field in graphene, where one can switch between conducting and insulating regimes of a two-dimensional material, either by a chemical process (e.g. oxidation or hydrogenation) or by applying an external electric field [11].
The opening of a gap can be observed experimentally either by a direct measurement of the density of states (e.g., by scanning tunneling microscopy [12]) or indirectly by measuring transport properties. In the gapless case we observe a metallic conductivity σ ∝ ρD, where D is the diffusion coefficient (which is proportional to the scattering time) and ρ is the density of states (DOS). This gives typically a conductivity of the order of e 2 /h. The gapped case, on the other hand, has a strongly temperature-dependent conductivity due to thermal activation of charge carriers [13] σ(T ) = σ 0 e −T 0 /T (1) with some characteristic temperature scale T 0 which depends on the underlying model. A different behavior was found experimentally in the insulating phase of graphane [7]: which is known as 2D variable-range hopping [14]. This behavior indicates the existence of well-separated localized states, even at the charge-neutrality point, where the parameter T 0 depends on the DOS at the Fermi energy E F as T 0 ∝ 1/ρ(E F ). The experimental observation of a metal-insulator transition in graphane raises two questions: (i) what are the details that describe the opening of a gap and (ii) what is the DOS in the insulating phase? In this paper we will focus on the mechanism of the gap opening due to a SBP in MLG and BLG. It is crucial for our study that the SBP is not uniform in the realistic two-dimensional material. One reason for the latter is the fact that graphene is not flat but forms ripples [15,16,17]. Another reason is the incomplete coverage of a graphene layer with hydrogen atoms in the case of graphane [7]. The spatially fluctuating SBP leads to interesting effects, including a second-order phase transition due to spontaneous breaking of a discrete symmetry and the formation of Lifshitz tails.

Model
Quasiparticles in MLG or in BLG are described in tight-binding approximation by a nearest-neighbor hopping Hamiltonian where c † r (c r ) are fermionic creation (annihilation) operators at lattice site r. The underlying lattice structure is either a honeycomb lattice (MLG) or two honeycomb lattices with Bernal stacking (BLG) [11,18]. We have an intralayer hopping rate t and an interlayer hopping rate t ⊥ for BLG. There are different forms of the potential V r , depending on whether we consider MLG or BLG. Here we begin with potentials that are uniform on each sublattice, whereas random fluctuations are considered in subsection 2.4.

MLG
V r is a staggered potential with V r = m on sublattice A and V r = −m on sublattice B. This potential obviously breaks the sublattice symmetry of MLG. Such a staggered potential can be the result of chemical absorption of non-carbon atoms in MLG (e.g. oxygen or hydrogen [7]). A consequence of the symmetry breaking is the formation of a gap ∆ g = m: The spectrum of MLG consists of two bands with dispersion where for lattice spacing a = 1.

BLG
V r is a biased gate potential that is V r = m (V r = −m) on the upper (lower) graphene sheet. The potential in BLG has been realized as an external gate voltage, applied to the two layers of BLG [8]. The spectrum of BLG consists of four bands [11] with two low-energy bands where ǫ k is the monolayer dispersion of Eq. (5), and two high-energy bands The spectrum of the low-energy bands has nodes for m = 0 where E − k (0) vanishes in a (k − K) 2 manner, where K is the position of the nodes, which are the same as those of a single layer. For small m ≪ t ⊥ , a mexican hat structure develops around k = K, with local extremum in the low-energy band at E − k (m) = ±m, and a global minimum/maximum in the upper/lower low energy band at E − k (m) = mt ⊥ / t 2 ⊥ + 4m 2 . For small gating potential V r = ±m we can expand E − k (m) under the square root near the nodes and get t ⊥ apparently reduces the gap. Very close to the nodes we can approximate the factor in front of m 2 by 1 and obtain an expression similar to the dispersion of MLG: Here we notice the absence of the mexican hat structure in this approximation. The resulting spectra for MLG and BLG are shown in Fig. 1.

Low-energy approximation
The two bands in MLG and the two low-energy bands in BLG represent a spinor-1/2 wave function. This allows us to expand the corresponding Hamiltonian in terms of Pauli matrices σ j as Near each node the coefficients h j read in low-energy approximation [19] where (∇ 1 , ∇ 2 ) is the 2D gradient.

Random fluctuations
In a realistic situation the potential V r is not uniform, neither in MLG nor in BLG, as discussed in the Introduction. As a result, electrons experience a randomly varying potential V r along each graphene sheet, and m in the Hamiltonian of Eq. (9) becomes a random variable in space as well. For BLG it is assumed that the gate voltage is adjusted at the charge-neutrality point such that in average m r is exactly antisymmetric with respect to the two layers: At first glance, the Hamiltonian in Eq. (3) is a standard hopping Hamiltonian with random potential V r . This is a model frequently used to study the generic case of Anderson localization [20]. The dispersion, however, is special in the case of graphene due to the honeycomb lattice: at low energies it consists of two nodes (or valleys) K and K ′ [17,19]. It is assumed here that randomness scatters only at small momentum such that intervalley scattering, which requires large momentum at least near the nodal points (NP), is not relevant and can be treated as a perturbation. Then each valley contributes separately to the DOS, and the contribution of the two valleys to the DOS ρ is additive: ρ = ρ K + ρ K ′ . This allows us to consider the low-energy Hamiltonian in Eqs. (9), (10), even in the presence of randomness for each valley separately. Within this approximation the term m r is a random variable with mean value m r m =m and variance (m r −m)(m r ′ −m) m = gδ r,r ′ . The following analytic calculations will be based entirely on the Hamiltonian of Eqs. (9),(10) and the numerical calculations on the lattice Hamiltonian of Eq. (3). In particular, the average Hamiltonian H m can be diagonalized by Fourier transformation and is for MLG with eigenvalues E k = ± √m 2 + k 2 . For BGL the average Hamiltonian is with eigenvalues E k = ± √m 2 + k 4 .

Symmetries
Low-energy properties are controlled by the symmetry of the Hamiltonian and of the corresponding one-particle Green's function G(iǫ) = (H + iǫ) −1 . In the absence of sublattice-symmetry breaking (i.e. for m = 0), the Hamiltonian H = h 1 σ 1 + h 2 σ 2 has a continuous chiral symmetry with a continuous parameter α, since H anticommutes with σ 3 . The term mσ 3 breaks the continuous chiral symmetry. However, the behavior under transposition h T j = −h j for MLG and h T j = h j for BLG in Eq. (10) provides a discrete symmetry: where n = 1 for MLG and n = 2 for BLG. This symmetry is broken for the oneparticle Green's function G(iǫ) by the iǫ term. To see whether or not the symmetry is restored in the limit ǫ → 0, the difference of G(iǫ) and the transformed Green's function −σ n G T (iǫ)σ n must be evaluated: For the diagonal elements this is the DOS at the NP ρ(E = 0) ≡ ρ 0 in the limit ǫ → 0. Thus the order parameter for spontaneous symmetry breaking is ρ 0 . According to the theory of phase transitions, the transition from a nonzero ρ 0 (spontaneously broken symmetry) to ρ 0 = 0 (symmetric phase) is a second-order phase transition, and should be accompanied by a divergent correlation length at the transition point. Since our symmetry is discrete, such a phase transition can exists in d = 2 and should be of Ising type. A calculation, using the SCBA of ρ 0 , gives indeed a second-order transition at the point where ρ 0 vanishes with a divergent correlation length ξ for the DOS fluctuations . Whether or not this transition is an artefact of the SCBA or represents a physical effect due to the appearence of two types of spectra (localized for vanishing SCBA-DOS and delocalized for nonzero SCBA-DOS) is not obvious here and requires further studies.

Density of states
Our focus in the subsequent calculation is on the DOS of MLG and BLG. In the absence of disorder, the DOS of 2D Dirac fermions opens a gap ∆ ∝m as soon as a nonzero termm appears in the Hamiltonian of Eq. (9), since the low-energy dispersion is E k = ± √m 2 + k 2 for MLG and E k = ± √m 2 + k 4 for BLG, respectively (cf Fig. 2). Here we evaluate the DOS of MLG and BLG in the presence of a uniform gap. Given the energy spectrum, the DOS is defined as By using the MLG dispersion, this reduces to where Θ(x) is the Heaviside function. For BLG, this gives  which are shown in Fig. 2. By retaining the full low-energy spectrum for BLG, E − k , the DOS can still be evaluated in closed form, with the result In the limit of t ⊥ ≫ (E, m), this reduces to Eq. (18) after dividing it by t ⊥ , which was set to 1 in the low-energy approximation, and the DOS saturates to a constant value after the initial divergence. For finite t ⊥ , however, the Dirac nature of the spectrum appears again, and the high energy DOS increases linearly even for the BLG, similarly to the MLG case. For m = 0, and E ≪ t ⊥ , this lengthy expression gives An interesting question, from the theoretical as well as from the experimental point of view, appears here: What is the effect of random fluctuations aroundm? Previous calculations, based on the self-consistent Born approximation (SCBA), have revealed that those fluctuations can close the gap again, even for an average SBP termm = 0 [22]. Only ifm exceeds a critical value m c (which depends on the strength of the fluctuations), an open gap was found in these calculations. This describes a special transition from metallic to insulating behavior. In particular, the DOS at the Dirac point ρ 0 vanishes withm like a power law The exponent 1/2 of the power law is probably an artefact of the SCBA, similar to the critical exponent in mean-field approximations.

Self-consistent Born approximation
The average one-particle Green's function can be calculated from the average Hamiltonian H m by employing the self-consistent Born approximation (SCBA) [23,24,25] The SCBA is also known as the self-consistent non-crossing approximation in the Kondo and superconducting community. The self-energy Σ is a 2 × 2 tensor due to the spinor structure of the quasiparticles: Σ = −(iησ 0 + m s σ 3 )/2. Scattering by the random SBP produces an imaginary part of the self-energy η (i.e. a one-particle scattering rate) and a shift m s of the average SBPm (i.e.,m → m ′ ≡m + m s ). Σ is determined by the self-consistent equation The symmetry in Eq. (14) implies that with Σ also is a solution (i.e. m s → −m s creates a second solution).
The average DOS at the NP is proportional to the scattering rate: ρ 0 = η/2gπ. This reflects that scattering by the random SBP creates a nonzero DOS at the NP if η > 0. Now we assume that the parameters η and m s are uniform in space. Then Eq. (23) can be written in terms of two equations, one for the one-particle scattering rate η and another for the shift of the SBP m s , as η = gIη, m s = −mgI/(1 + gI) .
I is a function ofm and η and also depends on the Hamiltonian. For MLG it reads with momentum cutoff λ and for BLG A nonzero solution η requires gI = 1 in the first part of Eq. (25), such that m s = −m/2 from the second part. Since the integrals I are monotonically decreasing functions for largem, a real solution with gI = 1 exists only for |m| ≤ m c . For both, MLG and BLG, the solutions read where the model dependence enters only through the critical average SBP m c : m c is much bigger for BGL, a result which indicates that the effect of disorder is much stronger in BLG. This is also reflected by the scattering rate atm = 0 which is η = m c /2. A central assumption of the SCBA is a uniform self-energy Σ. The imaginary part of Σ is the scattering rate η, created by the random fluctuations. Therefore, a uniform η means that effectively random fluctuations are densely filling the lattice. If the distribution of the fluctuations is too dilute, however, there is no uniform nonzero solution of Eq. (23). Nevertheless, a dilute distribution can still create a nonzero DOS, as we will discuss in the following: we study contributions to the DOS due to rare events, leading to Lifshitz tails.

Lifshitz tails
In the system with uniform SBP the gap can be destroyed locally by a local change of the SBP m → m + δm r due to the creation of a bound state. We start with a translationalinvariant system and add δm r on site r. To evaluate the corresponding DOS from the Green's function G = (H + iǫ + δmσ 3 ) −1 , using the Green's function G 0 = (H + iǫ) −1 with uniform m, we employ the lattice version of the Lippmann-Schwinger equation [26] with the 2 × 2 scattering matrix In the case of MLG we have (remark: the DOS of BLG has the same form.) Then the imaginary part of the Green's function reads Thus the DOS is the sum of two Dirac delta peaks The Dirac delta peak appears with probability ∝ exp(−(g 0 ± g 3 ) 2 /g) for a Gaussian distribution. This calculation can easily be generalized to δm r on a set of several sites r [26]. Then the appearance of the several such Dirac delta peaks decreases exponentially. Moreover, these contributions are local and form localized states. For stronger fluctuations δm r (i.e., for increasing g) the localized states can start to overlap. This is a quantum analogue of classical percolation. The localized states in the Lifshitz tails can be taken into account by a generalization of the SCBA to non-uniform self-energies. The main idea is to search for space-dependent solutions Σ r of Eq. (23). In general, this is a diffult problem. However, we have found that this problem simplifies essentially when we study it in terms of a 1/m expansion. Using a Gaussian distribution, this method gives Lifshitz tails of the form [27]:

Numerical approach
To understand to behavior of random gap fluctuations in graphene, and also the limitations of the SCBA, we carried out extensive numerical simulations on the honeycomb lattice, allowing for various random gap fluctuations on top of a uniform gap m. These fluctuations are simulated by box and Gaussian distributions. From the SCBA, the emergence of a second-order phase transition at a critical mean m c is obvious for a given variance. This is best manifested in the behavior of the DOS, which stays finite for m < m c , and vanishes afterwards, and serves as an order parameter. Does this picture indeed survive, when higher order corrections in the fluctuations are taken into account?
To start with, we take a fix random mass configuration with a given variance and the honeycomb lattice (HCL) with the conventional hoppings (t), represented by H 0 . Then, we take a separate Hamiltonian, responsible for the uniform, non-fluctuating gaps, denoted by H gap , and study the evolution of the eigenvalues of H 0 + mH gap by varying m for a 600x600 lattice. By using Lanczos diagonalization, we focus our attention only to the 200 eigenvalues closest to the NP. Their evolution is shown in Fig. 4. This supports the existence of a finite m c , but since it originates from a single random disorder configuration, rare events can alter the result. As a possible definition of the rigid gap, we also show the maximum of the energy level spacing for these eigenvalues as a function of m. As seen, it starts to increase abruptly at a given value of m, which can define m c .  To investigate whether a finite critical m c survives, we take smaller systems and evaluate the averaged DOS directly from many disorder realizations. To achieve this, we take a 200x200 HCL, and evaluate the 200 closest eigenvalues to the NP, and count their number in a given small interval, ∆E (smaller than the maximal eigenvalues) around zero. This method was found to be efficient in studying other types of randomness [28]. We mention that large values of ∆E take contribution from higher energy states into account, while too small values are sensitive to the discrete lattice and consequently the The resulting DOS is plotted in Figs. 5 and 6 for Gaussian (with variance g) and box distribution (within [−W..W ], variance g = W 2 /3). This does not indicate a sharp threshold, but rather the development of long Lifshitz tails due to randomness, as we already predicted in the previous section. To analyze them, we fitted the numerical data by assuming exponential tails of the form for a Gaussian and for a box distribution, as suggested by Ref. [29]. The obtained c values are visualized in the insets of Figs. 5 and 6. Given the good agreement, we believe that the DOS at the NP is made of states that are localized in a Lifshitz tail. We mention that these results are not sensitive to finite size scaling at these values of the disorder and uniform gap, only smaller systems (like the 30x30 HCL) require more averages (∼ 10 4 ), whereas for larger ones (such as the 200x200 with 400 averages) fewer averages are sufficient. In Fig. 7, the energy dependent DOS is shown for Gaussian random mass with g = 1 and for several uniform gap values. With increasing m, the DOS dimishes rapidly at low energies, and develops a pseudogap. The logarithmic singularity at E = t is washed out for g = 1. We also show the inverse of the DOS, proportional to T 0 , the characteristic temperature scale of variable range hopping as a function of the carrier density (which is proportional to E 2 ).

Discussion
MLG and BLG consist of two bands that touch each other at two nodal points (or valleys). Near the nodes the spectrum of MLG is linear (Dirac-like) and the spectrum of BLG is quadratic. The application of a uniform SBP opens a gap in the DOS for both cases. For a random SPB, however, the situation is less obvious. First of all, it is clear that randomness leads to a broadening of the bands. If we have two separate bands due to a small uniform SPB, randomness can close the gap again due to broadening (cf. Fig. 2a). The broadening of the bands depends on the strength of the fluctuations of the random SBP. In the case of a Gaussian distribution there are energy tails for all energies.
Now we focus on the NP, i.e. we consider E = 0 and ρ 0 . Then we have two parameters in order to change the gap structure: the average SBP m ≡m and the variance g.m allows us to broaden the gap and g has the effect of closing it due to broadening of the two subbands. Previous calculations have shown that at the critical value m c (g) of Eq. (29) the metallic behavior breaks down form > m c (g) [22]. On the other hand, Gaussian randomness creates tails at all energies. Consequently, there are localized states for |m| ≥ m c (g) at the NP, and there are delocalized states for |m| < m c (g) at the NP. The localized states in the tails are described, for instance, by the Lippmann-Schwinger equation (30) . The SCBA with uniform self-energy is not able to produce the localized tails. An extension of the SCBA with non-uniform self-energies provides localized tails though, as an approximation for largem has shown [27]. This is also in good agreement with our exact diagonalization of finite systems up to 200 × 200 size.
A possible interpretation of these results is that there are two different types of spectra. In a special realization of m r the tails of the DOS represent localized states. On the other hand, the DOS at the NP E = 0, obtained from the SCBA with uniform self-energy, comes from extended states [22]. The localized and the delocalized spectrum separate at the critical value m c , undergoing an Anderson transition.
Conductivity: Transport, i.e. the metallic regime, is related to the DOS trough the Einstein relation σ ∝ ρD, where D is the diffusion coefficient. The latter was found in Ref. [22] for E ∼ 0 as where a = 1 (a = 2) for MLG (BLG). Together with the DOS ρ 0 = η/2gπ and the scattering rate η in Eq. (28), the Einstein relation gives us at the NP In the localized regime (i.e. for |m| ≥ m c ) the conductivity is nonzero only for positive temperatures T > 0. Then we can apply the formula for variable-range hopping in Eq. (2), which fits well the experimental result in graphane of Ref. [7]. The parameter T 0 is related to the DOS at the Fermi level as [14] k where ξ is the localization length. T 0 has its maximum at the NP E F = 0, as shown in Fig. 7 and decreases monotonically with increasing carrier density, as in the experiment on graphane [7].
In conclusion, we have studied the density of states in MLG and BLG at low energies in the presence of a random symmetry-breaking potential. While a uniform symmetrybreaking potential opens a uniform gap, a random symmetry-breaking potential also creates tails in the density of states. The latter can close the gap again, preventing the system to become an insulator at the nodes. However, for a sufficiently large gap the tails contain localized states with nonzero density of states. These localized states allow the system to conduct at nonzero temperature via variable-range hopping. This result is in agreement with recent experimental observations [7].