# Graphene under spatially varying external potentials: Landau levels, magnetotransport, and topological modes

## Abstract

Superlattices (SLs) in monolayer and bilayer graphene, formed by spatially periodic potential variations, lead to a modified bandstructure with extra finite-energy and zero-energy Dirac fermions with tunable anisotropic velocities. We theoretically show that transport in a weak perpendicular (orbital) magnetic field allows one to not only probe the number of emergent Dirac points but also yields further information about their dispersion. or monolayer graphene, we find that a moderate magnetic field can lead to a strong reversal of the transport anisotropy imposed by the SL potential, an effect which arises due to the SL induced dispersion of the zero energy Landau levels. This effect may find useful applications in switching or other devices. For bilayer graphene, we discuss the structure of Landau level wave functions and local density of states in the presence of a uniform bias, as well as in the presence of a kink in the bias which leads to topologically bound ‘edge states’. We consider implications of these results for scanning tunneling spectroscopy measurements, valley filtering, and impurity induced breakdown of the quantum Hall effect in bilayer graphene.

###### pacs:

Graphene, a ‘Dirac material’, exhibits such novel phenomena as an anomalous integer quantum Hall effect, weak anti-localization, and ‘zitterbewegung’, which stem from its low energy excitations being massless chiral Dirac fermions. (1); (2) In addition, it holds great potential for applications due to its high mobility and the experimental ability to tune its properties via gating. (2) Early experiments aiming to characterize the electronic properties of graphene used quantum Hall measurements to probe its underlying structure and unequivocally demonstrate the relativistic nature of the charge carriers (3); (4). In these studies, the Hall conductivity displayed plateaus at atypical values of , where is an integer. The origin of these plateaus is well known to be a consequence of the electrons in graphene being governed by a relativistic Dirac Hamiltonian, together with the inherent valley and spin degeneracy. The chiral symmetry of the Dirac Hamiltonian generates a particle-hole symmetric Landau level (LL) spectrum – every positive energy LL has a conjugate negative energy LL that together are responsible for the positive and negative conductivity plateaus. The Dirac Hamiltonian also supports an additional zero energy LL, and leads to the ‘half-step’ offset in the first conductivity plateaus. Finally, the step size between each plateau can be attributed to the presence of four degenerate Dirac cones labelled by different spin/valley indices. Further information about the velocity of the Dirac cone, , can also be obtained by measuring the energy gap between the LL, which scales as . Thus, information about LL spectrum and Hall conductivity of a system can provide direct evidence for the existence of Dirac-like quasiparticles and can be used to probe the degeneracy and velocity of the Dirac cones.

Inspired by studies of superlattices (SLs) in semiconductor systems (5), such SLs have been proposed as a route to band structure engineering in graphene, paving the way to further new physics and applications. Recent studies of one dimensional (1D) superlattices (SLs) in monolayer graphene (MLG) have demonstrated their remarkable ability to not only regulate the Fermi velocity but also generate new Dirac cones. (6); (7); (8); (9); (10); (11); (12); (13); (14); (15); (16); (17); (18); (19); (20) In seminal work, Park et al., (6) showed that a one dimensional (1D) SL in MLG could lead to a strong Fermi velocity renormalization in the direction transverse to the SL modulation, leading to strong transport anisotropies. Later, it was shown that SLs can even lead to nonzero energy Dirac points (8) and multiple zero energy Dirac points near each valley. (9); (10) These effects stem from the chirality of the low energy quasiparticles in MLG. The new Dirac points generated by the SL lead to extra zero energy Landau levels for SLs in a weak orbital magnetic field. (9)

For bilayer graphene (BLG), a minimal model of the low energy physics has quadratic band touching points which can be gapped out by an electric field perpendicular to the layers.(21); (22) 1D SLs in BLG are also predicted to display remarkable restructuring of the band dispersion, leading to new Dirac points in the spectrum which have anisotropic velocities determined by the SL potential, and are perturbatively stable against interactions (unlike the quadratic band touching point). For bilayer systems, there are two possible types of electrostatic SLs: (i) a chemical potential modulation where the potential profiles are the same on both layers, or (ii) an interlayer bias modulation where the layer potentials are of opposite sign leading to a periodically modulated electric field.

The effect of a chemical potential SL in BLG is highly dependent on the the strength of the modulation — the band structure goes through a series of stages as the strength of the modulation is swept from weak to strong.(16); (23) Starting from a weak potential, the quadratic band touching point splits into two anisotropic Dirac cones. This happens while preserving the net “pseudospin winding number”.(24) These Dirac points continually move out towards the mini-Brillouin zone (MBZ) before becoming gapped at a certain critical SL amplitude. (16); (23) For even stronger SLs potentials, the widening of the band gap reverts and eventually closes by forming a series of new zero-energy Dirac points. The protection of the Dirac points has been shown to be directly related to the chirality of the low energy BLG quasiparticles.(16); (23)

In the case of symmetric 1D electric field SLs in BLG, four anisotropic Dirac cones are generated, two at zero energy and two at finite energy. The existence of these Dirac points can be seen (23) to arise from a series of coupled ‘topological’ modes (25); (15); (23); (26) that emerge along the zero-line where the parity of interlayer bias reverses. Interestingly, these Dirac points are exceptionally robust to the size of the bias modulation. In the limit of a long period SL, these Dirac points evolve into independent topologically confined one-dimensional chiral modes.

In this paper, we investigate the electronic properties of both bilayer and monolayer graphene 1D SLs in the presence of a uniform magnetic field perpendicular to the plane. For the bilayer, we consider two types of SLs: a chemical potential modulation and an electric field modulation. Similar to early magnetic field studies of graphene and graphene SLs (3); (27); (9), careful analysis of the Landau levels spectrum and Hall conductivity proffers many potential experimental signatures of the underlying Dirac cones generated by the SL. Special attention is also made to electric field modulations where the period length is large and decoupled ‘topological’ modes emerge. Although our SL results are most relevant to graphene with chiral low energy excitations, the effects explored in this paper may also find interesting counterparts in low density two dimensional electron gas in semiconductor heterostructures, in which a honeycomb pattern has been successfully imprinted, with a possibility of other kinds of SL patterns imprinted in the future.(28)

