{'fraction_non_alphanumeric': 0.03222408907390225, 'fraction_numerical': 0.014971035935501668, 'mean_word_length': 4.310693834065788, 'pattern_counts': {'":': 0, '<': 12, '<?xml version=': 0, '>': 5, 'https://': 0, 'lorem ipsum': 0, 'www.': 0, 'xml': 0}, 'pii_count': 0, 'substrings_counts': 0, 'word_list_counts': {'cursed_substrings.json': 2, 'profanity_word_list.json': 0, 'sexual_word_list.json': 0, 'zh_pornsignals.json': 0}}
Landau level splitting due to graphene superlattices

The Landau level spectrum of graphene superlattices is studied using a tight-binding approach. We consider noninteracting particles moving on a hexagonal lattice with an additional one-dimensional superlattice made up of periodic square potential barriers, which are oriented along the zigzag or along the armchair directions of graphene. In the presence of a perpendicular magnetic field, such systems can be described by a set of one-dimensional tight-binding equations, the Harper equations. The qualitative behavior of the energy spectrum with respect to the strength of the superlattice potential depends on the relation between the superlattice period and the magnetic length. When the potential barriers are oriented along the armchair direction of graphene, we find for strong magnetic fields that the zeroth Landau level of graphene splits into two well-separated sublevels, if the width of the barriers is smaller than the magnetic length. In this situation, which persists even in the presence of disorder, a plateau with zero Hall conductivity can be observed around the Dirac point. This Landau level splitting is a true lattice effect that cannot be obtained from the generally used continuum Dirac-fermion model.


I. INTRODUCTION
Recently, the electronic and transport properties of graphene superlattices have been the subject of intense investigation. Theoretically, it was shown that the presence of periodic electrostatic [1][2][3][4][5][6][7][8][9][10][11][12][13][14] or vector [15][16][17][18][19][20] potentials, and also of periodic arrays of corrugations [21][22][23][24][25] tailors the graphene properties in a unique way, leading to novel features and interesting physics. In onedimensional superlattices, i.e., two-dimensional (2D) superlattice potentials depending on only one spatial direction, the Dirac cones of graphene are distorted, and hence the velocity of a particle moving parallel to the potential steps is reduced. 2 Moreover, for certain superlattice parameters, this component of the velocity is suppressed, and the carriers move only perpendicular to the potential steps of the superlattice. Furthermore, for other specific superlattice parameters, extra Dirac cones [5][6][7] and even Dirac lines 8 appear in the energy spectrum of graphene besides the usual K or K ′ Dirac points that exist in the continuum model of graphene at the neutrality point in the absence of a superlattice.
Interestingly enough, the emergence of the new Dirac points is controlled by the ratio between the potential amplitude and the superlattice period, irrespective of the superlattice profile, e.g., a cosine 5 or a Kronig-Penney 6,7 type, as long as the period is much larger than graphene's lattice constant. The extra Dirac points and their associated zero-energy modes 5 drastically affect the transport properties 5,8-10 of the system and also the Landau level sequence, 6 and hence the plateaus in the quantum Hall conductivity when a magnetic field is applied.
The implementation of two-dimensional superlattices is another route to modulate the electronic properties of graphene. For example, for rectangular superlattices, the velocity of carriers is also anisotropic and, depending on the Fermi level, the charge carriers are electrons, holes, or a mixture of both. 2 Recently, it has been shown that for two-dimensional rectangular superlattices the conductivity is unchanged from the result of pristine graphene, even if the velocity renormalization induced by the superlattice is quite large. 13 Also, new Dirac points with and without energy gaps can emerge at high-symmetry points in the Brillouin zone in two-dimensional triangular superlattices. 3,12 Ab-initio studies of the electronic and magnetic properties of graphene-graphane superlattices have also been reported. 26,27 For instance, it was shown that the zig-zag or arm-chair orientation of the graphene-graphane interface has a significant impact on the electronic properties of the system. 28 Experimentally, there are different possibilities to fabricate graphene superlattices. For example, it is possible to imprint superlattice patterns with periodicity as small as 5 nm using electron-beam induced hydrocarbon lithography on graphene membranes. 29 Graphene grown epitaxially on Ru(0001) [30][31][32][33] or Ir(111) [34][35][36] surfaces, and also on SiC 37,38 shows two-dimensional superlattice patterns with lattice period of a few nanometers and potential strength in the range of few tenths of an electron volt. In suspended graphene, the existence of periodic ripples has been recently demonstrated. 39 Another possibility to make superlattices with controlled potential amplitude is to fabricate periodically patterned gates: p-n and p-n-p junctions in graphene [40][41][42][43][44] have already been realized.
In the present work, we study the evolution of the Landau levels of graphene appearing in the presence of a one-dimensional superlattice and a strong magnetic field applied perpendicular to the graphene plane. For the superlattice, we assume a Kronig-Penney type of electrostatic potential with alternating barriers of +V and −V potential strengths and barrier width L a . Here, two cases are considered, one with the potential barriers oriented along the zig-zag (zz) and the other oriented along the arm-chair (ac) directions of graphene. We use a tightbinding approach with nearest-neighbor hopping. The magnetic field is incorporated in the hopping integral through the Peierls phases as usual. Since the superlattice depends only on one direction, we use plane waves in the other direction. Thus, the two-dimensional problem is reduced to a set of one-dimensional equations, known as the Harper equations. 45 There are two characteristic lengths in the system: one is the superlattice barrier width L a and the other is the magnetic length l B . The behavior of the energy spectrum is governed by the relation between L a and l B . We find that if L a is greater than l B , the Landau levels acquire a finite broadening (irrespective of disorder) when the superlattice potential strength increases. In this case, the Landau level sequence disappears already for small values of V . In the other case, when L a is smaller than l B , the Landau levels picture survives for much higher values of the potential strength. Also, two or more Landau levels may merge together when V is increased. Consequently, the plateaus in the Hall conductivity show an unconventional sequence that may be controlled by either the external electrostatic potential or the applied magnetic field. When the barriers of the superlattice potential are along the arm-chair direction of graphene, and when L a < l B , we find that the zeroth Landau level splits into two sublevels, a novel signature of the interplay between the magnetic field and superlattice potential orientation. We show that the splitting is robust against different types of disorder affecting the superlattice potential. This holds also for a two-dimensional chess-board type superlattice.
The paper is organized as follows. In Section II we introduce the theoretical model that we use in order to investigate graphene superlattices in the presence of a strong magnetic field. In Section III we present the results concerning the energy spectrum of the system for various superlattice parameters and magnetic flux density strengths. A summary and concluding remarks are given in Section IV.

