High-energy Landau levels in graphene beyond nearest-neighbor hopping processes: Corrections to the effective Dirac Hamiltonian

We study the Landau level spectrum of bulk graphene monolayers beyond the Dirac Hamiltonian with linear dispersion. We consider an effective Wannier-like tight-binding model obtained from ab initio calculations, that includes long-range electronic hopping integral terms. We employ the Haydock-Heine-Kelly recursive method to numerically compute the Landau level spectrum of bulk graphene in the quantum Hall regime and demonstrate that this method is both accurate and computationally much faster than the standard numerical approaches used for this kind of study. The Landau level energies are also obtained analytically for an effective Hamiltonian that accounts for up to third nearest neighbor hopping processes. We find an excellent agreement between both approaches. We also study the effect of disorder on the electronic spectrum. Our analysis helps to elucidate the discrepancy between theory and experiment for the high-energy Landau levels energies.


I. INTRODUCTION
The pioneering experimental reports [1,2] on the unique features of the quantum Hall (QH) effect in graphene systems played a key role to establish the remarkable massless Dirac electronic properties of this material [3][4][5].Subsequent investigations on QH effects in graphene continued to produce fascinating results, such as Klein tunneling [6,7], fractional quantum Hall effect [8,9], Hofstadter butterflies [10,11], to name a few, attracting a lot of attention to the field (see, for instance, Ref. [12] for a recent review).
In addition to electronic transport properties, graphene under strong magnetic fields also shows unique spectral properties under strong magnetic fields.The graphene Landau levels have been theoretically predicted [4,13] to follow where N is the Landau level index, ω c = v F 2eB/ is the cyclotron frequency, B is the magnetic field strength, and v F stand for the electron Fermi velocity for B = 0.The LLs spectrum have been experimentally measured by transmission [14,15] and scanning tunneling spectroscopy [16,17] for both exfoliated and epithaxial graphene.These studies have verified the √ BN dispersion to a good approximation.A systematic transmission spectroscopy study [18] focused on the high-energy LL showed deviations from N predicted by Eq. (1).The puzzle is that the disagreement with the experimental data persists even when one improves the theoretical description by accounting next-to-nearest neighbor hopping processes.
In should be stressed that, while the massless Dirac Hamiltonian has been shown to be very effective in describing a variety of low energy electronic graphene properties, the theoretical modeling for higher energies is not unique.The standard approach uses density functional theory (DFT) calculations to obtain a tightbinding model based on Wannier orbitals, that contain an arbitrary range of hopping terms that fit the nonlinear features of the dispersion relation [19][20][21][22][23][24].
Another possibility for the discrepancy between theory and experiment is disorder, which is ubiquitous in graphene [25].Numerical investigations have studied the broadening the Landau subbands [26], but there is no systematic study of the corresponding disorder-induced peak shifts.We examine this issue with emphasis on the large |N | limit.
Our analysis uses the Haydock-Heine-Kelly (HHK) recursive method [27][28][29], also called the Haydock method, an order N [30] real-space computational approach developed to study local spectral functions.It transforms an arbitrary sparse Hamiltonian matrix in a tridiagonal form and evaluates the diagonal Green's function by a continued fraction expansion, avoiding the need of solving the full eigenvalue problem.The HHK method has been successfully used to compute the LDOS of different compounds [31,32] and more recently in the study of carbon nanotubes [33,34], and disordered graphene systems [35][36][37].
In addition to its efficiency, another attractive feature of the HHK method is that, since it relies on the concept of nearsightedness [38], it does not use periodic boundary conditions [39].Hence, it can be applied to study local spectral properties of disordered systems, quasicrystals [40,41], and systems with very large primitive unit cells, such as twisted graphene layers and 2D systems under realistic magnetic fields B. Surprising the HHK method has only been employed once in a case where B = 0, namely, in the investigation of Hoftstadter butterfly energy gaps in square lattices [42].To the best of our knowledge, so far no one has realized that the method is also fast and very accurate (as we show) for the study of discrete QH spectra.
Regarding our results, we numerically show that the large-|N | LL spectrum analytical solution of the contin-uum (long wavelength) effective graphene Hamiltonian is very accurate up to |N | ≤ 25 and B = 25 T. We include DFT-fitted third nearest hopping matrix elements into the graphene tight-binding model Hamiltonian and show that agreement between theory and experiment is significantly improved.
This paper is organized as follows.In Sec.II, we present the model Hamiltonian we use to describe the electronic properties in graphene, namely, a tight-binding Hamiltonian that accounts for up to third nearest neighbor hopping processes.Next, we discuss disorder effects and review the SCBA predictions for the shifts in the LL subband peak energies and widths.Finally, we briefly present the main steps of the implementation of the HHK method, discuss its computational cost and benchmark its accuracy.In Sec.III we present our results.We expand previous analytical results for N and show that, despite the approximations involved (discussed in App.A), the agreement with the numerical values is remarkable.We further show that the inclusion of third nearest neighbor matrix elements helps to improve the agreement between theory and experiment and that disorder plays a minor role.We summarize our conclusions in Sec.IV.