We obtain the following key results for 1D SLs in MLG and BLG in a perpendicular (orbital) magnetic field. (1) We show that transport in a weak magnetic field can not only be used to probe for the presence of extra Dirac points, but may also be used, in clean systems, to obtain further information about their velocities via the energy and degeneracy of higher Landau levels (LLs). (2) We show that a moderate magnetic field in monolayer graphene leads to a dramatic reversal of the transport anisotropy generated by the SL potential alone, an effect arising from the SL induced dispersion of the zero energy LL. This field tunable transport anisotropy may find useful applications in monolayer graphene SLs. For example, one can use the change of the anisotropy as an on/off switch and even perform bit or gate operations with these SLs. However, this field tunable transport anisotropy is found to be absent in BLG. (3) We consider the Landau levels in BLG in the presence of a uniform bias, and a kink in the bias which leads to 1D topologically bound states, and the coupling between such modes in an array of kinks. We also discuss the real space structure of the LL wave functions and the local density of states which is expected to be relevant to scanning tunneling spectroscopy (STS) studies such as those described in Refs. (29); (30). (4) Finally, we consider possible implications of these results for valley filtering and breakdown of the quantum Hall effect.

For realistic SL period, with being the nearest neighbor lattice constant, the field strengths corresponding to weak and intermediate magnetic field are about 0.1T and 8T, respectively. Also, for this SL period, the minimum SL potential required to generate three Dirac points in one valley is about 200meV. All these values are experimentally achievable in order to see interesting physics described in points (1) and (2) above.

The organization of this paper is as follows: Section I concerns 1D superlattices in monolayer graphene in a magnetic field. We begin this section with a brief description of the band dispersion of 1D superlattices in the absence of a magnetic field. In the next part, Section I.2, we introduce the low energy effective Hamiltonian and then identify the following three magnetic field regimes that lead to very different electronic dispersions: (i) weak, (ii) intermediate and (iii) high fields. We present and the Landau level dispersions of these three regimes in Section I.2.1 - I.2.3, respectively. We then present the diagonal and Hall conductivity in Section I.3. We find that the diagonal conductivity shows strong anisotropy reversal and Hall conductivity no longer exhibits plateaus, when magnetic field is tuned from weak to intermediate strength.

Section II of this paper examines 1D superlattices in bilayer graphene and begins with description of the low energy theory that describes both chemical potential and electric field SLs. We discuss the LL dispersion and transport for 1D chemical potential and electric field SLs in Sections II.2.1 and II.2.2, respectively.

Then in the Section III, we consider the real space picture of interlayer modulations with dilute kinks and antikinks that form decoupled soliton modes. We first provide an overview of the LL spectrum of uniformly biased BLG in Section III.1 and then compare it to the spectrum generated when there is a single isolated antikink in the interlayer bias in the beginning of Section III.2. Also in this section, we present the tunneling current and provide key signatures that could identify the presence of the kink in STS measurements and suggest how the occurrence of disorder induced kinks may contribute to the breakdown of the quantum Hall effect. Following the discussion of a single isolated interlayer bias kink, we discuss a periodic array of kink and antikinks in Section III.3. We argue that for , it is only necessary to consider a single kink-antikink pair and observe that regardless of the magnetic field strength, the low energy physics is governed by the soliton modes which remain robust. Employing a simple low energy effective theory to describe the soliton modes in a magnetic field, we gain further insights into the low energy physics of the SL and discuss why the zero energy LLs are robust from this perspective. Finally, in Section III.4 we consider the coupling between the soliton states generated by a kink-antikink pair and how this leads to a single 1D one-way conducting band with definite valley index in each kink separately.

## I Monolayer graphene superlattices

### i.1 Band structure of monolayer graphene superlattices

It is useful to briefly review the effect of the SL potential (9); (10); (13); (14) at zero magnetic field on the low energy spectrum near before we present the SL spectrum in a nonzero magnetic field (see Fig. 1). For , where is the SL strength, the spectrum consists of a single Dirac point at , with an isotropic velocity at low energy. The first effect of turning on a nonzero is to make this Dirac cone anisotropic, with , but . For , with as the SL period, one finds . Further increasing , with , leads to three zero energy Dirac points: one of these continues to be located at , while two new Dirac points emerge which are symmetrically split off from the point along the -direction. All three Dirac cones have anisotropic velocities.(13) It is clearly convenient to measure the strength of the SL potential in terms of a dimensionless parameter , defined via , so that tuning with leads to the different interesting spectra discussed above. A further increase in leads to a more complex spectrum with even more Dirac points.

### i.2 Landau Level spectrum

Ignoring spin, the low energy Hamiltonian of pristine MLG is given by a matrix at each valley, , where pseudospin labels the two sublattices, while the two (decoupled) valleys at are labelled by . Here, is the isotropic Fermi velocity, with and eV being the nearest neighbor carbon-carbon distance and transfer integral respectively. (We set for convenience.) In a uniform perpendicular magnetic field, ; for , and in the Landau gauge, the vector potential .

Diagonalizing , we find energies , where , with . For (i.e., at valley ), the eigenfunctions are given by

(1) |

where and are the system length and electron momentum deviation from , both along the -direction, while for ,

(2) |

Here, is the n-th eigenstate of a (shifted) 1D harmonic oscillator,

(3) | |||||

centered at , and are Hermite polynomials. For (i.e., at ), the eigenfunctions are given by . The full low energy LLs of MLG are thus .

We now turn to the effect of a periodic 1D chemical potential modulation , with period , on these Landau levels. The set of eigenfunctions , with , form a convenient basis to study the SL Hamiltonian in a magnetic field. (This basis choice is different from the one used by Park, et al, and allows us to numerically access a wide range of magnetic fields.) Due to momentum conservation along the -direction, the SL Hamiltonian is diagonal in . Further, for , intervalley scattering is strongly suppressed. We will therefore assume that the two valleys stay completely decoupled. (We focus below on valley with ; we expect identical physics around valley .) With this approximation, the only effect of the SL potential is, thus, to induce Landau level mixing.

To proceed, we need to choose a concrete form for the SL potential. For simplicity, we set , although our results can be easily generalized to other SL potentials which are smooth on the atomic scale by including multiple Fourier components. For numerical computations, the following integral identity (31) for harmonic oscillator states proves useful:

(4) | |||||

where , and is the generalized Laguerre polynomial. Using this identity, we simplify and numerically compute the matrix elements of the full Hamiltonian, retaining up to 3000 Landau levels, and diagonalize this to obtain the spectrum of the 1D SL in a magnetic field.

In order to study the effect of the magnetic field on the 1D SL in graphene, with , it is useful to consider three regimes for the magnetic field.

(i) Weak field: This regime corresponds to having , where the Landau level spacing is much smaller than the SL amplitude, so that . In this regime, the magnetic field may be viewed as effectively ‘probing’ the zero field SL excitations.

(ii) Intermediate field: In this regime, , which means , so that the SL potential and the magnetic field have to be treated on equal footing.

(iii) Strong field: Here, or, equivalently, . In this regime, the SL potential only weakly perturbs the Landau levels of pristine graphene.