II. HARPER EQUATION FOR GRAPHENE RIBBONS
The Harper equations of graphene ribbons in a uniform magnetic field reduce the two-dimensional tightbinding problem to two one-dimensional equations due to the translational invariance along one axis. Derivations of the Harper equations can be found in Ref. 46 for a hexagonal lattice with zig-zag edges and in Ref. 47 for bricklayer lattices with zig-zag and with arm-chair edges. Here, we briefly derive Harper equations in the presence of a superlattice potential. From these, we calculate the electronic energy spectra of monolayer graphene ribbons with oriented edges in a perpendicular magnetic field and a superlattice potential.
We consider a one-band tight-binding model with near- est neighbor hopping on a hexagonal lattice. The structure of graphene ribbons of length L with zig-zag (zz) and arm-chair (ac) edges is shown in Fig. 1. The two interpenetrating triangular sublattices are denoted by A and B. The ribbon width W ≪ L (x direction) is defined by the number N zz of zz lines for the zz ribbon and by the number N ac of dimer lines for the ac ribbon where a is the distance between two neighboring carbon atoms (a = 1.41Å).
With ψ(r i ) the wavefunction amplitude on site r i , the Schrödinger equation reads where ε is the associated eigenvalue and V the electrostatic potential (superlattice). The sum runs over the nearest neighbors of atom i and t ij is the transfer integral between atoms i and j. The magnetic field perpendicular to the graphene plane is incorporated in the hopping term by means of the Peierls phase where t is the hopping parameter (t =2.7 eV) and γ ij is given by the line integral of the vector potential from site r i to site r j , and φ 0 = h/e is the magnetic flux quantum We use the Landau gauge A = (0, Bx, 0) and the translational invariant direction of each ribbon is taken as the y axis. With this particular choice of the gauge, the line integral of the vector potential between two nearest neighbors i and j becomes The magnetic flux φ = BS through the area S = 3a 2 √ 3/2 of a hexagon is taken to be a rational multiple of the flux quantum φ = (p/q)φ 0 , hence the magnetic flux density B and the magnetic length are set by the integers p and q which are chosen to be mutually prime: .
The one-dimensional superlattice potential is taken to be periodic along the x direction, with periodicity 2L a (barrier width L a ) and constant in the y direction We are ready now to write the Harper equations for the zz and the ac ribbon.