A. Model Hamiltonian
The tight-binding Hamiltonian that describes the electronic structure of graphene monolayers reads [3,20] where t ij stands for the hopping matrix element between the Wannier electronic orbitals centered at the carbon sites i and j.Most studies consider only first nearest neighbor hopping processes, a simple model that is able to describe the low energy properties of bulk graphene [3,25].Tight-binding parameterizations based on density functional theory (DFT) [19,21,22,24] show the necessity to including hopping terms beyond first nearest neighbors for a more accurate modeling of the electronic dispersion, particularly when addressing higher energies.
Here, we consider first t (1) , second t (2) , and third t (3)  nearest neighbor hopping terms.Within this approximation it is convenient to write the graphene Hamiltonian in a sublattice matrix representation that in reciprocal space reads [4,43] with and where a 1 = √ 3a 0 êx and a 2 = √ 3a 0 /2 êx + √ 3ê y are the primitive vectors honeycomb lattice and a 0 = 1.4A is the carbon-carbon bond length [3].The first and third nearest neighbor hopping terms, that connect different sublattices [3,4], correspond to the off-diagonal matrix elements, while the second-nearest hoppings are related to the diagonal ones.This apparent correspondence between even-odd nearest neighbors and diagonal off-diagonal matrix elements breaks down for fourth nearest neighbors and beyond [23].In Eq. ( 3) we neglect a constant diagonal term that shifts the energy spectrum by −3t (2) .
The energy dispersion reads where λ labels the valence (λ = −1) and conduction bands (λ = +1).Note that the second-neighbors hopping contributions do not depend of λ and, thus, break the electron-hole symmetry.
The presence of an external magnetic field B can be accounted by the Peierls substitution [44,45], that is, by the transformation where R n is the lattice vector associated with the site n, e is the electron charge.Here, we choose the vector potential A(r) = (0, Bx, 0) that gives a magnetic field B = ∇ × A = Bê z perpendicular to the graphene layer.