Fig. 2 shows the spectrum of the graphene SL in different field regimes for SL strengths (or ) and (or ). This allows us to contrast the behaviour of the spectrum of the SL in a magnetic field without or with extra Dirac points being present at zero field, and to explore consequences for quantum Hall physics and transport.

#### Weak field regime

When the magnetic field is weak, (top panels in Fig. 2), we find that the energy spectrum barely depends on the value of , or equivalently, . This is due to the fact that when magnetic length is larger than the SL period , the matrix elements of the Hamiltonian do not depend on the center of the LL wavefunctions, which yields flat bands. Equivalently, in this regime, the magnetic field may be viewed as effectively ‘probing’ the structure of the zero field SL dispersion leading to Landau levels which depend on the nature of the Dirac spectrum at low energy.

For , the low energy spectrum of the SL contains a single anisotropic Dirac point at zero energy. For an anisotropic Dirac cone described by an effective Hamiltonian , the LLs are given by . Since the SL renormalizes , but leaves , the Landau levels at weak field resemble those of pristine graphene, but with a renormalized lower effective velocity .

For , the low energy spectrum of the SL contains three anisotropic Dirac points at zero energy, so that the zero energy Landau level has three times the degeneracy of the case with . Further, the Dirac cone centred at has a slightly different average velocity compared with the two cones which are symmetrically split off from along . This degeneracy breaking results in the Landau levels at nonzero energy becoming weakly split, as is most clearly seen for the first two excited Landau levels (at positive or negative energy, i.e., with ). We have numerically determined and for each of the three Dirac points and found good agreement between the energy levels obtained on this basis of having Dirac fermions with two different average velocities, and that obtained directly numerically.

At higher energies, for or for , the spectrum begins to deviate from this simple behavior expected for a linear Dirac spectrum. This deviation results from curvature in the dispersion, which appears upon going beyond the linearized approximation.

#### Intermediate field regime

At intermediate fields, for , the spectrum at low energy is most simply understood as arising from the SL potential inducing a strong dispersion to the Landau levels. In simple terms, if we assume that the state labelled by momentum , or equivalently position , have an energy which is modulated by the SL potential, we expect a periodic modulation of this energy with period and amplitude proportional to the SL amplitude . The behaviour of the low energy Landau levels, , as seen from the lower panels in Fig.2, is consistent with this scenario, with the modulation following the form of the SL potential and the modulation for being roughly thrice as strong as the modulation for . We can also see that the low energy Landau levels when overlap with each other. This will have nontrivial effect on the dc conductivity, as shown in the following subsection. For higher energy Landau levels, the energy spectrum still has a periodic modulation but no longer retains the simple form of cosine function. This is due to the fact that as the energy gets higher, the distribution of Landau levels becomes more dense and the energy difference between two adjacent levels is now comparable to the matrix element of SL potential. Therefore, a simple first order perturbation correction is not enough to account for the dispersion and second order perturbation from adjacent levels must be taken in account, which causes the Landau level loses the simple cosine form.

#### High field regime

For very strong magnetic field, the Landau level structure of pristine graphene is recovered. Here, only one zero energy level exists at the Dirac point, and other energy levels follow the square root relation. This is simply because in such a strong magnetic field, the SL is just a perturbation and can only give rise to a small modulation of the LLs following our argument at intermediate field. From a perturbative point of view, the energy corrections up to first order to the LL energies are given by

(5) |

which, upon using the identity (4), gives a sinusoidal dependence on the center position of LL wavefunctions. Thus, even in a strong magnetic field, the energy spectrum is not dispersionless but has a spatial modulation following the SL; however the ratio of the amplitude of this modulation to the Landau level spacing, , is extremely small in the high field regime. This dispersion, though small, can give rise to interesting magnetoresistance oscillation known as Weiss oscillation, on top of the usual Shubnikov-de Hass oscillation. (32) It was shown that, compared to two-dimensional electron gas with parabolic dispersion relation, Weiss oscillation in graphene SL is more pronounced and is more robust against temperature damping in small field region. This is a consequence of the different Fermi velocities of Dirac and normal electrons at same chemical potential. (32) Similarly, Weiss oscillation in bilayer graphene chemical potential SLs has also been discussed. (33)

### i.3 Diagonal and Hall Conductivity in MLG SLs

Once we have the eigenvalues and eigenfunctions for the superlattice in a perpendicular magnetic field, both ac and dc conductivities can be calculated directly by Kubo formula, (34)

(6) | |||||

Here, we have set as the Landau level broadening and measure the conductivities in the unit of , and are the -th eigenvalue and the corresponding eigenstate of the system which can be expanded in the basis of , with and being the LL wavefunctions for pristine graphene. is the Fermi-Dirac distribution, and is the velocity operator in -direction and , where is the Pauli matrix. Using Eq. (1), (2), we can get the expression for the matrix elements of velocity operators from

(7) | |||||

and

(8) | |||||

by expanding into . Note that is always true for any state. This is because in the direction, the wavefunction is always localized.

Results for dc diagonal conductivities as a function of chemical potential are shown in Fig. 3. This can be done by setting the frequency to zero in Eq. (6), and only the real part of the conductivity tensor is nonzero. In weak magnetic field, the conductivities show strong anisotropy, with larger than , which is a direct consequence of the Fermi velocity renormalization in the absence of magnetic field (see Fig. 3 (a) and (c)). Besides, since and due to the flat band structure, the major contribution to the diagonal conductivities comes from off-diagonal matrix elements, with . Since the matrix elements of are always larger than those of due to the anisotropic dispersion, this gives rise to the anisotropy in the weak field. In intermediate magnetic field, conductivities still show anisotropy, but with significantly larger than (see Fig. 3 (b) and (d)). This is because has acquired diagonal matrix element, since the energy spectrum is dispersive, while still lacks this contribution. Notice the positions of the conductivity peaks of exactly correspond to the minimum and maximum of the energy band, where the density of states is the largest. For , however, the conductivity is minimum at the band edge, since the average of the velocity operator, , is zero at these energies. Therefore, the intra-LL contribution to is the smallest at the band edge and thus leads to dips in it. For weak SL potential, can drop to zero when there is no overlapping LLs, while in a strong SL, always show dispersive transport property. In strong magnetic field (not shown here), where the Landau levels structure of pristine graphene is recovered, the conductivities become isotropic (see, for example, Ref. (35)). In this case, the SL is merely a perturbation to the magnetic field and thus should have minor effect on determining the magnetotransport properties.