A. Harper equations for the zz ribbon
In the following, we label the atoms as shown in Fig. 1 (b) and use m to index the zz chains, where m takes values between 1 and N zz . Then, the Schrödinger equations for the atoms 1A and 2B are From Eq. (6), the Peierls phases are m and so the Harper equations take the simple form where we defined Here, we use x A m = (3m/2 − 1)a and x B m = (3/2)(m − 1)a. Equation (12) does not contain any boundary conditions yet. In the case of zz ribbons with finite N zz , Dirichlet boundary conditions are imposed by setting ψ For periodic boundary conditions (torus), one sets ψ Nzz and Eq. (12) becomes However, a more careful treatment is needed at the edge sites for periodic boundary conditions, due to the requirement for commensurability between the magnetic and the lattice lengths: the Peierls phases must be periodic with the ribbon width. This implies that and one finds the condition where n can be any nonzero positive integer. This means that the parameters for the magnetic flux density (p/q) and the ribbon width (given by N zz ) are not independent but must be chosen in such a way that Eq. (17) holds. In order to obtain the energy spectrum of the system, we diagonalize numerically the RHS of Eq. (12), i.e., a 2N zz × 2N zz matrix with the appropriate boundary conditions.

B. Harper equations for the ac ribbon
To obtain the Harper equations for ac ribbons, we proceed in a similar way as in the case of zz edges. Here, m indexes now the dimer lines of the ac ribbon, and takes values between 1 and N ac . The Schrödinger equations for the atoms 1A and 2B of Fig. 1 with the Peierls phases Then the Harper equations take the simple form where we defined . In the case of ac ribbons with periodic boundary conditions, commensurability between the lattice and the magnetic field requires that The energy spectrum of an ac ribbon is obtained by diagonalizing the RHS of Eq. (20), i.e., a 2N ac × 2N ac matrix.

III. RESULTS AND DISCUSSIONS
In this section, we present results obtained by solving the Harper equations corresponding to graphene superlattices under perpendicular magnetic fields. The superlattice potential is given by a Kronig-Penney function periodic along the x direction, with barrier width L a , and with alternating ±V barrier heights. The perpendicular magnetic field is set by the parameters p and q, according to Eq. (7). For simplicity, we take in the following p = 1.
As a test, we consider first the cases with no superlattice potential or no magnetic field. Then we show that, for V = 0 and B = 0, the energy spectrum and hence the system properties strongly depend on the relation between the superlattice barrier width L a and the magnetic length l B , as well as on the orientation of the superlattice barriers with respect to the zz or ac directions of graphene.
For infinite 2D systems, when the superlattice potential is not present (V = 0) and for high magnetic fields, the energy eigenvalues are the usual Landau level (LL) sequence which, according to the continuum model, occur at 48,49 where n is the LL index, and for the last equality we have used Eq. (7) and v F = 3/2ta.
The left panel of Fig. 2 (a) shows the energy spectrum in the 1 st Brillouin zone of a graphene ribbon in the absence of the superlattice potential. The low lying LL energies (shown in green, independent of k y ) appear according to Eq. (23). Higher energies (not shown) may deviate from the sequence because Eq. (23) is valid for an infinite system. We treat here ribbons with a finite width and are only interested in the low-energy domain. For Dirichlet boundary conditions, the edge states (k y dependent) are shown in red. Between two LLs, the Hall conductivity σ xy is equal to e 2 /h×(number of edge states) and experiences a jump every time the Fermi energy coincides with the LL energy. The corresponding Hall conductivity is shown in the right panel of Fig. 2 (a), with the quantized values of the integer quantum Hall effect of graphene. 50,51 In the absence of the magnetic field, and depending on the superlattice parameters, the energy spectrum of a graphene zz superlattice may exhibit additional Dirac points for E = 0. For a Kronig-Penney superlattice with barrier width L a and height V , the number of Dirac points increases by two (×4 considering valley and pseudospin degeneracy) whenever the potential amplitude exceeds a value of with n a positive integer. 6 Figure 2  One superlattice barrier contains 60 zz lines, which corresponds to a barrier width of L a = 12.7 nm. One clearly sees that when the barrier height V exceeds integer multiples of U 0 = 0.052 t (= 0.14 eV), then additional zeroenergy modes appear in the energy spectrum. Our spec- The electronic properties of graphene nanoribbons with ac edges strongly depend on their size. It is known that such systems are metallic when (N ac − 2)/3 ∈ Z and semiconductor otherwise. 52 Figure 3 (a) shows the energy bands of two typical ac ribbons, a semiconducting one with N ac = 600 that corresponds to W 600 ac = 73.6 nm and an energy gap E g = 0.0162 eV, and a metallic one with W 620 ac = 76.12 nm. A superlattice potential with L a = 1.22 nm (10 ac dimer lines under each barrier) with barriers along the ac edges is imposed on the systems and the evolution of the energy gap around E = 0 is calculated as a function of the potential strength. The results are shown in Fig. 3 (b). For the semiconduc- tor ac ribbon (N ac = 600), the energy gap is reduced from E g = 0.0162 eV when V = 0 to E g = 0 when V = 0.39 eV, and then increases again for higher values of the potential strength. For the metallic ac ribbon (N ac = 620) the superlattice potential opens a spectral gap that reaches a maximum value of E g = 0.022 eV for V = 0.45 eV and decreases for higher V . The respective minimal and maximal energy gaps E g depend both on N ac and L a .