B. Numerical method
We compute the Landau level spectra using the Haydock-Heine-Kelly (HHK) recursion technique [27][28][29].The latter has been developed to calculate the local properties of electronic systems represented on a basis of localized (orthogonal) states |i , like the one used in the tight-binding Hamiltonian of Eq. ( 2).The HHK method provides a very efficient O(N ) recursive procedure to transform a given Hamiltonian matrix into a tridiagonal one, that is much more amenable for numerical calculation.As mentioned introduction, the HHK method has been successfully used to compute the LDOS of continuous spectra of several systems [31][32][33][34][35][36][37].Here, we use it to compute the LDOS of a Hamiltonian with discrete eigenvalues.
Let us quickly review the main ingredients of the HHK method before discussing its use in computing the graphene Landau level spectrum.
The recursion method starts by targeting a given state |0} = |j .The method generates a hierarchy of states |n} based on the three-term recursion relations [27][28][29] with the recursive coefficients and the orthogonal basis element where b 0 = 0 and | − 1} = 0. Thus, by construction, the Hamiltonian matrix in the orthogonal basis {|n}} is tridiagonal.In turn, the basis functions |n} can be expressed in terms of Wannier-like states |j , Here, we assume that the electronic wave functions are a superposition of P states centered at the atomic sites i and the atomic orbitals are orthogonal to each other, in line with the tight binding model of Eq. ( 2).
The diagonal Green's function for the seed state |0} is given by continued fraction [27][28][29] expressed in terms of the matrix elements of the tridiagonal Hamiltonian in the basis |n}.The LDOS at any site j can be written as In practice, a finite η serves as a convenient regularization parameter.For continuous spectra, setting η ≈ 2D/M , where D is the bandwidth, guarantees a nice smooth approximation to LDOS(j, ) [27].
For pristine systems, due to translational symmetry, the LDOS(j, ) at any j is proportional to the (total) density of states ρ( ) [46].In our calculations we fix j at the center of the honeycomb lattice of size P .
In previous applications [31,32,42] it has been observed that for a sufficiently large M , the recursive coefficients converge towards their asymptotic values, namely, a M → a ∞ and b M → b ∞ .The asymptotic value of a ∞ is associated with the center of the energy band 0 , namely, 0 = a ∞ /2 [27,28].In graphene systems, a ∞ /2 corresponds to the Dirac point energy.For the nearestneighbor tight-binding Hamiltonian, the center band energy is 0 = 0 and all a n coefficients are zero.In this case, by inspecting Eq. ( 13), one immediately finds that G 00 ( ) = −G 00 (− ).This implies that the LDOS(j, ) has electron-hole symmetry for any j and ρ( 0 ) = 0, a condition that defines the so-called Dirac points [3].When t (2) = 0, the coefficient a n are no longer zero, the Dirac points are energy shifted, and the electronhole symmetry broken.These simple properties are in line with well established literature results [3], as they should.
Let us now apply the HHK method to graphene systems.We begin showing results of the LDOS for bulk graphene in the absence of external magnetic fields.Throughout this paper we use the tight-binding parameters given in Table I.
The optimal number of iterations M depends on P , the size of the system chosen to represent the bulk, as well as on the desired accuracy.A detailed analysis about the dependence of the HHK method accuracy on M for continuum spectra can be found in Ref. [28].For graphene in the absence of magnetic field, we find that M ≈ √ P guarantees a good accuracy.At the end of this section we study the accuracy for the case where B = 0.

Parameterization
t (1) (eV) -3.0933 0.19915 -0.16214 Figure 1 shows the LDOS of a graphene monolayer obtained from the HHK method for P = 2.5 × 10 7 carbon atoms and M = 5000 iterations.We have contrasted the the latter with ρ( ) computed by a direct numerical evaluation of ρ( ) = 1 ABZ k δ( − kλ ), where kλ is given by Eq. ( 6) and the wavevectors k are sampled over the Brillouin zone of area A BZ .By accounting for the normalization factor that relates LDOS(j, ) with ρ( ), we observe a nice agreement within the numerical precision imposed by the regularization parameter η.We find that differences are only appreciable, as expected, at the band edges as well as at the van Hove singularities.
Let us now consider the case of a pristine graphene sheet under an external perpendicular magnetic field B. The electronic spectrum becomes discrete and strongly degenerate forming a sequence of Landau levels.Let us now benchmark the HHK results against another well-established method to compute the LL energies [47].The most straightforward numerical solution solves the eigenvalue problem for supercells with periodic boundary conditions whose size are dictated by the magnetic field strength and the phases of the Peierls substitution [48].Here, the computation of the LL spectra involves the diagonalization of matrices of dimension ∼ 9.07 × 10 4 /B[T] [49], that is computationally very costly for realistic values of the magnetic field.Alternatively, Ref. [47] proposed to consider systems with an infinite ribbon geometry where, by a suitable gauge choice, the magnetic field is encoded in the hopping terms transversal to the ribbon preserving a "free" dynamics in the longitudinal direction.In this way, the primitive unit cell has N W sites that depend linearly on the ribbon width W .The k-space domains where the bands are flat, corresponding to energies GNR N , are associated with bulk Landau levels, namely, N ≈ GNR N .The accuracy relies on W/ B 1, where B = /eB ≈ 26 B[T] nm.We consider graphene nanoribbons with zigzag edges for which W = (3N W /4 − 1)a 0 .Figure 3 shows the convergence of N as a function of N W .In this approach, an accurate assessment of GNR N , with negligible inter edge hybridization, requires increasing W (or N W ) for larger |N |, with a computational cost that scales with N 3 W .We note that by finding an optimized truncated basis set with N opt < N W one can significantly reduce computational time, but the number of operations of such procedure still scales as N 3 opt .Thus, by comparing Figs.2(b) and 3 and recalling that the HHK as an O(M ) method, one concludes that the HHK method is rather unexpensive.

