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

Kevin J. U. Vidarte  Instituto de Física, Universidade Federal do Rio de
Janeiro, 21941-972 Rio de Janeiro - RJ, Brazil  Caio Lewenkopf  Instituto de
Física, Universidade Federal Fluminense, 24210-346 Niterói - RJ, Brazil  Max
Planck Institute for the Physics of Complex Systems, 01187 Dresden, Germany

###### Abstract

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 [13, 4] to follow

$\epsilon_{N}={\rm sgn}(N)\hbar\omega_{c}\sqrt{|N|}$ (1)

where $N$ is the Landau level index, $\omega_{c}=v_{F}\sqrt{2eB/\hbar}$ 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 $\sqrt{BN}$ dispersion to
a good approximation. A systematic transmission spectroscopy study [18]
focused on the high-energy LL showed deviations from $\epsilon_{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 tight-binding 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 ${\cal 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\neq 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 continuum (long wavelength) effective graphene
Hamiltonian is very accurate up to $|N|\leq 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
$\epsilon_{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.

##  II Theory and methods

###  II.1 Model Hamiltonian

The tight-binding Hamiltonian that describes the electronic structure of
graphene monolayers reads [20, 3]

$H=\sum_{i,j}\left(t_{ij}c^{\dagger}_{i}c_{j}+{\rm H.c}\right),$ (2)

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, 22,
24, 21] 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 [43, 4]

$H_{\textbf{k}}\equiv\left(\begin{array}[]{cc}t^{(2)}|\gamma_{\textbf{k}}|^{2}&t^{(1)}\gamma_{\textbf{k}}^{*}+t^{(3)}\gamma_{\textbf{k}}^{\prime}\\\
t^{(1)}\gamma_{\textbf{k}}+t^{(3)}\gamma_{\textbf{k}}^{\prime*}&t^{(2)}|\gamma_{\textbf{k}}|^{2}\end{array}\right),$
(3)

with

$\gamma_{\textbf{k}}\equiv
1+e^{i\textbf{k}\cdot\textbf{a}_{2}}+e^{i\textbf{k}\cdot(\textbf{a}_{2}-\textbf{a}_{1})}$
(4)

and

$\gamma^{\prime}_{\textbf{k}}\equiv
1+e^{i2\textbf{k}\cdot\textbf{a}_{2}}+e^{i2\textbf{k}\cdot(\textbf{a}_{2}-\textbf{a}_{1})},$
(5)

where $\textbf{a}_{1}=\sqrt{3}a_{0}\hat{\bf e}_{x}$ and
$\textbf{a}_{2}=\sqrt{3}a_{0}/2\left(\hat{\bf e}_{x}+\sqrt{3}\hat{\bf
e}_{y}\right)$ are the primitive vectors honeycomb lattice and
$a_{0}=1.4\,$\mathrm{\SIUnitSymbolAngstrom}$$ 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

$\epsilon_{{\bf
k}\lambda}=t^{(2)}|\gamma_{\textbf{k}}|^{2}+\lambda|t^{(1)}\gamma_{\textbf{k}}+t^{(3)}\gamma_{\textbf{k}}^{\prime*}|,$
(6)

where $\lambda$ labels the valence ($\lambda=-1$) and conduction bands
($\lambda=+1$). Note that the second-neighbors hopping contributions do not
depend of $\lambda$ 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

$t_{ij}\longrightarrow
t_{ij}\,\exp\left[i\dfrac{e}{\hslash}\int_{\textbf{R}_{i}}^{\textbf{R}_{j}}d\textbf{r}\cdot\textbf{A}(\textbf{r})\right],$
(7)

where $\textbf{R}_{n}$ is the lattice vector associated with the site $n$, $e$
is the electron charge. Here, we choose the vector potential
$\textbf{A}(\textbf{r})=\left(0,Bx,0\right)$ that gives a magnetic field
$\textbf{B}=\nabla\times\textbf{A}=B\hat{\bf e}_{z}$ perpendicular to the
graphene layer.

###  II.2 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\rangle}$, like the one used in the tight-binding
Hamiltonian of Eq. (2). The HHK method provides a very efficient $O({\cal 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 [32, 31, 35, 36, 33, 34, 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\rangle$. The
method generates a hierarchy of states $|n\\}$ based on the three-term
recursion relations [27, 28, 29]

$\displaystyle H|n\\}=a_{n}|n\\}+b_{n}|n-1\\}+b_{n+1}|n+1\\},$ (8)

with the recursive coefficients

$\displaystyle a_{n}=\\{n|H|n\\},$ (9) $\displaystyle
b_{n+1}=\|(H-a_{n})|n\\}-b_{n}|n-1\\}\|,$ (10)

and the orthogonal basis element

$\displaystyle|n+1\\}=\frac{1}{b_{n+1}}\left[(H-a_{n})|n\\}-b_{n}|n-1\\}\right],$
(11)

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\rangle$,

$\displaystyle|n\\}=\sum_{i=1}^{P}A_{ni}|i\rangle.$ (12)

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]

$\displaystyle\begin{split}G_{00}(\epsilon)&=\\{0|\dfrac{1}{\epsilon-H}|0\\}\\\
&=\dfrac{1}{\epsilon-a_{0}-{\textstyle\dfrac{b_{1}^{2}}{\epsilon-
a_{1}-\dfrac{b_{2}^{2}}{\epsilon-
a_{2}-\dfrac{b_{3}^{2}}{\ddots}}}}},\end{split}$ (13)

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

$\displaystyle{\rm LDOS}(j,\epsilon)$ $\displaystyle=-\frac{1}{\pi}{\rm
Im}\;G^{r}_{jj}(\epsilon)$ (14)
$\displaystyle\equiv-\frac{1}{\pi}\lim_{\eta\rightarrow 0^{+}}\left[{\rm
Im}\;G_{jj}(\epsilon+i\eta)\right].$

In practice, a finite $\eta$ serves as a convenient regularization parameter.
For continuous spectra, setting $\eta\approx 2D/M$, where $D$ is the
bandwidth, guarantees a nice smooth approximation to LDOS$(j,\epsilon)$ [27].