The dc Hall conductivity is shown in Fig. 4. For weak magnetic field (Fig. 4 (a) and (c)), the Hall conductivity shows well-defined plateaus, as a consequence of nearly flat energy bands. The values of Hall conductivity around Dirac points, taking into account of spin and valley degeneracies, are in weak SL () and in strong SL (). This result resembles the anomalous half integer quantum Hall effect in pristine graphene and the Hall conductivity triples due to the existence of three Dirac points in a strong SL. Moving away from the Dirac point, we can observe quantum Hall plateaus with higher conductivities, and the value increases by 1 each time the chemical potential crosses an LL. For intermediate magnetic field (Fig. 4 (b) and (d)), there is no longer well defined plateaus due to the dispersive energy spectrum. However, for weak SL, the LLs are not overlapped with each other. If chemical potential falls between two LLs, a small plateau can still show up, with the value expected from Dirac physics. When the magnetic field becomes strong enough as the LL structure for pristine graphene is restored, Hall conductivity will show anomalous half integer quantum Hall plateaus.

Fig. 5 shows the ac conductivities of graphene SLs in an intermediate magnetic field. For weak and strong magnetic fields, the results resemble those of pristine graphene,(35) since in both cases the LLs are nearly flat and the real part of the conductivities show strong peaks when photon frequencies exactly correspond to the energy differences between two LLs. In an intermediate magnetic field, the result is complicated by the dispersion of LLs. At low frequencies, there can be optical transitions in a range of photon energies, and the real part of diagonal conductivities is maximum at the band edge where the DOS is also maximum. At high frequencies, the LLs become less dispersive and peaks will show up. These results can be linked with graphene’s unusual magneto-optical properties, for example, giant Faraday rotation.(35); (36) While the anisotropy in the diagonal conductivities can lead to anisotropic rotation angles for incident waves with different polarization plane, this effect is actually quite small and hard to observe experimentally.

## Ii Bilayer Graphene Superlattices

### ii.1 Band structure of bilayer graphene superlattices

Before presenting the results for BLG SLs in a perpendicular magnetic field, we briefly review the band structures of BLG SLs in the absence of a magnetic field. In BLG, there are two distinct types of SLs, namely, chemical potential and electric field SLs. Other SLs can be thought as a superposition of these two fundamental ones, and can be similarly studied.

In a chemical potential SL, the SL potentials are exactly the same on each layer. As has been shown by us in previous work, (15) the electron and hole state are completely decoupled along the SL direction. Therefore, the SL correction to the energy spectrum pushes electron (hole) states down (up) and causes the two quadratic bands cross each other inside the mini Brillouin zone (MBZ). However, in other directions, this decoupling is absent and level anti-crossing gives rise to a gapped spectrum which increasing linearly as we increase the angle with respect to the SL modulation direction. This induces the appearance of two Dirac cones. With a further increase in the SL strength, the two band crossing points move outside MBZ and the spectrum gets gapped out (see top figures of Fig. 6).

In an electric field SL, the situation is drastically different. Now, the SL potentials on the two layers are exactly opposite in sign, which can effectively be viewed as a periodic arrangement of potential kinks and anti-kinks. Earlier results have indicated that two topological 1D modes will be confined to each kink or anti-kink.(23); (25); (26) These topological modes propagate perpendicular to the SL direction and have opposite velocities at a kink and an anti-kink. When one has such modes arranged periodically, they will couple each other and zero energy and finite energy Dirac cones will appear, as explained by an effective “wire” model in LABEL:Killi1 (see bottom figures of Fig. 6).

### ii.2 Magnetic properties of bilayer graphene superlattices

From the previous section, we can see that the phenomenon resulting from the existence of Dirac points is most manifest when magnetic field is weak compared to the SL strength. A very special feature of bilayer graphene SL is that anisotropic Dirac points can be generated in the energy spectrum, where originally only quadratic band touching points is present. Our numerical results also indicate that in intermediate magnetic field, the LLs are only slightly dispersive, compared to monolayer graphene. Therefore, in this section, we will focus on the properties of bilayer graphene SLs in weak magnetic fields using a two-band effective Hamiltonian. We hope that our results can provide means to the determination and characterization of these anisotropic Dirac points.

The low energy physics of bilayer graphene near the Brillouin zone corners at is, in the presence of a perpendicular magnetic field, described by the following two-band effective Hamiltonian,

(9) |

with corresponding to as in the case of monolayer graphene. Similar to the single layer case, we have replaced the momentum operator with its canonical counterpart to take into account of the magnetic field effect. In the following, we will work with the same gauge choice, . For , the eigenvalues and the corresponding eigenvectors of the above Hamiltonian are

(10) |

with . In addition, there are two zero energy solutions,

(11) |

For , the corresponding eigenvectors are given by . The full low energy LL wavefunctions thus take the form . and these serve as a good basis to study the magnetic field effect of bilayer graphene SLs.

The SL can be modelled by assuming spatially varying potential profiles on each layer,

(12) |

where different combination of and gives rise to different types of SL, e.g., for chemical potential SL and for electric field SL. In following, we will only consider 1D SLs, where for proper choice of SL strength, various anisotropic Dirac points can be generated in the spectrum.

#### Chemical potential superlattice

For 1D chemical potential SLs, we will take . In previous work, (15) we have shown that for 1D chemical potential SL, two anisotropic Dirac points will derive from the original quadratic band touching points as long as the SL strength is smaller than a critical value. This results from the chiral symmetry protected level crossing in the SL direction and level anti-crossing in other directions. Once SL strength increases beyond this critical value, Bragg reflection will open up a gap in the energy spectrum on the mini Brillouin zone boundary. This behavior has also been found by Tan et al. (16)

Results for energy spectrum of chemical potential SL in a weak magnetic field are shown in Fig. 7. In the left panel, we have chosen , where eV is the intralayer nearest neighbor transfer integral. For this value of SL potential, our previous calculation has shown that there will be two Dirac points generated in the mini Brillouin zone. Here, similar to the single layer case, the energy spectrum shows flat band structure. At the Dirac point, there are two nearly degenerate zero energy levels. Right above and below the Dirac point, we can also observe two sets of doubly degenerate energy levels. These are the LLs of each anisotropic Dirac cone. Moving further away from the Dirac point, no such degenerate levels are present. This is because these Dirac points only have linear dispersion at rather low energy. In the right panel, , energy gap opens up and there exist nonzero energy levels near the Dirac point. Energy bands in intermediate field regime (not shown here) are only slightly dispersive, where the amplitudes of dispersion are extremely small compared to the LL spacings.