C. Disorder
Several studies have addressed the effects of disorder in graphene in the quantum Hall regime [26,[50][51][52][53]. Here, our goal is to obtain insight on disorder effects in the density of states (DOS) of large |N | Landau subbands, particularly on their energy peaks, rather than studying specifics of a given experimental setup.For this reason we consider a simple model that captures the main features of disorder in graphene, namely, the Anderson disorder model.The on-site energy u i is assumed to be uniformly distributed between −U and +U , that is, with zero mean u i = 0 and variance u 2 i = U 2 /3.This model describes the short range disorder regime [50,52,53] where the impurities scattering range is much smaller than the lattice constant and, thus, impurity scattering mixes the ξ = 1 (or K) and the ξ = −1 (or K ) valley spinors.
We compare our numerical simulation with analytical results.In a pioneering work, Shon and Ando [50] used the self-consistent Born approximation (SCBA) to calculate self-energy impurity average Σ in graphene for different kinds of disorder.For short range impurities [50], where is the parameter characterizing the disorder scattering strength and n imp denotes the impurity concentration.The cutoff N c is introduced, to account for the finite bandwidth.
Contrary to the constant spacing of LL in twodimensional electron gas systems, in graphene the spacing decreases with ω c / |N |.This feature makes the evaluation of Σ( ) somewhat more involved than in the standard case [53].In the limit of weak scattering strength γ 2  1 and considering the usual inequalities | − N |, |Σ| N and | − Σ| < ω c / |N | that correspond to isolated LL subbands, the density of states is given by [53] DOS where and D = ω c √ N c is a Debye-like energy cutoff.The DOS is a sequence of semi-circles (or better, elliptical shapes) with maxima at N (1 − χ) and full width at half maximum Γ N = γ 3 (1 − χ).Hence, disorder induces a ∆ N = −χ N shift in the LL subband peaks.Note that by taking χ = 0, one recovers the simple single-isolated subband approximation found in Ref. [50].

III. RESULTS
In this section we use the HHK method to study the LL energies N , as a function of |N |, the external magnetic field B, and disorder.We also compare the results with the experimental data [18].

A. Landau levels in pristine graphene
The analysis of the LL spectra measured in Ref. [18] reports a discrepancy between the analytical approximation for N and the corresponding experimental values.The problem has been attributed to the model Hamiltonian truncation to second nearest neighbor hopping terms.Here, we expand the original theoretical analysis by accounting for third nearest neighbor terms and compare our results both with numerical calculations and with the experimental data.
The dispersion relation of graphene can be expressed as a power series in q, defined as q = k ∓ K, where ±K correspond to γ ±K = γ ±K = 0 that defines the so-called Dirac points.The external magnetic field is accounted for by minimal substitution, Π = q + eA.
By using canonical quantization and the approximations suggested in Ref. [4] (see App.A for details), we write the energies of the LLs in the large |N | limit up to third-nearest neighbors hopping as where the first two terms (1) (2) were obtained in Refs.[4,18].Our derivation expands the latter results to account for t (3) , namely In what follows we show that Eqs. ( 20) to (22) give an excellent approximation to the LLs energies obtained from the tight-binding Hamiltonian of Eqs. ( 2) and (7).
Figure 4(a) compares the analytical N with the numerical HHK results to q 3 that have been neglected in the momentum expansion, see App. A.
Figure 4(b) contrasts the analytical N second nearest neighbor contribution to the LL spectrum with its numerical counterpart (2) N .The later is defined as N , where N and (1) N are computed numerically.Both are obtained using the tight-binding parameter set B, but for the calculation of  (2) is small discrepancy between the constant diagonal term that shifts the energy spectrum by −3t (2)  (t (2) = 0.3 eV, see Table I) in the analytical approach and the HHK results where the shift is −0.899 eV.
Finally, Fig. 4 (c) compares N with (3) .Here, we define N , where we use the HHK method to compute N with the parameter set C and Let us now compare our results with the experimental data [18].Figure 5 (a) shows that, since to leading order N ∼ |N | 1/2 , all tight-binding parameterizations correctly capture the main trend of N .In Fig. 5  for the conduction and the valence bands.We find that the parameterization B (up to second nearest neighbor hopping), considered in Ref. [18] shows the most agreement with the experimental data, with . This explains the effectiveness on the extra arbitrary "W" parameter used in Ref. [4,18] to obtain a good fit.The parameterization C (up to third nearest neighbor hopping) gives the best results, but still shows some deviations from the experiments.We stress that the tight-binding parameters are fits of DFT calculations to k for B = 0 and not to the LL energies N .

B. Disordered graphene
Let us now investigate the effect of disorder in the LL spectrum of bulk graphene.We briefly address the disorder induced LL broadening Γ N but our main interest is the energy shift ∆ N .
Disorder breaks translational invariance and the sites are no longer equivalent.By involving the ergodic hypothesis one identifies the ensemble average of the LDOS at a given system site with the DOS of the disordered system.In this study, the ensemble averages involve 10 4 disorder realizations.We use η = 0.5 meV as regularization parameter.This choice guarantees that η Γ N for the disorder strengths we study and, thus, has a negligible effect on the results.
Since the short-range disorder effects on the DOS are not expected to depend on the tight-binding parameters, here we only consider the parameterization A. The SCBA equation, Eq. ( 15), can be numerically solved by a simple iteration starting with a given initial value of the selfenergy.We have introduced a Debye-like cutoff N c = 1342 and start the self-consistent loop with an initial selfenergy of Im Σ( ) = 1 meV.The iteration converges very rapidly.
Figure 6 compares the HHK with the SCBA results for the DOS( ) of graphene at B = 25 T for different disorder strengths, namely, u 2 i ≈ 0.02 and 0.08.The Landau subbands are broadened and shifted toward the zeroenergy in the presence of on-site disorder with the exception of the zero-energy subband which is only broadened,  I) and B = 25 T. The experimental data are obtained from Ref. [18].Deviation between the computed and the experimental LL energies, HHK Figures 6(c) to 6(e) compare the HHK with the SCBA density of states of different LL subbands.The SCBA predicts a semi-circular (or, more precisely, elliptical) shaped DOS( ) for all N .For IQH systems with quadratic dispersion, it is well established that the DOS of the N = 1 subband has Gaussian-like shape [54] and the semi-circular DOS is expected for N 1.Our results show large deviations from the SCBA independent of N .This observation is consistent with previous numerical investigations of the DOS in graphene at the QH regime [26] and deserve further investigation.
Figure 7 shows the disorder renormalization of the Landau level subbands peak energies ∆ N .We define ∆ HHK N as the difference between the ensemble average HHK N and the pristine value HHK N .The ∆ SCBA N can be obtained from the numerical solution of Eq. ( 15) or evaluated from the approximation given by Eq. ( 17).The agreement is very good.
The results show that disorder contributes to increase  the deviation between theory and experiment.However, for realistic disorder strengths, ∆ N is very small.The black and red symbols correspond to the numerical results of HHK and SCBA, respectively, and the blue symbol corresponds to the analytical prediction of Ref. [53].