C. Superlattice parallel to zig-zag edges
We now consider Kronig-Penney superlattices with potential barriers oriented along the zz direction of graphene in a strong perpendicular magnetic field. In our calculations we use graphene ribbons with N zz = 12000, which corresponds to a ribbon width of 2555 nm. The Kronig-Penney superlattice parameters are chosen in such a way that the number of barriers is even. Also, the number of magnetic flux quanta per graphene plaquette is fixed at p/q = 1/6000, which gives B = 13.15 T for the magnetic field and l B = 7.07 nm for the magnetic length. These settings allow us to study infinitely long (in the y direction) ribbons with Dirichlet or periodic boundary conditions in the transverse direction.
First, we discuss the case when the barrier width is larger than the magnetic length. Figure 4 shows the results for L a = 21.3 nm (100 zz lines per barrier) and l B = 7.07 nm. Typical energy spectra as function of the wavevector k y in the first Brillouin zone are given in Fig. 4 (a) for two different values of the potential strength V = 0.02 eV (left) and V = 0.05 eV (right). The difference between periodic and Dirichlet boundary conditions consists in the appearance of edge states when applying the latter. The oscillations seen in the LL energies correspond to mini-Brillouin zones imposed by the superlattice. The presence of the superlattice modifies the energies of the LLs, which acquire now a finite width depending on the value of V . Note that the finite bandwidth of the LLs is a consequence of applying a superlattice and not because of disorder, which is always present in experimental situations and induces additional broadening of the LLs. The width of the LLs increases with increasing V , and for larger values of the potential strength the Landau bands merge together in the sense that there remain no energy gaps between them. For example, in Fig. 4 (a) and for V = 0.05 eV, this is the case for the ±1 st and the ±2 nd LLs. Interestingly enough, the broadening of the LLs is not the same for all the LLs. For the parameters used in Fig. 4 the width of the ±3 rd LL is smaller than that of the other LLs, and this particular LL will merge with the others only for higher values of V . The origin of this individual broadening is discussed in detail below.
In Fig. 4 (b) we show the energies of the LLs as a function of the superlattice potential strength. The energy gaps between the LLs dissappear altogether in this case when V 0.07 eV. Also, for V < 0.07 eV we find that the bandwidth of the 0 th LL is equal to the value of V . When increasing the superlattice barrier width while L a ≫ l B , the energy bandwidth of the higher LLs also increases with the strength of the potential, and the merging of the LLs occurs for even smaller values of V .
The unusual structure of the energy bands is most directly reflected in the plateau sequence of the quantum Hall conductivity. Figure 4 (c) schematically shows the Hall conductivity versus the Fermi energy corresponding to the parameters used in panel (a), i.e., L a = 21.3 nm, l B = 7.07 nm. Comparing the cases V = 0.02 eV and 0.05 eV, one sees that with increasing V the reduction of the energy gaps between the LLs leads to a decrease of the plateau widths. Moreover, for V = 0.05 eV the plateaus at ±(3/2)(4e 2 /h) are not present, and the Hall conductivity has an unconventional step size. This is because the ±1 st and the ±2 nd LLs are now merged together, and there is no energy gap between them. Of course, the plateaus in the energy-dependent Hall conductance must not be confused with the experimentally observed Hall plateaus, which originate from localization due to the intrinsic disorder.