For pristine systems, due to translational symmetry, the LDOS$(j,\epsilon)$ at
any $j$ is proportional to the (total) density of states $\rho(\epsilon)$
111Here, LDOS$(j,\epsilon)$ is normalized to unit and $\rho(\epsilon)$ to
$2{\cal A}_{BZ}$ where ${\cal A}_{BZ}$ is the Brillouin zone area.. In our
calculations we fix $j$ at the center of the honeycomb lattice of size $P$.

In previous applications [42, 31, 32] it has been observed that for a
sufficiently large $M$, the recursive coefficients converge towards their
asymptotic values, namely, $a_{M}\rightarrow a_{\infty}$ and $b_{M}\rightarrow
b_{\infty}$. The asymptotic value of $a_{\infty}$ is associated with the
center of the energy band $\epsilon_{0}$, namely, $\epsilon_{0}=a_{\infty}/2$
[27, 28]. In graphene systems, $a_{\infty}/2$ corresponds to the Dirac point
energy. For the nearest-neighbor tight-binding Hamiltonian, the center band
energy is $\epsilon_{0}=0$ and all $a_{n}$ coefficients are zero. In this
case, by inspecting Eq. (13), one immediately finds that
$G_{00}(\epsilon)=-G_{00}(-\epsilon)$. This implies that the
LDOS$(j,\epsilon)$ has electron-hole symmetry for any $j$ and
$\rho(\epsilon_{0})=0$, a condition that defines the so-called Dirac points
[3]. When $t^{(2)}\neq 0$, the coefficient $a_{n}$ are no longer zero, the
Dirac points are energy shifted, and the electron-hole 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 1.

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\approx\sqrt{P}$ guarantees a good accuracy. At the end
of this section we study the accuracy for the case where $B\neq 0$.

Table 1: Tight-binding hopping parameter sets in eV. Parameterization |  $t^{(1)}$ (eV) |  $t^{(2)}$ (eV) |  $t^{(3)}$ (eV)
---|---|---|---
A [3] | -2.7 |  |
B [19] | -3.0 | 0.3 |
C [22] | -3.0933 | 0.19915 | -0.16214

Figure 1 shows the LDOS of a graphene monolayer obtained from the HHK method
for $P=2.5\times 10^{7}$ carbon atoms and $M=5000$ iterations. We have
contrasted the the latter with $\rho(\epsilon)$ computed by a direct numerical
evaluation of $\rho(\epsilon)=\frac{1}{{\cal A}_{\rm
BZ}}\sum_{\textbf{k}}\delta(\epsilon-\epsilon_{\textbf{k}\lambda})$, where
$\epsilon_{\textbf{k}\lambda}$ is given by Eq. (6) and the wavevectors ${\bf
k}$ are sampled over the Brillouin zone of area ${\cal A}_{\rm BZ}$. By
accounting for the normalization factor that relates LDOS$(j,\epsilon)$ with
$\rho(\epsilon)$, we observe a nice agreement within the numerical precision
imposed by the regularization parameter $\eta$. We find that differences are
only appreciable, as expected, at the band edges as well as at the van Hove
singularities.

Figure 1: Graphene LDOS (in arbitrary units) as a function of the energy
$\epsilon$ (in eV) calculated with the HHK method for the tight-binding
Hamiltonian parameterizations A, B, and C. We set $\epsilon_{0}=0$ for better
visualization.

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. 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\times 10^{6}$, $M=1500$, and
$\eta=0.1$ meV. Due to the finite $\eta$, the LL are broadened and become
Lorentzian distributions, whose energy peaks (obtained by fitting) are
associated with the LL energies $\epsilon_{N}$. Here, we chose
$\eta\ll|\epsilon_{N}-\epsilon_{N-1}|$, guaranteeing that the $|N|\leq 30$
lowest LL peaks are nicely resolved.