Fig. 8 shows the dc diagonal conductivity of chemical potential SLs. Similar to the single layer case, we can observe the anisotropy in the conductivity. However, there is no anisotropy reversal as magnetic field strength is tuned. In weak field, the transport anisotropy is determined by the anisotropic Dirac cones. From Ref. (15), for emergent Dirac cones. Therefore, the conductivity in the direction, , should be larger than in a weak magnetic field. For intermediate magnetic field, due to the dispersion of the energy bands, the average velocity in the direction is not zero, . On the other hand, is always equal to zero. This means that will acquire intra-LL contributions, while is mainly determined by inter-LL contributions and is thus small compared to .

We can further demonstrate how the flat energy levels evolve as the strength of the SL potential varies, as shown in Fig. 9. Here, we have fixed the magnetic field strength by setting and also fix the center position of the wavefunctions to be at 0. As the SL potential tends to zero, where the problem reduces to pristine bilayer graphene in a perpendicular magnetic field, we can see the energies of LLs will follow the pattern of Eq. (10). As SL potential increases, the physics gradually becomes dominated by anisotropic Dirac cones, which can be seen from the appearance of doubly degenerate levels at nonzero energies. This crossover from non-relativistic to relativistic behavior can be qualitatively understood by considering the competition between the characteristic energy scales in these two regimes. In pristine bilayer graphene, the low energy excitation are massive electron with an effective mass . In a magnetic field , the characteristic energy scale is the cyclotron frequency, . On the other hand, the anisotropic Dirac points generated by SL have anisotropic Fermi velocities and , where . (15) Therefore, the characteristic energy scale associated with these Dirac points is . We expect to see a smooth crossover as these two energy scales are comparable to each other. A rough estimate shows that the crossover should occur around , which is quite close to the value that we have observed.

Further increasing the SL potential, the doubly degenerate zero energy levels become gapped and all levels are pushed away from Dirac point. Surprisingly, at rather strong SL potential, , zero energy LLs appear again, and all of the higher energy levels become doubly degenerate. This phenomenon can be understood from the result of Tan et al. (16) As it has been shown, for a chemical potential SL, when SL potential is strong enough, anisotropic Dirac cones will show up again in the energy spectrum, which naturally leads to the zero energy LL in the presence of a magnetic field. According to Ref. (16), for certain values of SL potential, there can be four Dirac points in the spectrum, which is not present in our calculation. Of course, for strong SL potential, the two-band description of bilayer graphene is no longer valid. However, we have reproduced this result in a four-band effective Hamiltonian approach. In strong chemical potential SL, there can exist four-fold degenerate zero-energy LLs, and the degeneracy reduces to two upon further increasing the SL strength.

#### Electric field superlattice

For 1D electric field SLs, we will take . Following the same procedure as in the previous section, we have calculated the energy spectrum and the dc conductivities for an electric field SL in a magnetic field.

In a 1D symmetric electric field SL, there are always two zero energy Dirac points present in the spectrum, which results from the coupling of 1D zero modes at kink/antikink of the SL potential profile. Earlier, we have developed an effective “wire” model to describe these Dirac points and demonstrated that they are topological stable if a generalized “inversion” symmetry, (i.e., flipping the electric field and shifting by ), is preserved.(15) This implies when a weak magnetic field is applied, doubly degenerate zero energy levels should appear at the Dirac point. Indeed, as it can be seen from Fig. 10 and 11, these two levels always stay at zero energy, independent of the SL potential. Moreover, for the SL potential we choose, , nearly degenerate levels appear even at nonzero energies (Fig. 10). These correspond to the LLs derived from anisotropic Dirac points, up to . From Fig. 11, it is more clear that at strong SL potential, physics is strongly dominated by the Dirac points, where higher energy levels become doubly degenerate and resemble the higher LLs of the Dirac cones. When SL potential is weak, equally spaced LLs are recovered, as in the chemical potential SL, which also indicates a non relativistic to relativistic crossover at certain SL strength. Remarkably, and different from the chemical potential SL case, the relativistic behavior survives to higher energies as SL potential increases, which means the linear approximation description of Dirac cones works in a larger energy range. This is consistent with earlier result. (23) Consider a potential profile with a kink of bias reversal. When the bias is large, the conduction and valence bands are widely separated and a pair of linearly dispersed zero modes will traverse the gap. Now, when these kinks and antikinks are arranged periodically, Dirac cones will result from the coupling of these modes, and the Fermi velocity along SL direction inherits its value from the freestanding zero mode. (15) Therefore, as the SL potential increases, the energy range where the Dirac cone approximation is valid also increases, which leads to the robust relativistic physics at large SL potential.

## Iii Real space picture for interlayer bias modulations in bilayer graphene

As suggested in an early study by the present authors (15), it is extremely useful to view an interlayer bias modulation as establishing a series of ‘kink’ and ‘antikink’ modes along the zero-lines where the interlayer bias reverses polarity. Within the effective low-energy description of BLG, the interlayer bias can be regarded as a mass generating term and the zero-lines represent boundaries where the sign of the mass changes. As consequence, along a single isolated kink (antikink), a pair of two right (left) moving ‘topological’ modes at and corresponding two left (right) moving modes at the other (25); (37); (38); (15). The typical length scale over which these kink-antikink modes are confined, , varies from for to lattice sites for , where is the length the primitive lattice vector. For the remainder of this section we set .

Accordingly, an interlayer bias modulation can be thought of as a periodic array of coupled kink-antikink solitons. With this viewpoint, and using the symmetry properties of the kink and antikink modes, Killi et al. (15) derived an effective low-energy hamiltonian describing how anisotropic Dirac cones precipitate precisely at the band crossing point between the kink and antikink modes when the soliton modes overlap and couple. With this same perspective, further intuition into the quantum Hall effect in the SL system can be obtained by studying the properties of 1D kink states in a magnetic field.

Aside from its relevance to SLs, a study into the properties of the kink/antikink states also provides valuable insights into disorder BLG samples in the quantum Hall regime. Various sources of disorder are known to generate a random electrostatic potential landscape. Throughout these samples, 1D kink-states are expected to percolate along precipitous fluctuations that reverse the parity of interlayer bias. If the percolation networks are well extended, these kink states would contribute to the breakdown of the quantum Hall effect (29). Recent STS measurements of BLG samples in the quantum Hall regime also indicate that even when the sample is uniformly biased by external gates, the disorder is still strong enough to locally reverse the polarity of the interlayer bias (30); (39). These recent measurements provide an additional motivation to study the properties of kink states in a magnetic field and, specifically, to examine how they modify the tunneling current.