Landau level broadening
The different broadening of the distinct LLs labeled n can be explained using perturbation theory. We consider the superlattice potential as a perturbation to the graphene wavefunctions ψ 0 n belonging to a given energy. The integral gives the energy corrections up to first order to one of the n th LL energies as a function of k y . The values of |κ n | for different LLs provide a quantitative measure of the influence of the superlattice on the LL spectrum. The energies of the LLs with a larger |κ n | are more spread out than the energy levels corresponding to a smaller |κ n |, as illustrated in Fig. 5. Here, we consider a system with N zz = 120, p/q = 1/60, V = 0.135 eV and L a = 15 a. The energies of the lowest four LLs given in Fig. 5 (a) show that the 0th LL (E 0 ∼ 0 eV) is most broadened by the superlattice, and the broadening decreases successively for the 1 st LL (E 1 ∼ 1.13 eV), 2 nd (E 2 ∼ 1.56 eV) and 3 rd (E 3 ∼ 1.98 eV). The amplitudes of κ n reflect directly the broadenings of different LLs, as shown in Fig. 5 (b). For example, κ 1 has the largest amplitude in comparison with κ 2 and κ 3 , the latter has the smallest amplitude. Correspondingly, in the presence of a superlattice potential, the 1 st LL exhibits a larger broadening than the 2 nd , and the 3 rd LL shows the smallest broadening.
To understand the oscillations of the κ n integrals with k y , we analyze the spatial representation of the wavefunctions ψ 0 n of graphene in a perpendicular magnetic field. Figure 5 (c) shows |ψ 0 n (x)| 2 for n = 1 (left) and n = 2 (right) and for two values of k y , which correspond to zero (upper part) and maximum (lower part) values of the respective |κ n | integrals. Note that changing k y only shifts the wavefunctions with respect to the x axis but does not alter the structure of |ψ 0 n | 2 . The wavefunctions consist of two symmetrical contributions from the A and the B sublattices, and the two parts have a reflection symmetry axis (dashed line). The position of the reflection symmetry axis of the wavefunctions with respect to the superlattice potential barriers is crucial in determining the amplitude of the respective κ n integrals. In the upper parts of Fig. 5 (c), for k y = 2mπ/(6a √ 3), m = 0, 1, 2, . . . the symmetry axis of both the 1 st and the 2 nd LLs coincides with the superlattice barrier edge, where the potential changes sign from +V to −V . In this case, the contribution to the κ n integral of the wavefunctions situated to the left and to the right of the reflection symmetry axis cancel each other and κ n takes on a minimal value. On the contrary, for k y = (2m + 1)π/(6a √ 3) the reflection symmetry axis coincides with the middle of a superlattice potential barrier, the wavefunctions from the left and the right sides of the reflection symmetry axis contribute the same amount when taking the integral Eq. (25), leading to a maximal κ n . For all other k y , the values of κ n fall in between these two limiting cases.
The different amplitudes of κ n for different LLs can be also explained from the lower parts of Fig. 5 (c), by carefully examining the spatial distribution of the wavefunction with respect to the barrier width. For the 1 st LL, the contributions left and right of the reflection symmetry axis extend mainly over a single superlattice potential barrier, in this case +V . In the case of the 2 nd LL, the wavefunction extends over three barriers. When calculating κ 2 , one subtracts from the main contribution under the +V barrier the part of the wavefunction that are extended over the −V barriers. Hence, the amplitude of κ 1 is bigger than the amplitude of κ 2 and, correspondingly, the broadening of the 1 st LL is larger than the broadening of the 2 nd LL.
This approach was carried out for several system parameters to check its validity. We conclude that such an analysis provides an easy way to find out a-priori, starting only from the unperturbed wavefunctions of graphene in a magnetic field, which of the LLs will exhibit a large or a small broadening when a superlattice potential is switched on.
When L a ≈ l B , the lowest LLs survive for much higher values of the superlattice potential strength, as can be seen in Fig. 6 (left). Here, L a = 6.4 nm (30 zz lines per barrier) and l B = 7.07 nm. The merging of the 0 th and the ±1 st LLs occurs when V takes the value 3U 0 /2 (= 0.42 eV for L a = 6.4 nm), which depends only on the inverse of the superlattice barrier width according to (24). For V = 3U 0 /2 the zero-energy LLs becomes three-fold degenerate and, as the Fermi energy is scanned from negative to positive energies, the Hall conductivity presents a step size of 3(4e 2 /h). 6 The right side of Fig. 6 shows the evolution of the LLs as V is increased for a system with L a < l B . The width of the low energy LLs, although finite, is very small, and depends weakly on V . A general tendency is the bend- ing of the LL energies towards E = 0. Again, when V = 3U 0 /2 (= 1.05 eV for L a = 2.5 nm, which are 12 zz lines per barrier), the 0 th and the ±1 st LLs merge and the zero-energy LL becomes three-fold degenerate. The small widths of the LLs in the L a < l B case can be explained using a similar analysis of the κ n integrals from above. In this case, the wavefunctions of the LLs spread over many superlattice barriers. The contributions of the wavefunctions under adjacent barriers cancels out when calculating the κ n integrals, which leads to a very small broadening of the LLs. Our results obtained within the Harper equation method are consistent with the results from a perturbative approach starting from the continuum model for graphene. 53 There, it was also found that in the case of weak fields (i.e., L a < l B ) the matrix elements of the perturbation Hamiltonian do neither depend on the center nor on the spread of the LL wavefunctions, which leads to flat bands (i.e., small widths of the LLs).
We have performed several calculations for different magnetic field strengths and superlattice parameters, and the results shown here are most illustrative. The qualitative behavior of the LL energy spectrum in the three regimes, namely L a > l B , L a ≈ l B and L a < l B is robust against changing the system parameters, with the prerequisite that both the superlattice barrier width and the magnetic length are much larger than the graphene lattice constant.