IV. CONCLUSIONS AND DISCUSSION
In this paper we have studied the large |N | Landau level spectrum of graphene monolayers beyond the firstnearest neighbor tight-binding approximation and the effective low-energy Dirac Hamiltonian.The study has to goals, namely, to discuss the effect of different kinds of hopping processes to the high-energy LLs in graphene and to introduce an efficient numerical method for this analysis.
Regarding the methodology, our results show that the HHK method [27][28][29] is very efficient for the computation of the LDOS of graphene in the QH regime.This method has a long record of success, but to the best of our knowledge has not been used to study discrete spectra.Our paper shows that the HHK method is very accurate and, since it is an O(N ) approach, it is, computationally much faster then other methods found in the literature [47,49].
As for the large-|N | LL spectrum, we show that the analytical solution of the continuum (long wavelength) effective graphene Hamiltonian is very accurate up to |N | ≤ 25 and B = 25 T. We find that by including third nearest hopping processes into the tight-binding model Hamiltonian the agreement between theory and experiment is significantly improved.Although the tightbinding parameters can vary depending on the methodology employed to extract them from DFT calculations, the variations are small and we expect our results to be robust.
We also analyzed the effect of disorder.We find a very good agreement between numerical simulations and the SCBA for the disorder induced renormalization of the LL subbands energy peaks ∆ N .For realistic disorder strengths, ∆ N gives a small correction to N and, thus, does not impact our conclusions.Interestingly, our simulations indicate that even for large-|N | the subband DOS does not show a semi-circular shape, at odds with the theoretical studies of 2D electron gas systems in the QH regime [55][56][57][58].We believe this issue deserves further investigation.
We expect the HHK method to be usefull for the analysis of the LDOS of other QH 2D systems and, since it does not rely on periodic boundary conditions, it can also used in the study of quasicrystals [40,41] and fractal lattices [59].
Figure 2(a) shows the graphene LDOS for the tight-binding parameterizations A, B, and C computed with the HHK method for B = 25 T, P = 2.25 × 10 6 , M = 1500, and η = 0.1 meV.Due to the finite η, the LL are broadened and become Lorentzian distributions, whose energy peaks (obtained by fitting) are associated with the LL energies N .Here, we chose η | N − N −1 |, guaranteeing that the |N | ≤ 30 lowest LL peaks are nicely resolved.