In this section, we consider interlayer bias modulations in quantum Hall regime from the perspective of the kink-antikink modes. We study the dilute limit where the kink and anti-kink states are well separated and confined to the zero-lines (i.e. ). Within this regime, the low energy states of the system are completely dominated by midgap soliton modes. We start by first reviewing the effects of a magnetic field on uniformly biased BLG in section III.1 and then non-uniformly biased BLG with a single isolated anti-kink in the potential profile in section III.2 . In addition to studying the dispersion relation, we examine how the tunneling current is modified in the proximately of a kink. Next, we consider an array of decoupled kink and anti-kink modes in Section III.3 . For magnetic fields such that , it is only necessary to understand the dispersion relation of an isolated kink/anti-kink pair as this is sufficient to describe the entire band structure. This allows us to consider a zigzag ribbon of width with one period of modulation (i.e. a pair of kink anti-kink states) and implement a Peierl’s substitution in the full tight-binding Hamiltonian. For weaker magnetic fields such that , it becomes necessary to include a larger number of kink/anti-kink states to describe the high energy features of this system. Despite this, further insights into the low energy modes can be made by combining the knowledge obtained from the study a single kink/anti-kink pair and by using a simple low energy effective theory (see Ref. (40)). Aside from this, in section III.4 we demonstrate that a coupled pair of kink/anti-kink states opens asymmetric bandgaps in each valley that can be precisely controlled using a magnetic field. This effect may provide a viable route to developing a switchable valley filter.

### iii.1 Uniform Bias

Although it has been discussed extensively in the literature, we review and emphasize a few important points about the LL structure in uniformly biased BLG. Details about the LL structure in biased BLG can be found in Ref. (21); (22); (41); (42); (43); (44); (45); (46); here we only briefly summarize some of the features as they pertain to our current discussion.

In the Landau gauge that preserves translation symmetry along the zigzag direction, inspection of Fig. 12 (a) shows that there are two distinct energy level spectrums in the dispersion relation, one about each valley. Each energy level consists of two well-defined regions: a flat, non-dispersive region consisting of bulk LLs and dispersive edge states (44); (45) (note, the edge state dispersion for the and are unique in that they contained non-dispersive regions, as seen in Fig. 12 (a,c)). The wavefunctions of the non-dispersive LL states are well confined within the bulk and, starting from states on the far left side of the level and ending on the right side, move from one edge of the sample to the other.

A comparison of the spectrum about the K-points clearly shows that the interlayer bias breaks the valley degeneracy of the LL spectrums (22); (43). Although the valley degeneracy is broken, the states in each valley are connected by a symmetry that relates the low energy Hamiltonian of opposite -point. If we define the operator to correspond to the interchanging of sublattices combined with layer exchange (i.e. and ), it is a simple matter to show that . Thus, upon interchanging the sites labelling, the quasiparticles in each valley are governed by the same low energy Hamiltonian but with the effective relative parity of the interlayer bias reversed, leading to the observed valley degeneracy breaking.

Inspection of the band structure shown in Fig. 12 (c) reveals that in the presence of an interlayer bias, the zero energy and LLs are shifted to finite energy and their degeneracy is lifted. Note, the energy shift of these two these states is positive in one valley an negative in the other. A second crucial observation about the and Landau states (and to a lesser extent the higher energy states) is that the wavefunctions are strongly localized to one of the two layers determined by its valley index. Such layer polarization is directly observed in STS measurements and information about the LL spectrum, valley index and local interlayer bias can be obtained (30). Below, we discuss further how the presence of a kink and anti-kink affects the tunneling current and give rise to clear signatures observable in STS measurements.

### iii.2 Single Kink or Antikink in the Bias

In the absence of a magnetic field, an isolated kink in the interlayer bias generates two unidirectional dispersing subgap bands in one valley and two oppositely dispersing bands in the other valley related by time reversal. An anti-kink generates similar low energy bands, although the velocity of these modes is reversed in each valley (see Fig. 12 and Ref. (25) for details).

When a large magnetic field is introduced, the energy spectrum resembles that of the LL structure of uniformly bias BLG, but with further distinctive features. Consider the energy spectrum about a single valley. Comparison of Fig. 12 (c) and (d) shows that instead of each energy level having a single non-dispersive bulk LL, each LL breaks into two flat regions shifted in energy, which are connected through a dispersive mode. States in these two flat regions correspond trivially to the bulk Landau states that would have been generated by having a uniform interlayer bias of the same parity as the local bias. Hence, these two regions are composed entirely of Landau states that are well localized in the bulk and respond only the local interlayer bias. As in the previous case, far way from the -point the bulk states evolve into edge states, and are again sensitive only to the local interlayer bias at that edge.

However, it is for the and LLs that the presence of the kink states are directly observed. Starting from a bulk state in the () Landau level on one side of the interface and tracking the states to the other side of interface (moving left to right in momentum space), we observe the following sequence. The bulk () LL states begin to continuously evolve into dispersive kink states as the guiding center approaches the zero-line and then emerge from the interface on the other side as bulk states in the () LL.

Although other dispersive states connecting the high energy bulk LL on either side of the interface are also present, they are not of topological origin. These states are simply LL states that becomes dispersive as their wavefunction begins to overlap with the region of spatially varying potential. Hence, the properties of these states are inherently sensitive to the the magnetic field and the details of potential profile. In contrast, the kink-states are of topological origin and so persist even in the absence of a magnetic field. Due the strong confined of these states, the properties of these states are remarkably robust and insensitive to the magnetic field (see Zarenia et al. (47) for more details).

The tunneling current of the sites in the vicinity of a kink is shown in Fig. 13. It was computed by allowing electrons to tunnel within an energy window between the chemical potential and the sample-to-tip bias . Away from the interface, sharp LL plateaus easily distinguish the LL spectrum. The marked suppression of the tunneling current on the other side of the interface where the parity of the bias is reversed is indicative of the strong layer polarizability discussed previously. Only when the lies above the LL is there any weight on the sublattice in this region, consistent with the low energy theory.

When lies between the chemical potential and the LL, there is notable enhancement of the tunneling current close to the interface. A cut along shown in Fig. 13 (b) shows signatures that this enhancement is due to the presence of strongly confined kink states. Namely, the confinement modes to the zero-line manifests as a transfer of weight from the bulk to the interface and results in the combined downturn in the bulk and sudden enhancement in the tunneling current close to the interface. Evidence of the strong localized states can also be seen when lies just above the LL and lies just below. Here, the tunneling into the kink states is suppressed by Pauli blocking, leading to a sharp reduction in the tunneling current close to the interface. In contrast, at higher energies the reduction in the tunneling current is a very gradual rolloff across the interface, showing no indication of localized interface modes.

Evidence demonstrating the robustness of these modes in the presence of magnetic field suggest these modes may play an important role in establishing percolation networks responsible for the break down of the quantum Hall effect.(29)

### iii.3 Array of Kinks and Antikinks