D. Superlattice parallel to arm-chair edges
In the following, we discuss the case of graphene superlattices in a strong magnetic field with the potential barriers oriented along the ac direction. When no superlattice is present, the LL energy spectra for ribbons with zz or ac edges are identical if periodic boundary conditions are applied. That is because the LLs are a property of the bulk of the graphene ribbon and not of the edges. However, when a one-dimensional superlattice potential is switched on, then the system can be considered to consists of many 'sub-ribbons' with internal edges between them. Although we still investigate superlattice barrier widths much larger than the graphene lattice constant, we show below that the orientation of the sub-ribbon edges with respect to the zz or ac directions of graphene plays an essential role when the magnetic length l B is larger than the sublattice barrier width L a . Figure 7 shows the energy of the LLs of ac superlattices as a function of the potential strength for two different regimes: L a > l B (left) and L a ≈ l B (right). Here, p/q = 1/6000 (B = 13.15 T and l B = 7.07 nm), L a = 24.6 nm corresponds to 200 dimer lines per barrier (left) and L a = 7.3 nm corresponds to 60 dimer lines per barrier (right). The behavior of the LLs is qualitatively similar to the case of zz superlattices: for L a > l B the LLs acquire a large broadening and merge together when increasing the potential strength, and the LL sequence breaks down for even smaller values of V . Also, when L a ≈ l B , the LLs bend towards zero energy, and their broadening is not so strong, so that the LL picture survives for higher values of the superlattice potential strength. However, there are some significant quantitative differences between the ac and the zz cases. For instance, the merging of the 0 th and the ±1 st LLs does not occur at V = 3U 0 /2 any more, but at values of V which are more complicated to predict from a continuum model, as they do not depend only on the inverse of the superlattice barrier width.
For ac graphene superlattices with strong magnetic field, the L a < l B regime is the most interesting one. In this case, the electrons travel, in classical terms, over several potential barriers before closing a cyclotron radius, and the ac edge of each 'sub-ribbon' has an unexpected influence on the energy spectrum of the LLs. We find a splitting of the 0 th LL into two sub-bands when the potential barrier is increased. Figure 8 (a) illustrates this effect for an ac superlattice with L a = 1.47 nm (12 dimer lines per barrier) and l B = 7.07 nm. The 0 th LL splits as soon as V is turned on, and the energy difference between the two subbands increases continuously with V . For very large V , a splitting of the higher Landau bands was observed too (not shown). Note also that the higher LLs bend again towards zero energy, and their broadening is very weak. In the presence of Dirichlet boundary conditions, no edge states do appear within the energy range between the split Landau level.