Figure 2(b) shows the convergence of the Landau level energies $\epsilon_{N}$
as a function of $M$ for a lattice size of $P=2.25\times 10^{6}$, $B=25$ T,
$\eta=0.1$ meV and the tight-binding parameterization A. Throughout this paper
we set the absolute numerical precision to $\delta=10^{-8}$ eV. We find that
by setting $M^{*}=1500$, the accuracy of $\epsilon_{N}$ is better than
$\delta$ for $|N|\leq 30$. Figure 2(b) shows that the LL energies
$\varepsilon_{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 $\delta$ scales as
$M\sim\sqrt{|N|}$.

Figure 2: (a) LDOS (in arbitrary units) versus $\epsilon$ (in eV) for $B=25$
T, $\eta=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
$\Delta\epsilon_{N}=|\epsilon_{N}^{M}-\epsilon_{N}^{M^{*}}|$ as a function of
the number of iterations $M$ for the tight-binding parameterization A. Here,
$\epsilon_{N}^{M^{*}}$ stands for the energy of the $N$th Landau level
converged to a precision better than $\delta=10^{-8}$ eV.

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 $\sim 9.07\times 10^{4}/B[{\rm
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 $\epsilon_{N}^{\rm
GNR}$, are associated with bulk Landau levels, namely,
$\epsilon_{N}\approx\epsilon_{N}^{\rm GNR}$. The accuracy relies on
$W/\ell_{B}\gg 1$, where $\ell_{B}=\sqrt{\hbar/eB}\approx 26\sqrt{B[{\rm T}]}$
nm.

We consider graphene nanoribbons with zigzag edges for which
$W=(3N_{W}/4-1)a_{0}$. Figure 3 shows the convergence of $\epsilon_{N}$ as a
function of $N_{W}$. In this approach, an accurate assessment of
$\epsilon_{N}^{\rm GNR}$, with negligible inter edge hybridization, requires
increasing $W$ (or $N_{W}$) for larger $|N|$, with a computational cost that
scales with $N_{W}^{3}$. We note that by finding an optimized truncated basis
set with $N_{\rm opt}<N_{W}$ one can significantly reduce computational time,
but the number of operations of such procedure still scales as $N_{\rm
opt}^{3}$. Thus, by comparing Figs. 2(b) and 3 and recalling that the HHK as
an ${\cal O}(M)$ method, one concludes that the HHK method is rather
unexpensive.

Figure 3: Accuracy of the LLs energies $\Delta\epsilon_{N}=|\epsilon_{N}^{\rm
GRN}-\epsilon_{N}^{M^{*}}|$, for graphene nanoribbons with zigzag edges as a
function of their width $W=(3N_{W}/4-1)a_{0}$.

###  II.3 Disorder

Several studies have addressed the effects of disorder in graphene in the
quantum Hall regime [50, 51, 52, 53, 26]. 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 $\left\langle
u_{i}\right\rangle=0$ and variance $\langle\langle
u_{i}^{2}\rangle\rangle=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
$\xi=1$ (or $K$) and the $\xi=-1$ (or $K^{\prime}$) 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 $\Sigma$ in graphene for different
kinds of disorder. For short range impurities [50],

$\displaystyle\Sigma(\epsilon)=\dfrac{\gamma^{2}}{4}\sum^{N_{c}}_{N=-N_{c}}\dfrac{1}{\epsilon-\epsilon_{N}-\Sigma(\epsilon)},$
(15)

where

$\displaystyle\gamma^{2}=\dfrac{\langle\langle u_{i}^{2}\rangle\rangle n_{\rm
imp}}{\pi\ell_{B}^{2}}$ (16)

is the parameter characterizing the disorder scattering strength and $n_{\rm
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 two-dimensional electron gas
systems, in graphene the spacing decreases with
$\hslash\omega_{c}/\sqrt{|N|}$. This feature makes the evaluation of
$\Sigma(\epsilon)$ somewhat more involved than in the standard case [53]. In
the limit of weak scattering strength $\gamma^{2}\ll 1$ and considering the
usual inequalities $|\epsilon-\epsilon_{N}|$, $|\Sigma|\ll\epsilon_{N}$ and
$|\epsilon-\Sigma|<\hslash\omega_{c}/\sqrt{|N|}$ that correspond to isolated
LL subbands, the density of states is given by [53]

$\displaystyle{\rm
DOS}(\epsilon)=\dfrac{2}{\pi\gamma^{2}}\sqrt{\gamma^{2}\left(1-\chi\right)-\left[\epsilon-\epsilon_{N}\left(1-\chi\right)\right]^{2}},$
(17)

where

$\displaystyle\chi=\dfrac{\gamma^{2}}{\left(\hslash\omega_{c}\right)^{2}}\ln{\left(D/\epsilon_{N}\right)},$
(18)

and $D=\hbar\omega_{c}\sqrt{N_{c}}$ is a Debye-like energy cutoff. The DOS is
a sequence of semi-circles (or better, elliptical shapes) with maxima at
$\epsilon_{N}\left(1-\chi\right)$ and full width at half maximum
$\Gamma_{N}=\gamma\sqrt{3\left(1-\chi\right)}$. Hence, disorder induces a
$\Delta_{N}=-\chi\epsilon_{N}$ shift in the LL subband peaks. Note that by
taking $\chi=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 $\epsilon_{N}$,
as a function of $|N|$, the external magnetic field $B$, and disorder. We also
compare the results with the experimental data [18].

###  III.1 Landau levels in pristine graphene

The analysis of the LL spectra measured in Ref. [18] reports a discrepancy
between the analytical approximation for $\epsilon_{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 $\bf
q$, defined as ${\bf q}={\bf k}\mp{\bf K}$, where $\pm{\bf K}$ correspond to
$\gamma_{\pm\textbf{K}}=\gamma_{\pm\textbf{K}}^{\prime}=0$ that defines the
so-called Dirac points. The external magnetic field is accounted for by
minimal substitution, ${\bm{\Pi}}=\hbar{\bf q}+e{\bf A}$.

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

$\widetilde{\epsilon}_{N}=\widetilde{\epsilon}_{N}^{(1)}+\widetilde{\epsilon}_{N}^{(2)}+\widetilde{\epsilon}_{N}^{(3)},$
(19)

where the first two terms

$\displaystyle\widetilde{\epsilon}_{N}^{(1)}$ $\displaystyle\approx$
$\displaystyle{\rm
sgn}(N)\hslash\omega_{c}|N|^{1/2}\left(1-\dfrac{3}{8}\dfrac{a_{0}^{2}}{\ell_{B}^{2}}|N|\right),$
(20) $\displaystyle\widetilde{\epsilon}_{N}^{(2)}$ $\displaystyle\approx$
$\displaystyle\hslash\omega_{c}\dfrac{t^{(2)}}{|t^{(1)}|}\dfrac{3}{\sqrt{2}}\dfrac{a_{0}}{\ell_{B}}|N|\left(1-\dfrac{3}{4}\dfrac{a_{0}^{2}}{\ell_{B}^{2}}|N|\right),$
(21)

were obtained in Refs. [18, 4]. Our derivation expands the latter results to
account for $t^{(3)}$, namely

$\widetilde{\epsilon}_{N}^{(3)}\approx-{\rm
sgn}(N)\hslash\omega_{c}\dfrac{2t^{(3)}}{t^{(1)}}|N|^{1/2}\left(1-\dfrac{t^{(3)}}{t^{(1)}}-\dfrac{59}{32}\dfrac{a_{0}^{2}}{\ell_{B}^{2}}|N|\right).$
(22)

Note the subleading terms in the above expansions give minor corrections,
since $a_{0}/\ell_{B}=5.4\times 10^{-3}\sqrt{B[{\rm T}]}$ is small for $B$
fields within the experimental reach.

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: Comparison between the analytical contributions to
$\widetilde{\epsilon}_{N}$ 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,
$|\Delta\epsilon_{N}/\epsilon_{N}^{{\rm HHK}}|$, see text for details.

Figure 4(a) compares the analytical $\widetilde{\epsilon}_{N}^{(1)}$ with the
numerical HHK results $\epsilon_{N}^{(1)}=\epsilon_{N}^{\rm HHK}$ for the
first nearest neighbor tight-binding Hamiltonian, namely, the parameterization
A. The insert shows the relative deviation
$|\Delta\epsilon_{N}/\epsilon_{N}^{\rm HHK}|$ versus $N$, where
$\Delta\epsilon_{N}=\epsilon_{N}^{\rm HHK}-\widetilde{\epsilon}_{N}^{(1)}$.
The agreement is very good with a relative deviation smaller than $2.3\times
10^{-4}$ for $|N|<25$. The results indicate that
$\Delta\epsilon_{N}\propto|N|^{5/2}$, a scaling corresponding to terms related
to $q^{3}$ that have been neglected in the momentum expansion, see App. A.

Figure 4(b) contrasts the analytical $\widetilde{\epsilon}_{N}^{(2)}$ second
nearest neighbor contribution to the LL spectrum with its numerical
counterpart $\epsilon_{N}^{(2)}$. The later is defined as
$\epsilon_{N}^{(2)}\equiv\epsilon_{N}-\epsilon_{N}^{(1)}$, where
$\epsilon_{N}$ and $\epsilon_{N}^{(1)}$ are computed numerically. Both are
obtained using the tight-binding parameter set B, but for the calculation of
$\epsilon_{N}^{(1)}$ we set $t^{(2)}=0$ to single-out the second-nearest
neighbor contribution. The insert shows $|\Delta\epsilon_{N}/\epsilon_{N}^{\rm
HHK}|$ versus $N$, $\Delta\epsilon_{N}=\epsilon_{N}^{\rm
HHK}-\widetilde{\epsilon}_{N}^{(1)}$. Here $\Delta\epsilon_{N}$ is
approximately $\sim 1.01\times 10^{-3}$ eV for $|N|\lesssim 25$. The origin of
$\Delta\epsilon^{(2)}$ is small discrepancy between the constant diagonal term
that shifts the energy spectrum by $-3t^{(2)}$ ($t^{(2)}=0.3$ eV, see Table 1)
in the analytical approach and the HHK results where the shift is $-0.899$ eV.

Finally, Fig. 4 (c) compares $\widetilde{\epsilon}_{N}^{(3)}$ with
$\epsilon^{(3)}$. Here, we define
$\epsilon_{N}^{(3)}=\epsilon_{N}-\epsilon_{N}^{(1)}-\epsilon_{N}^{(2)}$, where
we use the HHK method to compute $\epsilon_{N}$ with the parameter set C and
$\epsilon_{N}^{(1)}+\epsilon_{N}^{(2)}$ by taking $t^{(3)}=0$. The results
indicate that the third nearest neighbor contributions to $\epsilon_{N}$ are
nicely described by the leading $|N|^{1/2}$ term. Here, in distinction to the
first-nearest neighbors, the deviation is $\Delta\epsilon\sim|N|^{3/2}$ but
still small, $|\Delta\epsilon_{N}/\epsilon_{N}^{\rm HHK}|\leq 5\times 10^{-2}$
for $|N|\leq 25$.

Let us now compare our results with the experimental data[18]. Figure 5 (a)
shows that, since to leading order $\epsilon_{N}\sim|N|^{1/2}$, all tight-
binding parameterizations correctly capture the main trend of $\epsilon_{N}$.
In Fig. 5 (b) and 5 (c) we show $\epsilon_{N}^{\rm HHK}-\epsilon_{N}^{\rm
exp}$ 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
$|\epsilon_{N}^{\rm HHK}-\epsilon_{N}^{\rm exp}|\propto|N|^{3/2}$. This
explains the effectiveness on the extra arbitrary “$\mathcal{W}$” parameter
used in Ref.[18, 4] 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 $\epsilon_{\bf k}$ for $B=0$ and not to the LL
energies $\epsilon_{N}$.

Figure 5: (a) The Landau levels spectrum as a function of $\sqrt{|N|}$ for the
tight-binding parameterizations A, B and C (see Table 1) and $B=25$ T. The
experimental data are obtained from Ref. [18]. Deviation between the computed
and the experimental LL energies, $\epsilon_{N}^{\rm HHK}-\epsilon_{N}^{\rm
exp}$, for the (b) conduction and (c) valence bands.

###  III.2 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 $\Gamma_{N}$
but our main interest is the energy shift $\Delta_{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 $\eta=0.5$ meV as regularization parameter. This choice
guarantees that $\eta\ll\Gamma_{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 self-energy. We have introduced a
Debye-like cutoff $N_{c}=1342$ and start the self-consistent loop with an
initial self-energy of ${\rm Im}\;\Sigma(\epsilon)=1$ meV. The iteration
converges very rapidly.

Figure 6 compares the HHK with the SCBA results for the DOS($\epsilon$) of
graphene at $B=25$ T for different disorder strengths, namely, $\langle\langle
u_{i}^{2}\rangle\rangle\approx 0.02$ and $0.08$. The Landau subbands are
broadened and shifted toward the zero-energy in the presence of on-site
disorder with the exception of the zero-energy subband which is only
broadened, as nicely shown by Figs. 6(a) and 6(b).

Figure 6: DOS (in arbitrary units) versus $\epsilon$ (in eV) of disordered
graphene monolayers for the tight-binding parameterization A and $B=25$ T
obtained by the HHK method and the SCBA for (a) $\langle\langle
u_{i}^{2}\rangle\rangle\approx 0.02$ and (b) $0.08$. Blow up of the subband
DOS for (c) $\langle\langle u_{i}^{2}\rangle\rangle\approx 0.02$, $N=0$; (d)
$\langle\langle u_{i}^{2}\rangle\rangle\approx 0.02$, $N=10$; and (e)
$\langle\langle u_{i}^{2}\rangle\rangle\approx 0.08$, $N=0$.

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$(\epsilon)$ 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\gg 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 $\Delta_{N}$. We define $\Delta_{N}^{\rm HHK}$ as the difference
between the ensemble average $\langle\epsilon_{N}^{\rm HHK}\rangle$ and the
pristine value $\epsilon_{N}^{\rm HHK}$. The $\Delta_{N}^{\rm SCBA}$ 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, $\Delta_{N}$
is very small.

Figure 7: Renormalization of the LL subband peak energy $\Delta_{N}$ as
function of $N$ for $\langle\langle u_{i}^{2}\rangle\rangle\approx 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].

##  IV Conclusions and discussion

In this paper we have studied the large $|N|$ Landau level spectrum of
graphene monolayers beyond the first-nearest 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({\cal 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|\leq 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 tight-binding 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 $\Delta_{N}$. For realistic disorder strengths,
$\Delta_{N}$ gives a small correction to $\epsilon_{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].

###### Acknowledgements.

This work was partially supported by the Brazilian Institute of Science and
Technology (INCT) in Carbon Nanomaterials and the Brazilian agencies CAPES,
CNPq, FAPEMIG, and FAPERJ.

##  Appendix A Landau levels in graphene: Analytical approach

As mentioned i the main text, the low energy spectrum of graphene is obtained
by expanding the phase factors $\gamma$ and $\gamma^{\prime}$ in powers of q
up to third order in q. For the series expansion, it is convenient[43] to use
$\gamma\rightarrow\gamma\exp{\left(i{\bf k}\cdot\bm{\delta}_{3}\right)}$ and
$\gamma^{\prime}\rightarrow\gamma^{\prime}\exp{\left(i2{\bf
k}\cdot\bm{\delta}_{3}\right)}$ where
$\bm{\delta}_{3}=-a_{0}\hat{\textbf{e}}_{y}$. These changes do not affect the
dispersion relation and render more symmetric expressions to work with [43].
Thus, the expansions of the phase factors read

$\displaystyle\gamma^{\xi}_{\bf q}\approx$
$\displaystyle-\xi\dfrac{3a_{0}}{2}\left(q_{x}+i\xi q_{y}\right)$
$\displaystyle+\dfrac{3a^{2}_{0}}{8}\left(q_{x}-i\xi
q_{y}\right)^{2}+\xi\dfrac{3a^{3}_{0}}{16}|{\bf q}|^{2}\left(q_{x}+i\xi
q_{y}\right),$ (23) $\displaystyle\gamma^{\prime\xi}_{\bf q}\approx$
$\displaystyle-\xi 3a_{0}\left(q_{x}+i\xi q_{y}\right)$
$\displaystyle+\dfrac{3a^{2}_{0}}{2}\left(q_{x}-i\xi
q_{y}\right)^{2}+\xi\dfrac{3a^{3}_{0}}{2}|{\bf q}|^{2}\left(q_{x}+i\xi
q_{y}\right),$ (24)

where $\xi=\pm 1$ denotes the valley index and $|{\bf q}|a_{0}\ll 1$
guarantees the accuracy of the approximation.

The energy dispersion relation in power of $|{\bf q}|a_{0}$ is obtained by
combining the relations (A), (A) and (6). As a result, the energy dispersion
up to third-nearest neighbors hopping is conveniently written as

$\epsilon_{\lambda\xi\textbf{q}}=\epsilon_{\lambda\xi\textbf{q}}^{(1)}+\epsilon_{\xi\textbf{q}}^{(2)}+\epsilon_{\lambda\xi\textbf{q}}^{(3)},$
(25)

where

$\displaystyle\epsilon_{\lambda\xi\textbf{q}}^{(1)}$ $\displaystyle\approx$
$\displaystyle\lambda\hslash v^{(1)}_{F}|{\bf q}|\left[1-\xi\dfrac{a_{0}|{\bf
q}|}{4}\cos{\left(3\varphi_{\bf q}\right)}\right],$ (26)
$\displaystyle\epsilon_{\xi\textbf{q}}^{(2)}$ $\displaystyle\approx$
$\displaystyle\hslash
v^{(1)}_{F}\dfrac{t^{(2)}}{|t^{(1)}|}\dfrac{3a_{0}}{2}|{\bf
q}|^{2}\left[1-\xi\dfrac{a_{0}|{\bf q}|}{2}\cos{\left(3\varphi_{\bf
q}\right)}\right],\;\;\;\;$ (27)
$\displaystyle\epsilon_{\lambda\xi\textbf{q}}^{(3)}$ $\displaystyle\approx$
$\displaystyle\lambda\hslash v^{(1)}_{F}\dfrac{2t^{(3)}}{t^{(1)}}|{\bf
q}|\bigg{[}\cos{\left(2\varphi_{\bf q}\right)}+\dfrac{t^{(3)}}{t^{(1)}}$ (28)
$\displaystyle-\xi a_{0}|{\bf q}|\left(\dfrac{3}{4}\cos{\left(\varphi_{\bf
q}\right)}+\dfrac{t^{(3)}}{t^{(1)}}\cos{\left(3\varphi_{\bf
q}\right)}\right)\bigg{]}.\;\;\;\;$

The term $\cos(3\varphi_{\textbf{q}})$ is remnant of the symmetry of the
underlying lattice and gives rise to the so-called trigonal warping [3, 4]. We
note that the third nearest neighbors hopping terms modify the standard Fermi
velocity $v^{(1)}_{F}=3a_{0}t^{(1)}/2\hslash$.

We include the magnetic field in the Dirac Hamiltonian by minimal
substitution, that is, by replacing the canonical momentum
$\textbf{p}=\hslash{\bf q}$ by the kinetic momentum,
$\boldsymbol{\Pi}=\textbf{p}+e\textbf{A}(\textbf{r})$. The ansatz remains
accurate as along as the lattice spacing is much smaller than the magnetic
length, $\ell_{B}=\sqrt{\hslash/eB}$, a condition fulfilled even by the most
intense magnetic field currently produced in laboratory [4]. It is convenient
to express the operator $\boldsymbol{\Pi}$ in terms of harmonic-oscillator
ladder operator namely [18, 4]

$\widehat{a}=\dfrac{\ell_{B}}{\sqrt{2}\hslash}\left(\Pi_{x}-i\Pi_{y}\right)\;\;\;{\rm
and}\;\;\;\widehat{a}^{\dagger}=\dfrac{\ell_{B}}{\sqrt{2}\hslash}\left(\Pi_{x}+i\Pi_{y}\right)$
(29)

with $\left[\widehat{a},\widehat{a}^{\dagger}\right]=1$. With the help of
standard canonical quantization rules, one writes the Dirac Hamiltonian in the
presence of an external magnetic field $\textbf{B}=B\hat{\rm\textbf{e}}_{z}$
as [18]

$H_{B}^{\xi}\equiv\left(\begin{array}[]{cc}h^{\prime}_{\xi}&h_{\xi}^{\dagger}\\\
h_{\xi}&h^{\prime}_{\xi}\end{array}\right),$ (30)

where the diagonal element is

$\displaystyle h^{\prime}_{\xi}$ $\displaystyle\thickapprox$
$\displaystyle\hslash\omega_{c}\dfrac{t^{(2)}}{|t^{(1)}|}\dfrac{3}{\sqrt{2}}\dfrac{a_{0}}{\ell_{B}}\Bigg{[}\widehat{a}^{\dagger}\widehat{a}-\xi\dfrac{1}{2\sqrt{2}}\dfrac{a_{0}}{\ell_{B}}\left(\widehat{a}^{\dagger
3}+\widehat{a}^{3}\right)$ (31)
$\displaystyle\qquad-\>\dfrac{3}{8}\dfrac{a_{0}^{2}}{\ell_{B}^{2}}\left(\widehat{a}^{\dagger}\widehat{a}\right)^{2}\Bigg{]}$

and the off-diagonal ones are

$\displaystyle h_{\xi}$ $\displaystyle\thickapprox$
$\displaystyle\xi\hslash\omega_{c}\Bigg{[}\left(\widehat{a}^{\dagger}+\dfrac{2t^{(3)}}{t^{(1)}}\widehat{a}\right)-\xi\dfrac{a_{0}}{\sqrt{2}\ell_{B}}\left(\dfrac{\widehat{a}^{2}}{2}+\dfrac{2t^{(3)}}{t^{(1)}}\widehat{a}^{\dagger
2}\right)$ (32)
$\displaystyle\qquad-\>\dfrac{a_{0}^{2}}{\ell_{B}^{2}}\left(\dfrac{1}{4}\widehat{a}^{\dagger
2}\widehat{a}+\dfrac{2t^{(3)}}{t^{(1)}}\widehat{a}^{\dagger}\widehat{a}^{2}\right)\Bigg{]},$

where $\omega_{c}=\sqrt{2}v^{(1)}_{F}/\ell_{B}$ is the cyclotron frequency.

Here, we are interested in the $|N|\gg 1$ limit. Thus, to a good approximation
we can neglect $1/N$ corrections related to the ordering of the ladder
operators[4, 18] and approximate
$\left(\widehat{a}^{\dagger}\widehat{a}\right)^{2}\approx\widehat{a}^{\dagger
2}\widehat{a}^{2}\approx\widehat{a}^{2}\widehat{a}^{\dagger 2}$.

Hence, solving the eigenvalue equation in two-spinors
$H_{B}^{\xi}\psi_{N}=\epsilon_{N}\psi_{N}$, where
$\psi_{N}=\left(u_{N},v_{N}\right)^{T}$, the eigenvalue equation for the
second spinor component becomes [4]

$h_{\xi}h_{\xi}^{\dagger}v_{N}\simeq(\epsilon_{N}-h^{\prime}_{\xi})^{2}v_{N}.$
(33)

There are terms in $h^{\prime}_{\xi}$ and $h_{\xi}h_{\xi}^{\dagger}$ that do
not commute with $\widehat{N}=\widehat{a}^{\dagger}\widehat{a}$. One solves
(33) by perturbation theory, that is justified by noting that the corrections
are proportional to $a_{0}/\ell_{B}$ small parameter. With this treatment,
trigonal warping is taken into account at leading order. After some algebra,
we can write the energies of the LLs in the large $|N|$ limit up to third-
nearest neighbors hopping as

$\epsilon_{N}=\epsilon_{N}^{(1)}+\epsilon_{N}^{(2)}+\epsilon_{N}^{(3)}.$ (34)

Note that $\epsilon_{N}$ is independent of the valley $\xi$ index. The first
two terms

$\displaystyle\epsilon_{N}^{(1)}$ $\displaystyle\approx$ $\displaystyle{\rm
sgn}(N)\hslash\omega_{c}|N|^{1/2}\left(1-\dfrac{3}{8}\dfrac{a_{0}^{2}}{\ell_{B}^{2}}|N|\right),$
(35) $\displaystyle\epsilon_{N}^{(2)}$ $\displaystyle\approx$
$\displaystyle\hslash\omega_{c}\dfrac{t^{(2)}}{|t^{(1)}|}\dfrac{3}{\sqrt{2}}\dfrac{a_{0}}{\ell_{B}}|N|\left(1-\dfrac{3}{4}\dfrac{a_{0}^{2}}{\ell_{B}^{2}}|N|\right),$
(36)

were obtained in Refs. [18, 4]. Our derivation expands the latter results to
account for $t^{(3)}$, namely

$\epsilon_{N}^{(3)}\approx-{\rm
sgn}(N)\hslash\omega_{c}\dfrac{2t^{(3)}}{t^{(1)}}|N|^{1/2}\left(1-\dfrac{t^{(3)}}{t^{(1)}}-\dfrac{59}{32}\dfrac{a_{0}^{2}}{\ell_{B}^{2}}|N|\right).$
(37)

## References

  * Novoselov  _et al._ [2005] K. S. Novoselov, A. K. Geim, S. V. Morozov, D. Jiang, M. I. Katsnelson, I. V. Grigorieva, S. V. Dubonos, and A. A. Firsov, Two-dimensional gas of massless Dirac fermions in graphene, Nature 438, 197 (2005).
  * Zhang  _et al._ [2005] Y. Zhang, Y. Tan, H. Stormer, and K. P, Experimental observation of the quantum Hall effect and Berry’s phase in graphene, Nature 438, 201–204 (2005).
  * Castro Neto  _et al._ [2009] A. H. Castro Neto, F. Guinea, N. M. R. Peres, K. S. Novoselov, and A. K. Geim, The electronic properties of graphene, Rev. Mod. Phys. 81, 109 (2009).
  * Goerbig [2011] M. O. Goerbig, Electronic properties of graphene in a strong magnetic field, Rev. Mod. Phys. 83, 1193 (2011).
  * Yin  _et al._ [2017] L.-J. Yin, K.-K. Bai, W.-X. Wang, S.-Y. Li, and L. He, Landau quantization of dirac fermions in graphene and its multilayers, Front. Phys. 12, 127208 (2017).
  * Williams  _et al._ [2007] J. R. Williams, L. DiCarlo, and C. M. Marcus, Quantum hall effect in a gate-controlled $p-n$ junction of graphene, Science 317, 638 (2007).
  * Carmier  _et al._ [2010] P. Carmier, C. Lewenkopf, and D. Ullmo, Graphene $n\text{$-$}p$ junction in a strong magnetic field: A semiclassical study, Phys. Rev. B 81, 241406(R) (2010).
  * Du  _et al._ [2009] X. Du, I. Skachko, F. Duerr, A. Luican, and E. Y. Andrei, Fractional quantum Hall effect and insulating phase of Dirac electrons in graphene, Nature 462, 192–195 (2009).
  * Dean  _et al._ [2011] C. Dean, A. Young, P. Cadden-Zimansky, L. Wang, H. Ren, K. Watanabe, T. Taniguchi, P. Kim, J. Hone, and K. Shepard, Multicomponent fractional quantum Hall effect in graphene, Nat. Phys. 7, 693–696 (2011).
  * Dean  _et al._ [2013] C. R. Dean, L. Wang, P. Maher, C. Forsythe, F. Ghahari, Y. Gao, J. Katoch, M. Ishigami, P. Moon, M. Koshino, T. Taniguchi, K. Watanabe, K. L. Shepard, J. Hone, and P. Kim, Hofstadter’s butterfly and the fractal quantum Hall effect in moiré superlattices, Nature 497, 598–602 (2013).
  * Hunt  _et al._ [2013] B. Hunt, J. D. Sanchez-Yamagishi, A. F. Young, M. Yankowitz, B. J. LeRoy, K. Watanabe, T. Taniguchi, P. Moon, M. Koshino, P. Jarillo-Herrero, and R. C. Ashoori, Massive Dirac fermions and Hofstadter butterfly in a van der Waals heterostructure, Science 340, 1427 (2013).
  * Goerbig [2022] M. O. Goerbig, From the integer to the fractional quantum Hall effect in graphene, arXiv:2207.03322 (2022).
  * McClure [1956] J. W. McClure, Diamagnetism of graphite, Phys. Rev. 104, 666 (1956).
  * Sadowski  _et al._ [2006] M. L. Sadowski, G. Martinez, M. Potemski, C. Berger, and W. A. de Heer, Landau level spectroscopy of ultrathin graphite layers, Phys. Rev. Lett. 97, 266405 (2006).
  * Jiang  _et al._ [2007] Z. Jiang, E. A. Henriksen, L. C. Tung, Y.-J. Wang, M. E. Schwartz, M. Y. Han, P. Kim, and H. L. Stormer, Infrared spectroscopy of Landau levels of graphene, Phys. Rev. Lett. 98, 197403 (2007).
  * Li  _et al._ [2009] G. Li, A. Luican, and E. Y. Andrei, Scanning tunneling spectroscopy of graphene on graphite, Phys. Rev. Lett. 102, 176804 (2009).
  * Song  _et al._ [2010] Y. J. Song, A. F. Otte, Y. Kuk, Y. Hu, D. B. Torrance, P. N. First, W. A. de Heer, H. Min, S. Adam, M. D. Stiles, A. H. MacDonald, and J. A. Stroscio, High-resolution tunnelling spectroscopy of a graphene quartet, Nature 467, 185 (2010).
  * Plochocka  _et al._ [2008] P. Plochocka, C. Faugeras, M. Orlita, M. L. Sadowski, G. Martinez, M. Potemski, M. O. Goerbig, J. N. Fuchs, C. Berger, and W. A. de Heer, High-energy limit of massless Dirac fermions in multilayer graphene using magneto-optical transmission spectroscopy, Phys. Rev. Lett. 100, 087401 (2008).
  * Charlier  _et al._ [1991] J.-C. Charlier, J.-P. Michenaud, X. Gonze, and J.-P. Vigneron, Tight-binding model for the electronic properties of simple hexagonal graphite, Phys. Rev. B 44, 13237 (1991).
  * Reich  _et al._ [2002] S. Reich, J. Maultzsch, C. Thomsen, and P. Ordejón, Tight-binding description of graphene, Phys. Rev. B 66, 035412 (2002).
  * Urban  _et al._ [2011] A. Urban, M. Reese, M. Mrovec, C. Elsässer, and B. Meyer, Parameterization of tight-binding models from density functional theory calculations, Phys. Rev. B 84, 155119 (2011).
  * Lherbier  _et al._ [2012] A. Lherbier, S. M.-M. Dubois, X. Declerck, Y.-M. Niquet, S. Roche, and J.-C. Charlier, Transport properties of graphene containing structural defects, Phys. Rev. B 86, 075402 (2012).
  * Jung and MacDonald [2013] J. Jung and A. H. MacDonald, Tight-binding model for graphene $\pi$-bands from maximally localized Wannier functions, Phys. Rev. B 87, 195450 (2013).
  * Linhart  _et al._ [2018] L. Linhart, J. Burgdörfer, and F. Libisch, Accurate modeling of defects in graphene transport calculations, Phys. Rev. B 97, 035430 (2018).
  * Mucciolo and Lewenkopf [2010] E. R. Mucciolo and C. H. Lewenkopf, Disorder and electronic transport in graphene, J. Phys. Condens. Matter 22, 273201 (2010).
  * Zhu  _et al._ [2009] W. Zhu, Q. W. Shi, X. R. Wang, J. Chen, J. L. Yang, and J. G. Hou, Shape of disorder-broadened Landau subbands in graphene, Phys. Rev. Lett. 102, 056803 (2009).
  * Haydock  _et al._ [1972] R. Haydock, V. Heine, and M. J. Kelly, Electronic structure based on the local atomic environment for tight-binding bands, J. Phys. C: Solid State Phys. 5, 2845 (1972).
  * Haydock  _et al._ [1975] R. Haydock, V. Heine, and M. J. Kelly, Electronic structure based on the local atomic environment for tight-binding bands: II, J. Phys. C: Solid State Phys. 8, 2591 (1975).
  * Haydock [1980a] R. Haydock, The recursive solution of the Schrödinger equation, in  _Solid State Physics_, Vol. 35, edited by H. Ehrenreich, F. Seitz, and D. Turnbull (Academic Press, 1980) pp. 215–294.
  * Fan  _et al._ [2021] Z. Fan, J. H. Garcia, A. W. Cummings, J. E. Barrios-Vargas, M. Panhans, A. Harju, F. Ortmann, and S. Roche, Linear scaling quantum transport methodologies, Phys. Rep. 903, 1 (2021).
  * Hsiao and Peterson [1988] J. S. Hsiao and K. L. Peterson, Haydock recursion-method calculations of the local density of states of NiSi2 and ${\mathrm{ni}}_{3}$si, Phys. Rev. B 38, 10911 (1988).
  * Woodruff  _et al._ [1987] D. M. Woodruff, S. M. Anlage, and D. L. Smith, Density-of-states calculations within the recursion method, Phys. Rev. B 36, 1725 (1987).
  * Büsser  _et al._ [2013] C. A. Büsser, G. B. Martins, and A. E. Feiguin, Lanczos transformation for quantum impurity problems in $d$-dimensional lattices: Application to graphene nanoribbons, Phys. Rev. B 88, 245113 (2013).
  * Pincak and Smotlacha [2013] R. Pincak and J. Smotlacha, Analogies in electronic properties of graphene wormhole and perturbed nanocylinder, Eur. Phys. J. B 86, 480 (2013).
  * Shimomura  _et al._ [2011] Y. Shimomura, Y. Takane, and K. Wakabayashi, Electronic states and local density of states in graphene with a corner edge structure, J. Phys. Soc. Jpn. 80, 054710 (2011).
  * Smotlacha  _et al._ [2012] J. Smotlacha, R. Pincak, and M. Pudlak, Electronic structure of disordered graphene with Green’s function approach, Phys. Lett. A 376, 3256 (2012).
  * Eremenko  _et al._ [2016] V. Eremenko, V. Sirenko, I. Gospodarev, E. Syrkin, S. Feodosyev, I. Bondar, and K. Minakova, Anisotropic behavior and inhomogeneity of atomic local densities of states in graphene with vacancy groups, J. Sci. Adv. Mater. Devices 1, 167 (2016).
  * Prodan and Kohn [2005] E. Prodan and W. Kohn, Nearsightedness of electronic matter, Proc. Natl. Acad. Sci. U.S.A. 102, 11635 (2005).
  * Haydock [1980b] R. Haydock, The recursive solution of the Schrödinger equation, Comput. Phys. Commun. 20, 11 (1980b).
  * Ahn  _et al._ [2018] S. J. Ahn, P. Moon, T.-H. Kim, H.-W. Kim, H.-C. Shin, E. H. Kim, H. W. Cha, S.-J. Kahng, P. Kim, M. Koshino, Y.-W. Son, C.-W. Yang, and J. R. Ahn, Dirac electrons in a dodecagonal graphene quasicrystal, Science 361, 782 (2018).
  * Yao  _et al._ [2018] W. Yao, E. Wang, C. Bao, Y. Zhang, K. Zhang, K. Bao, C. K. Chan, C. Chen, J. Avila, M. C. Asensio, J. Zhu, and S. Zhou, Quasicrystalline 30o twisted bilayer graphene as an incommensurate superlattice with strong interlayer coupling, Proc. Natl. Acad. Sci. U.S.A. 115, 6928 (2018).
  * Czycholl and Ponischowski [1988] G. Czycholl and W. Ponischowski, The recursion method for a two-dimensional electron system in a strong magnetic field, Z. Phys. B 73, 343 (1988).
  * Bena and Montambaux [2009] C. Bena and G. Montambaux, Remarks on the tight-binding model of graphene, New J. Phys. 11, 095003 (2009).
  * Peierls [1933] R. Peierls, Zur Theorie des Diamagnetismus von Leitungselektronen, Z. Physik 80, 763 (1933).
  * Cresti [2021] A. Cresti, Convenient Peierls phase choice for periodic atomistic systems under magnetic field, Phys. Rev. B 103, 045402 (2021).
  * Note [1] Here, LDOS$(j,\epsilon)$ is normalized to unit and $\rho(\epsilon)$ to $2{\cal A}_{BZ}$ where ${\cal A}_{BZ}$ is the Brillouin zone area.
  * Lado and Fernández-Rossier [2016] J. L. Lado and J. Fernández-Rossier, Landau levels in 2D materials using Wannier Hamiltonians obtained by first principles, 2D Mater. 3, 35023 (2016).
  * Brown [1964] E. Brown, Bloch electrons in a uniform magnetic field, Phys. Rev. 133, A1038 (1964).
  * Pereira  _et al._ [2011] A. L. C. Pereira, C. H. Lewenkopf, and E. R. Mucciolo, Correlated random hopping disorder in graphene at high magnetic fields: Landau level broadening and localization properties, Phys. Rev. B 84, 165406 (2011).
  * Shon and Ando. [1998] N. H. Shon and T. Ando., Quantum transport in two-dimensional graphite system, J. Phys. Soc. Jpn 67, 2421 (1998).
  * Peres  _et al._ [2006] N. M. R. Peres, F. Guinea, and A. H. Castro Neto, Electronic properties of disordered two-dimensional carbon, Phys. Rev. B 73, 125411 (2006).
  * Dóra and Thalmeier [2007] B. Dóra and P. Thalmeier, Magnetotransport and thermoelectricity in Landau-quantized disordered graphene, Phys. Rev. B 76, 035402 (2007).
  * Ostrovsky  _et al._ [2008] P. M. Ostrovsky, I. V. Gornyi, and A. D. Mirlin, Theory of anomalous quantum Hall effects in graphene, Phys. Rev. B 77, 195430 (2008).
  * Wegner [1983] F. Wegner, Exact density of states for lowest Landau level in white noise potential superfield representation for interacting systems, Z. Phys. B. 51, 279 (1983).
  * Benedict and Chalker [1986] K. A. Benedict and J. T. Chalker, An exactly solvable model of the disordered two-dimensional electron gas in a strong magnetic field, J. Phys. C: Solid State Physics 19, 3587 (1986).
  * Carra  _et al._ [1989] P. Carra, J. Chalker, and K. A. Benedict, Electron motion in a random potential and strong magnetic field: Scattering and quantum interference in high Landau levels, Ann. Phys. (N. Y.) 194, 1 (1989).
  * Efetov and Marikhin [1989] K. B. Efetov and V. G. Marikhin, Effect of impurities on electron motion in two dimensions in a strong magnetic field, Phys. Rev. B 40, 12126 (1989).
  * Dmitriev  _et al._ [2012] I. A. Dmitriev, A. D. Mirlin, D. G. Polyakov, and M. A. Zudov, Nonequilibrium phenomena in high Landau levels, Rev. Mod. Phys. 84, 1709 (2012).
  * Kempkes  _et al._ [2019] S. Kempkes, M. R. Slot, S. Freeney, S. J. M. Zevenhuizen, D. Vanmaekelbergh, I. Swart, and C. Morais Smith, Design and characterization of electrons in a fractal geometry, Nat. Phys. 15, 127–131 (2019).