We begin by first considering the dilute limit where the kink and anti-kink soliton modes are well separated and decouple. In the the regime where , the LL states are well localized and respond only to their local environment. Hence, the dispersion relation is expected to consist of a series of bulk regions – identical to those generated by a uniform bias with alternating bias parity – connected by intermediate kink/anti-kink states. This viewpoint is confirmed by inspecting the dispersion relation of a single kink/anti-kink pair (which, neglecting the spurious edge states, compose the unit cell) shown in Fig. 14 (b).

In the opposite regime where , the dispersion relation in Fig. 14 (c) indicates the bulk LLs become dispersive in general. Here, the wavefunctions are spread over a wide region encompassing the kink/anti-kink of the potential profile, and thereby altering their properties. In contrast, with respect to the low energy kink/antikink states the magnetic field again does little more than shift their position in momentum space; the mode velocity is robust and the wavefunctions are rigid. Only at higher energies, close to the LL energy levels are the dispersive kink states modified.

In sum, the effect of a increasing magnetic field at low energy is to shift the kink and anti-kink modes relative to each other (the anti-kink modes move right the kink modes move left in momentum space). Since the mode velocity of the both the kink and anti-kink are reversed in each valley, this shift pushing the band crossing points up in energy in one valley and down in the other. In this way, a magnetic field can be used to control the band crossing points between modes of adjacent wires.

Further insights into the robustness of the low energy modes present in the continuum model of the 1D interlayer bias SL studied above can also be made by employing a simple low energy theory inspired by the kink/anti-kink viewpoint. In the dilute limit and in the absence of a magnetic field each kink (anti-kink) generates an identical copy of the low energy modes about each valley. At each K-point in Fig. 14 (a), the two modes for each kink (antikink) can be linearized about zero energy (see Ref. (23) for more details). An effective low energy hamiltonian of the SL about one of the K-points can be written as

(13) |

where is the velocity of the modes about zero momentum, and () are operators that create electrons on wire , with momentum about momentum close to one of the two zero energy crossing points labelled . (Note, the low energy theory of the opposite valley is identical except that the velocity of the kink/anti-kink modes are reversed.) Representing a magnetic field by the Landau gauge field , electrons hopping along a given wire see an average constant gauge field given by . In a magnetic field, the low energy effective Hamiltonian becomes

(14) |

These results demonstrate that the full dispersion of the SL consists of multiple copies a single kink-antikink pair dispersion repeated every in momentum space, consistent with the numerically computation of the dispersions relations provided above. Notice, that the modes in a given wire crosses modes in the two neighbouring wires symmetrically about zero energy. If the magnetic field is weak, the kink modes of wire cross with the anti-kink modes of wire at an energy of . If the modes couple, level repulsion will cause the kink/antikink modes of the neighbouring wire to mix and split, but because the bandcrossing is symmetric about zero energy and each the each mode couples to both neighbouring modes equally (for a symmetric SL), a zero energy band remains intact.

### iii.4 Valley Filter

In studying a kink/anti-kink single pair in the coupled regime, we observe a remarkable effect not observed in a previous study (47) of the kink states that used a low energy model that did not capture the valley degeneracy lifting. As the wavefunctions of the kink and antikink modes begin to overlap, a bandgap opens at the band crossing points. Since the size of the bandgap is determined by the degree of overlap between the wavefunctions, it can be narrowed by either increasing the kinks-antikink separation or by reducing the interlayer bias, which causes the wavefunctions to spread further from the interface. Furthermore, as the magnetic field strength increases, the bandgap shifts to positive energy in one valley and negative energy in the other due to the valley asymmetry in the band crossing point. Remarkably, for a single coupled kink/anti-kink pair an energy window opens where only unidirectional modes are present along each kink and anti-kink (the direction of flow can be flipped by reversing the magnetic field). Moreover, all the conducting modes have the same valley index and so acts as a valley filter similar to those suggested along edge states(48) and line defects (49); (50) when the chemical potential is tuned to lie in the bandgap.

## Iv Summary and discussion

In this paper, we studied magnetic properties of MLG and BLG under various types of 1D external potentials. For MLG, we identified three regimes of magnetic field strengths that generate distinctive features in the LL dispersion and transport properties. Under a weak magnetic field, MLG SLs exhibit zero energy LLs whose degeneracies are identical to the numbers of Dirac points present in the spectrum. At higher energies but still within the linear range of the spectrum, differences between the Dirac cones cause the LL degeneracy to be lifted. Therefore, measurements that carefully determine the degeneracy of the LLs can be used to probe and characterize the underlying anisotropic Dirac cones. We further showed that the diagonal conductivities show strong anisotropy, with conductivity in the SL direction larger than that in the transverse direction. Interestingly, the anisotropy can be reversed for an intermediate magnetic field, where the LLs become dispersive. This field tunable anisotropy may find interesting device applications such as in switching or in resistive bits.

We then considered two types of SLs for BLG and again showed that analysis of the LL spectrum provides clear signatures of the underlying Dirac cones. However, the effects of a magnetic field were shown to be rather different in the bilayer SL systems than on the monolayer SL system. Although the diagonal conductivity of both SLs also have strong anisotropy like in the monolayer, there is no field tunable anisotropy reversal. Furthermore, regardless of the type of SL there are distinct crossovers from nonrelativistic physics to relativistic Dirac physics as the SL strength is tuned. In the case of a chemical potential modulation, the zero energy LL vanishes at a critical SL strength where the Dirac cones becomes gapped. As for the electric field SL, two zero energy LLs are always present irrespective of the strength of the magnetic field SL strength or magnetic field. This demonstrates a remarkable robustness of the Dirac points inherited from their topological origin.

We have also studied how a magnetic field affects BLG subject to a uniform bias and with bias reversing kinks in the potential profile. The motivation here was to gain a real space perspective of the LL wavefunctions in the quantum hall regime. Together with a low energy model, this study further established an intuition into the properties of the electric field SL. Moreover, it provided valuable insights into the topological kink states present in disorder BLG and suggests how that they may contribute to the breakdown of the quantum Hall effect. Finally, we proposed that a pair of coupled kink-antikink modes subject to a magnetic field could serve as a possible route to fabricate a switchable one-way valley filter.

###### Acknowledgements.

We acknowledge Gil Refael, Jeil Jung, and Jin-Luo Cheng for useful disscussions. Financial support was provided by the Early Research Award, Government of Ontario, NSERC of Canada and a University of Waterloo start-up grant.### References