LFIG. 1 .
FIG.1.Graphene LDOS (in arbitrary units) as a function of the energy (in eV) calculated with the HHK method for the tight-binding Hamiltonian parameterizations A, B, and C. We set 0 = 0 for better visualization.

Figure 2 (
Figure 2(b) shows the convergence of the Landau level energies N as a function of M for a lattice size of P = 2.25 × 10 6 , B = 25 T, η = 0.1 meV and the tight-binding parameterization A. Throughout this paper we set the absolute numerical precision to δ = 10 −8 eV.We find that by setting M * = 1500, the accuracy of N is better than δ for |N | ≤ 30.Figure 2(b) shows that the LL energies ε N converge rapidly to their asymptotic values as the number of iterations M increases.The results indicate that the number of iterations M necessary to obtain a given precision δ scales as M ∼ |N |.Let us now benchmark the HHK results against another well-established method to compute the LL energies[47].The most straightforward numerical solution solves the eigenvalue problem for supercells with periodic boundary conditions whose size are dictated by the magnetic field strength and the phases of the Peierls substitution[48].Here, the computation of the LL spectra involves the diagonalization of matrices of dimension ∼ 9.07 × 10 4 /B[T][49], that is computationally very costly for realistic values of the magnetic field.Alternatively, Ref.[47] proposed to consider systems with an infinite ribbon geometry where, by a suitable gauge choice, the magnetic field is encoded in the hopping terms transversal to the ribbon preserving a "free" dynamics in the longitudinal direction.In this way, the primitive unit cell has N W sites that depend linearly on the ribbon width W .The k-space domains where the bands are flat, corresponding to energies GNR

MFIG. 2 .
FIG. 2. (a) LDOS (in arbitrary units) versus (in eV) for B = 25 T, η = 0.1 meV and M * = 1500 for the tight-binding parameterizations A (red), B (blue) and C (green line).(b) Accuracy of the LL energies ∆ N = | M N − M * N | as a function of the number of iterations M for the tight-binding parameterization A. Here, M * N stands for the energy of the N th Landau level converged to a precision better than δ = 10 −8 eV.
) Note the subleading terms in the above expansions give minor corrections, since a 0 / B = 5.4 × 10 −3 B[T] is small for B fields within the experimental reach.

FIG. 4 .
FIG. 4. Comparison between the analytical contributions toN and HHK results for a pristine graphene monolayer with B = 25 T. (a) first neighbor contribution, (b) second neighbor contribution and (c) third neighbor contribution.The insets show the deviation between the HHK results and the analytical approximation, |∆ N / HHK N |, see text for details.

2 EFIG. 5 .
FIG. 5. (a) The Landau levels spectrum as a function of |N | for the tight-binding parameterizations A, B and C (see TableI) and B = 25 T. The experimental data are obtained from Ref.[18].Deviation between the computed and the experimental LL energies, HHK the (b) conduction and (c) valence bands.as nicely shown by Figs.6(a) and 6(b).

FIG. 7 .
FIG.7.Renormalization of the LL subband peak energy ∆N as function of N for u 2 i ≈ 0.02, 0.08 and 0.19.The black and red symbols correspond to the numerical results of HHK and SCBA, respectively, and the blue symbol corresponds to the analytical prediction of Ref.[53].

TABLE I .
Tight-binding hopping parameter sets in eV.