E. Origin of Landau level splitting
A possible explanation for the splitting can be found by carefully examining the superlattice barrier potential step. For ac superlattices, the barrier step asymmetrically divides the graphene hexagons, with 4 atoms on one side and two next-neighbor atoms on the other side (see the inset of Fig. 8 (a)). In the case of zz superlattices, where the splitting does not occur, the superlattice barrier symmetrically divides the graphene hexagons, with 3 carbon atoms under one barrier and the other 3 atoms under the next barrier with opposite sign. We have considered zz superlattices with an artificially imposed asymmetry, realized by placing an atom from each divided hexagon under the adjacent barrier, so that a 4-2 asymmetry is created for all hexagons associated with a barrier potential step. This configuration is schematically shown in the upper inset of Fig. 8 (b). Now, the energy spectrum clearly shows the splitting of the 0 th LL which is, however, not as strong in magnitude as in the case of ac superlattices. Another configuration of an artificial superlattice with the barriers oriented along the zz direction that divides the step hexagons into 4 atoms on one side and a pair of next-neighbor atoms on the other side, is shown in the lower inset of Fig. 8 (b). Again, the 0 th LL is split into two sub-bands as the strength of the potential V is increased, with a splitting magnitude larger compared to the previous case.
These examples point already to the origin of this new effect. It is the absence of inversion symmetry due to the combined influence of magnetic field and superlattice potential that matters in the ac case. And in the zz situations discussed above, the artificially imposed asymmetry that is responsible for the splitting destroys the inversion symmetry as well. Once more, a closer look at the wavefunctions in the absence of the superlattice together with the application of a (degenerate) perturbation theory treatment up to second order in the superlattice potential provides a microscopic explanation for For comparison, the splitting of an ac superlattice with L ac a = 1.23 nm is also shown (blue data, labeled 2). The insets in Fig. 8 (b) schematically show the artificial barrier edges of two zz superlattices with broken symmetry that divide the graphene hexagons into 2 and 4 carbon atoms.
the LL splitting. In both the zz and ac situations, the Landau levels get broadened in first order perturbation theory due to the superlattice potential. A splitting of the LL, however, is seen in second order only in the ac case because of the lacking inversion symmetry. The latter is still present in the zz situation (without an artificial symmetry breaking) so that energy levels remain degenerate. Thus, the ac superlattice induced LL splitting is the generic case, whereas the particular inversion symmetry present in the perfect zz orientation turns out to be only a special situation, probably not met in real samples. We conclude that this LL splitting is a true lattice effect that was not seen before. The splitting is not to be found in the usual continuum model descriptions, since there, one normally cannot distinguish between ac or zz oriented superlattices. The splitting for ac superlattices is still present when changing the system parameters, i.e., B, L a , and V , as long as we are in the L a < l B regime. Interestingly enough, we find that the splitting of the 0 th LL is maximal when L a is chosen such that the number of ac dimer lines under each superlattice barrier is a multiple of 3. Moreover, in such cases the energies of the split 0 th LL do not change with the magnetic field and depend only on V and L a . For other values of L a , when the number of ac dimers lines under each barrier is not a multiple of 3, the splitting is still present, but the energy difference between the split sub-bands is one order of magnitude smaller than in the previous case and depends on the strength of the magnetic field as ∝ B 2 .
Experimentally, ideal rectangular superlattices with sharp potential jumps can hardly be fabricated. To verify that the splitting occurs also when the barrier edges are smooth, we have performed calculations considering a width w between two adjacent barriers where the superlattice potential changes linearly from +V to −V . We find that when increasing w, the energy difference between the split subbands of the 0 th LL decreases, but the splitting is still present even for w = L a /2.