- K.S. Novoselov, A.K. Geim, S.V. Morozov, D. Jiang, Y. Zhang, S.V. Dubonos, I.V. Gregorieva, and A.A. Firsov, Science 306, 666 (2004).
- A.H. Castro Neto, F. Guinea, N.M.R. Peres, K.S. Novoselov, and A.K. Geim, Rev. Mod. Phys. 81, 109 (2009).
- Y. Zhang, Y. W. Tan, H. L. Stormer, P. Kim, Nature 438, 201 (2005).
- Novoselov et al., Science 315, 5817 (2007).
- R. Tsu, Superlattice to Nanoelectronics (Elsevier, Oxford, UK, 2005).
- Cheol-Hwan Park, Li Yang, Young-Woo Son, Marvin L. Cohen, and Steven G. Louie, Nat. Phys. 4, 213 (2008).
- Cheol-Hwan Park, Young-Woo Son, Li Yang, Marvin L. Cohen, and Steven G. Louie, Nano Lett. 8, 2920 (2008).
- Cheol-Hwan Park, Li Yang, Young-Woo Son, Marvin L. Cohen, and Steven G. Louie, Phys. Rev. Lett. 101, 126804 (2008).
- Cheol-Hwan Park, Young-Woo Son, Li Yang, Marvin L. Cohen, and Steven G. Louie, Phys. Rev. Lett. 103, 046808 (2009).
- L. Brey and H.A. Fertig, Phys. Rev. Lett. 103, 046809 (2009).
- M. Barbier, F.M. Peeters, P. Vasilopoulos, J. Milton Pereira, Jr., Phys. Rev. B 77, 115446 (2008).
- Michaël Barbier, P. Vasilopoulos, F.M. Peeters, and J. Milton Pereira, Jr., Phys. Rev. B 79, 155402 (2009).
- M. Barbier, P. Vasilopoulos, and F.M. Peeters, Phys. Rev. B 81, 075438 (2010).
- M. Barbier, P. Vasilopoulos, and F.M. Peeters, Phil. Trans. R. Soc. A 368, 5499 (2010).
- Matthew Killi, Si Wu, and Arun Paramekanti, Phys. Rev. Lett. 107, 086801 (2011).
- Liang Z. Tan, Cheol-Hwan Park, and Steven G. Louie, Nano Lett. 11, 2596 (2011).
- D.P. Arovas, L. Brey, H.A. Fertig, Eun-Ah Kim, and K. Ziegler, New J. Phys. 12, 123020 (2011).
- P. Burset, A. Levy Yeyati, L. Brey, and H. A. Fertig, Phys. Rev. B 83, 195434 (2011).
- I. Pletkosić, M. Kralj, P. Pervan, R. Brako, J. Coraux, A.T. N’Diaye, C. Busse, and T. Michely, Phys. Rev. Lett. 102, 056808 (2009).
- S. Rusponi, M. Papagno, P. Moras, S. Vlaic, M. Etzkorn, P.M. Sheverdyaeva, D. Pacilé, H. Brune, and C. Carbone, Phys. Rev. Lett. 105, 246803 (2010).
- E. McCann and V. I. Falâko, Phys. Rev. Lett. 96, 086805 (2006).
- Edward McCann, Phys. Rev. B 74, 161403(R) (2006).
- M. Killi, T.-C. Wei, I. Affleck, and A. Paramekanti, Phys. Rev. Lett. 104, 216406 (2010).
- C.H. Park and N. Marzari, Phys. Rev. B 84, 205440 (2011).
- Ivar Martin, Ya. M. Blanter, and A.F. Morpurgo, Phys. Rev. Lett. 100, 036804 (2008).
- L.J.P. Xavier, J.M. Pereira, Jr., Andrey Chaves, G.A. Farias, and F.M. Peeters, App. Phys. Lett. 96, 212108 (2010).
- K. S. Novoselov et al., Nature Phys. 2, 177 (2006).
- A. Singha, M. Gibertini, B. Karmakar, S. Yuan, M. Polini, G. Vignale, M. I. Katsnelson, A. Pinczuk, L. N. Pfeiffer, K. W. West, and V. Pellegrini, Science 332, 1176 (2011).
- M. Connolly et al., arXiv:1201.4029v1 (2012).
- G. M. Rutter et al., Nature Phys. 7, 649 (2011).
- 18.17.25 and 18.17.26 from the Digital Library of Mathematical Functions at http://dlmf.nist.gov/.
- A. Matulis and F. M. Peeters, Phys. Rev. B 75, 125429 (2007).
- M. Tahir and K. Sabeeh, arXiv:0804.4087.
- Eric Akkermans and Gilles Montambaux, Mesoscopic Physics of Electrons and Photons (Cambridge University Press, New York, 2007).
- A. Ferreira, J. Viana-Gomes, Yu. V. Bludov, Vitor M. Pereira, N. M. R. Peres, and A. H. Castro Neto, Phys. Rev. B 84, 235410 (2011).
- Iris Crassee, Julien Levallois, Andrew L. Walter, Markus Ostler, Aaron Bostwick, Eli Rotenberg, Thomas Seyller, Dirk van der Marel, and Alexey B. Kuzmenko, Nature Physics 7, 48 (2011).
- J. Jung, F. Zhang, Z. Qiao and A. MacDonald, Phys. Rev. B 84, 075418 (2011).
- Z. Qiao, J. Jung, Q. Niu and A. MacDonald, Nano. Lett. 11 (8), (2011).
- D. S. L. Abergel, H. Min, E. H. Hwang, and S. Das Sarma, Phys. Rev. B 84,195423 (2011).
- M. Zarenia, J. M. Pereira Jr., G. A. Farias, and F. M. Peeters, Phys. Rev, B 84, 125451 (2011).
- E. V. Castro et al., Phys. Rev. Lett. 99, 216802 (2007).
- J. M. Pereira, Jr., F. M. Peeters and P. Vasilopoulos, Phys. Rev. B 76, 115419 (2007).
- M. Nakamura, E. V. Castro, and B. Dora, Phys. Rev. Lett. 103, 266804 (2009).
- V. Mazo, E. Shimshoni, and H. A. Fertig, Phys. Rev. B 84, 045405 (2011).
- M. Nakamura, and L. Hirasawa, and K.-I. Imura, Phys. Rev. B 78, 033403 (2008).
- M. Nakamura, and L. Hirasawa, Journal of Phys: Conf. Series 150, 022064 (2009).
- M. Zarenia, J. M. Pereira, F.M. Peeters, and G.A. Farias, Nano Research Lett. 6, 452 (2011).
- A. Rycerz, J. Tworzydo, and Carlo Beenakker, Nature Phys. 3, 172 (2007).
- D. Gunlycke and C. T. White, Phys. Rev. Lett. 106, 136806 (2011).
- A. R. Wright and T. Hyart, Appl. Phys. Lett. 98, 251902 (2011).