G. Influence of disorder
To make sure that the splitting is not induced by an artificial hidden symmetry of the Hamiltonian, we have considered different kinds of superlattice disorder in our calculations. First, we have studied the case of barrier width disorder, by allowing 10% of the barriers to have the width L a ± δL a . The results are shown in Fig. 9 (a) for a system with l B = 4.08 nm, L a = 1.47 nm, and V = 0.25 eV. Increasing the value of δL a from 0.12 nm to 0.36 nm results in a broadening of the split 0 th LL. When the value of δL a is further increased, the broadening overlaps the splitting. The next case of disorder considered is barrier height disorder. The potential strengths of the superlattice barriers are allowed to take random values in an interval δV around ±V , as is schematically shown in the inset of Fig. 9 (b). Again, we find that this type of disorder does not destroy the splitting, and only induces an additional broadening of the sublevels, which grows with increasing disorder strength. Figure 9 (b) shows the energy band in the first Brillouin zone of the split 0 th LL for a system with L a = 1.47 nm, l B = 7.07 nm and V = 0.3 eV. Here, the disorder strengths is δV = 0.03, 0.08 and 0.13 eV, respectively. Shown are the results for 10 different disorder realizations.
Finally, we have considered ac superlattices with rough edges, as schematically shown in the inset of Fig. 9 (c). In this case, the superlattice potential depends on the y coordinate, and the Harper equation approach cannot be used. Therefore, we have directly diagonalized the tight-binding Hamiltonian for graphene with superlattice potential along the ac direction and with rough edges. Figure 9 (c) shows the energy of the split 0 th LL for ac superlattices with rough edges, compared with the case of perfect edges. In this case, 1.25% of the atoms along one barrier edge are taken at random and forced to have a potential strength equal to the one of the next barrier. We show the results for 100 disorder realizations and conclude that the splitting is still present in the spectrum.
The splitting of the 0 th LL leads to the prediction of a σ xy = 0 plateau to occur in the quantum Hall conductiv- ity, which may be regarded as a hallmark of ac oriented graphene superlattices. Most interesting, the splitting occurs only when L a < l B . Because B ∼ l −2 B , this corresponds to magnetic fields that are strong enough to have LLs, but low enough to be in the L a < l B regime. When increasing B, the magnetic length is reduced and the splitting disappears for L a > l B . Here, we have presented the results for L a = 1.5 nm, V ranging between 0 and 0.85 eV, and l B = 7.0 nm which corresponds to B = 13.1 T and are realistic parameters. Experimentally achievable are L a between 2 − 10 nm, and l B = 5 − 20 nm (B = 2 − 25 T), and superlattice potential strengths of a few tenths of an electronvolt. Thus it should be possible to check our findings experimentally.

IV. SUMMARY
We have investigated the Landau level structure of single layer graphene in the presence of a one-dimensional electrostatic square-potential superlattice. The superlattice modulates the Landau level spectrum in a unique way that is most directly reflected in a peculiar band broadening, in band bendings, and in an unusual plateau sequence of the quantum Hall conductivity. We have shown that, depending on the magnitude of the superlattice barrier width L a with respect to the magnetic length l B , the energy band structure of the LLs changes dramatically. In general, when L a > l B , the LLs are quickly broadened and merge together for small values of the superlattice potential strengths. In the opposite regime, when L a < l B , the LLs survive for much higher values of V , with a general tendency of the higher-energy LLs to bend towards zero energy.
The orientation of the superlattice barriers with respect to the zz or ac directions of graphene also plays a crucial role when L a < l B . That is because the superlattice divides the system into many sub-ribbons of width L a that become decoupled when increasing the strength of the superlattice potential V . When the magnetic length is larger than the sub-ribbon width, an electron travels over many barriers before it closes a cyclotron radius. Therefore, the results differ depending on the sub-ribbon edge type the electron encounters at each barrier jump. When the sub-ribbons have ac edges, we found a novel effect that originates from the interplay between graphene's hexagonal lattice and the additional superlattice potential barriers and can, therefore, not be found in the usual Dirac-fermion continuum model description. The new observation is the splitting of the zeroth LL which occurs with increasing the superlattice potential strength. Alternatively, one can tune the magnetic flux density instead and keep the superlattice potential strength fixed. This intriguing effect is linked to the absence of inversion symmetry in the ac case due to the presence of both a superlattice potential and the magnetic field. The parameters that we used here are experimentally accessible, and the peculiar features of the electronic structure may be tested directly in transport or optical experiments.
Finally we mention that according to further calculations (results not shown), a splitting of the zeroth Landau level occurs also in the presence of truly two-dimensional (chess-board type) superlattices, which remains robust even in the presence of additional on-site disorder. This observation opens the route to consider the charge density fluctuations 54,55 occurring in real graphene samples on a length scale between 20 nm and 30 nm as a natural 2D-superlattice that replaces the artificial one investigated in our model calculations. We then suggest that if the magnetic field is tuned such that the size of the charge 'puddles' and the magnetic length had the required ratio, the resulting gap opening could be responsible for an insulating behavior and the diverging resistance at the Dirac point with the accompanying σ xy = 0 Hall plateau as has been observed previously in experiments by Checkelsky et al. 56,57