# Graphene: Superlattices, Topological Kinks, Landau Levels and Tunable Magnetotransport

###### Abstract

We review recent work on superlattices in monolayer and bilayer graphene. We highlight the role of the quasiparticle chirality in generating new Dirac fermion modes with tunable anisotropic velocities in one dimensional (1D) superlattices in both monolayer and bilayer graphene. We discuss the structure of the Landau levels and magnetotransport in such superlattices over a wide range of perpendicular (orbital) magnetic fields. In monolayer graphene, we show that an orbital magnetic field can reverse the anisotropy of the transport imposed by the superlattice potential, suggesting possible switching-type device applications. We also consider topological modes localized at a kink in an electric field applied perpendicular to bilayer graphene, and show how interactions convert these modes into a two-band Luttinger liquid with tunable Luttinger parameters. The band structures of electric field superlattices in bilayer graphene (with or without a magnetic field) are shown to arise naturally from a coupled array of such topological modes. We briefly review some bandstructure results for 2D superlattices. We conclude with a discussion of recent tunneling and transport experiments and point out open issues.

Received (Day Month Year)

Revised (Day Month Year)

Keywords: Graphene, Bilayer graphene, Superlattice, Band structure, Transport, Landau level, Quantum Hall effect, Luttinger liquid

## 1 Introduction

Graphene is a two-dimensional carbon crystal that exhibits novel physics and transport properties due to its excitations
resembling chiral relativistic massless Dirac fermions at low energy.^{1, 2, 3}
Its bilayer cousin, Bernal-stacked bilayer graphene (BLG) has also garnered much interest due to the possible novel broken
symmetry states it could potentially exhibit in the presence of interactions that destabilize the quadratic band touching point
present in a minimal tight-binding
model.^{4, 5, 6, 7, 8, 9, 10, 11, 12} Both graphene
and bilayer graphene are also widely regarded as viable materials for developing new types of device applications due to the
chiral nature of their low energy excitations. This is largely due to the possibility of opening tunable band gaps - by
engineering a relative potential difference on each sublattice ^{13, 14, 15}
or by strain engineering ^{16, 17} in monolayer
graphene or by applying an electric field perpendicular to the layers in
bilayer graphene.^{18, 19, 20}

In this review, we focus on new physics that develops in the presence of slow spatial potential variations in both monolayer
and bilayer graphene. While a spatially varying chemical potential could result in p-n junctions in either case,^{21, 22} a perpendicular
electric field that changes direction as a function of position, specific to bilayer graphene, leads to midgap domain wall
states (or ‘kink’ states) that have topological character.^{23, 24, 25} For an isolated “nanowire”,
electron interactions drive this 1D system into a tunable two-band Luttinger liquid.^{26} These tunable nanowires have also been shown
to act as controllable “electronic highways”.^{27} In both monolayer and bilayer graphene, profiles with periodic potential variations form superlattices that can lead to dramatic modifications in the bulk band structure. Such superlattice potentials have been shown to induce
new Dirac modes at zero and finite energy with tunable velocities and transport
properties.^{28, 29, 30, 31, 32, 33, 34, 35, 36, 37, 38, 39, 40, 43}

We review the physics of such superlattices as well as the effect of a perpendicular orbital magnetic field for three
experimentally accessible field regimes: weak, moderate, and strong fields.^{44} A weak orbital magnetic field
essentially acts as a ‘probe’ of the Dirac modes, while a strong magnetic field overwhelms the effect of the superlattice
potential leading to quantum Hall physics indistinguishable from that of pure graphene. At moderate magnetic fields,
however, we find dispersing Landau levels and interesting field tuning of transport properties. In addition, we discuss the
effect of a magnetic field on the so-called topological ‘kink’ states that form at the interface that separates two region with
opposite interlayer bias.^{45, 46, 44, 47} Throughout this article we discuss recent
developments, ongoing experimental efforts, and open issues.

## 2 Superlattices (SLs) in monolayer graphene (MLG)

### 2.1 Bandstructure of 1D superlattices

For pristine MLG, ignoring spin, the low energy Hamiltonian is given by a matrix at each valley, , where pseudospin labels the two trigonal 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, is the momentum measured from . (We set for convenience.) The quasiparticles in the vicinity of each valley then behave as massless linearly dispersing Dirac fermions, with an energy dispersion .

We focus here on the effect of a smooth SL potential, for which the period of the SL is significantly larger than the interatomic distance. This means we can safely ignore the intervalley scattering of electrons which involves large momentum transfer. We therefore use the above low energy Hamiltonian and focus on the electronic properties of SLs near a single valley. To this end, a 1D SL potential can be modelled as , where with being the SL period, and is the identity matrix in the pseudospin (sublattice) space. To gain a qualitative understanding of the nontrivial phenomena arising from such SLs, such as the anisotropic Fermi velocity renormalization or SL induced band gaps, it is not necessary to assume any specific form for .

The problem of finding the energy spectrum of has been extensively studied. It
was noticed by Park et. al.^{31} that the
chiral nature of Dirac fermions in graphene leads to the anisotropic renormalization of Fermi
velocity near Dirac point. Surprisingly, the Fermi velocity is
not renormalized in the SL direction, but is suppressed perpendicular to the modulation direction,
a counterintuitive effect that is deeply rooted in the chiral
nature of the Dirac fermions in MLG.

#### 2.1.1 Weak SL Potential

For a weak SL potential, the energy spectrum near Dirac point and Fermi velocity renormalization can be well understood from perturbation theory. By expanding the Hamiltonian in the chiral basis, , where and denotes electron and hole states, the kinetic energy part can be brought into diagonal form, while the matrix elements between and , with as the reciprocal lattice vector, is given by , where . Therefore, for states with momenta parallel to the SL direction, , we can show that the full Hamiltonian matrix will consists of two decoupled blocks,

(1) |

where are the electron (hole) energies. Applying the second order perturbation theory, the energy correction at momentum is given by

(2) |

where the sum is carried out in the same block and is understood as the corresponding electron or hole energies. This term is zero because of the linearity of the spectrum, i.e., . This means, along the direction of the superlattice, the Fermi velocity is not renormalized. Therefore, the absence of Fermi velocity renormalization in the SL direction is a consequence of the chiral nature of Dirac electrons and the linearity of the spectrum.

At the MBZ boundary, , the two blocks become exactly identical, which means the energy will be doubly degenerate and energy spectrum is gapless at this point. This result is exact and independent of perturbation theory, and this band touching point will always be present.

Once the momentum is no longer parallel to the reciprocal lattice vector , this nice decoupling will break down and Fermi velocity in the corresponding direction, , will inevitably become renormalized. Within second order perturbation approximation, the renormalization with respect to pristine MLG is^{31}

(3) |

where is the angle between and . Since the right hand side of Eq. (2.1.1) is always negative, except for , the Fermi velocity is decreased from pristine MLG value. In contrast, for artificial electrons with linear dispersion but no chirality in a weak 1D SL, the second order perturbation result for the Fermi velocity is given by

(4) |

which is isotropically decreased and is independent of the direction of . Therefore, the anisotropic renormalization of the Fermi velocity near Dirac cone is truly a signature of chiral low energy excitations in MLG. The left Figure 1 shows the dispersion for a weak 1D SL potential, and corroborates the above results.

^{38}. Copyright (2010) American Physical Society].

For a square barrier SL (Fig.2), the energy spectrum can be found exactly.^{38, 36} By making use of Bloch theorem and matching boundary condition, it can be shown that the energy spectrum can be obtained from the following transcendental equation,

(5) |

Here, we have used the following notation:

(6) |

For a symmetric SL, , or equivalently , Eq. (2.1.1) reduces to

(7) |

where and . For this symmetric case, the energy spectrum is particle-hole symmetric, which means Eq. (7) is invariant under the transformation .

To obtain the behavior near point, we can expand Eq. (7) in small and . The result is

(8) |

From this, we can see that the low energy spectrum is indeed described by an anisotropic Dirac cone, with and

(9) |

#### 2.1.2 Strong SL potential

Since the Fermi velocity perpendicular to the SL direction can be significantly renormalized and even brought to zero for a broad region in momentum space, the energy spectrum becomes dispersionless in this direction and the electrons can be collimated in the SL direction.^{32} Moreover, extra Dirac points can be generated in the energy spectrum for an even stronger 1D SL, which are shown in the right panel of Fig. 1.^{38, 36, 33, 30}

To determine the condition for the emergence of extra Dirac points and also their locations, we consider a symmetric SL and assume and in Eq. (7). Then Eq. (7) reduces to

(10) |

which can be solved by either or . The former condition gives , which is just the original Dirac point. The latter condition leads to with being a nonzero integer. Then, the position of new Dirac points are subsequently found at

(11) |

For asymmetric SL, the energy spectrum is no longer particle hole symmetric and extra Dirac points will appear with nonzero energies.

At the induced Dirac points, the Fermi velocities behave differently from the original Dirac point.^{38} To see this, we can expand Eq. (7) in up to second order and obtain

(12) |

with . Then, at the th extra Dirac point, the Fermi velocities along and directions are given by

(13) |

In contrast, at the original Dirac point, and .

### 2.2 Zero field transport

By assuming a constant relaxation time at the Fermi energy , the dc conductivity for MLG SL can be calculated by^{38}

(14) |

where is the average velocity in the -th direction for -th energy band, is the Fermi-Dirac distribution with .

The results for the various conductivities as a function of Fermi energy are shown in Fig. 3. Fig. 3(a) and 3(b) show and , respectively, for an SL with and . Here, red and dashed blue curve correspond to symmetric () and asymmetric () SL, respectively. From Fig. 3(a), notice that is oscillating when the Fermi energy is below the barrier but increases on the average almost linearly when the Fermi energy is above the barrier. On the other hand, always increases on the average with the Fermi energy. Both and show oscillating behavior and are symmetric (asymmetric) for symmetric (asymmetric) SLs. For , however, there is a dip at the crossing energies of those mini bands. From Eq. (2.1.1), it can be shown that crossing energies occur at at .

Fig. 3(c) and 3(d) show and respectively for a symmetric SL with different SL potential strength. Notice that, at low energies, is smaller than its value in the absence of an SL. Upon increasing SL potential strength, also increases due to the appearance of extra Dirac points. Also, we can notice that is always smaller than at low energies, since near the Dirac point.

Now let us consider conductivities for symmetric SLs at zero temperature and charge neutrality (i.e., and ) with only one Dirac point in the spectrum. As the SL potential is not very strong, the low energy Hamiltonian can be written as , with anisotropic Fermi velocities. It can be shown that the conductivity along and perpendicular to the SL direction is given by

(15) |

where is the universal conductivity of an isotropic Dirac cone. The value of depends on the order of different limits taken, such as vanishing temperature and Fermi energy. However, the form of the result for an anisotropic Dirac cone does not depends how is calculated.

When there are extra Dirac points in the spectrum, by assuming their independence and using Eq. (15), the conductivities now are

(16) |

where counts pairs of extra Dirac points. From Eq. (16) we can see, every time a new pair of Dirac points is generated in the spectrum with and an integer, the conductivity parallel to the SL will shows a dip, while conductivity in the perpendicular direction will diverge.

These results can be confirmed by numerical calculation, using Landauer-Büttiker formalism.^{40} Fig. 4 shows the corresponding numerical results. In the left, the conductivity parallel to the SL direction (top panel) agrees with Eq. (15) when is smaller than the critical value where extra Dirac points are generated. The corresponding Fano factor is 1/3 (bottom panel), which agrees with pseuodiffusive character of transport. Once exceeds the critical value, new Dirac points will emerge and provide new transmission channels in the SL direction. We can see that the conductivity qualitatively agrees with Eq. (16). Since Eq. (16) is based on the assumption that all the Dirac points are independent and have linear dispersion, which is valid in a very small energy region, it is not surprising to see the numerical results depends on both SL potential strength and SL period. Also, the Fano factor is larger than 1/3, indicating that the transport is no longer pseudodiffusive.

^{40}. Copyright (2011) American Physical Society.]

In contrast, conductivity perpendicular to the SL direction is shown in the right figure of Fig. 4. Again, when SL potential strength is smaller than the critical value, the conductivity is well described by the simple picture of Eq. (15). When approaches the critical value, the conductivity shows a peak. For even stronger SL potential, the numerical result agrees with Eq. (16) quite well, which may suggest that the approximations adopted for Eq. (16) are appropriate in the perpendicular direction. Also, the Fano factor is 1/3 for almost all SL potential strengths, except for those critical values where new Dirac points are generated.

### 2.3 Landau levels

In a uniform perpendicular magnetic field, the eigenenergy for pristine MLG in the absence of SL is , where , with . For (i.e., at valley ), the eigenfunctions are given by

(17) |

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

(18) |

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

(19) |

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 at low energy. Recent work has shown these results to be generally
consistent with solving the Harper equation using the full tight-binding model.^{41}
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,^{48} and allows us to numerically access a wide range of magnetic fields. ^{42}
In the weak field regime, our results are consistent with Ref. ^{48}.)
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 (e.g., step-like) SL potentials by including multiple Fourier components. We can then expand the Hamiltonian in the above basis, retaining up to 3000 Landau levels, and diagonalize it 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. 5 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.

#### 2.3.1 Weak field regime

When the magnetic field is weak, (top panels in Fig. 5), 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.

#### 2.3.2 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.5, 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 to lose its simple cosine form.

#### 2.3.3 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

(20) |

which 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. ^{49} 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. ^{49}

### 2.4 Magnetotransport

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,

(21) |

Here, we have set as the Landau level broadening, and are the -th eigenvalue and the corresponding eigenstate of the system which can be expanded in the basis of , where is the LL wavefunctions for pristine graphene. is the velocity operator in -direction and , where is the Pauli matrix. Note that is always true for any state.

Results for dc diagonal conductivities as function of chemical potential are shown in Fig. 6. This can be done by setting the frequency to zero in Eq. (21), 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 consequence of the Fermi velocity renormalization in the absence of magnetic field (see Fig. 6 (a) and (c)). Since and because of the flat band structure, the major contribution to the diagonal conductivities comes from off-diagonal matrix elements, with . Numerically, we have observed that matrix elements of is always larger than those of , which gives rise to the anisotropy in the weak field. In intermediate magnetic field, conductivities still show anisotropy, but with significantly larger than (see Fig. 6 (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. Therefore, the intra-LL contribution to is the smallest at the band edge. 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. ^{50}). 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. 7. For weak magnetic field (Fig. 7 (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 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. 7 (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 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. 8 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,^{50}
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.^{50, 51}
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.

### 2.5 Bandstructure of 2D superlattices

In 2D SLs, the Fermi velocity near the Dirac point is anisotropically renormalized along every direction. Due to the chiral nature of low energy excitation, there are still energy band crossing at the MBZ boundary.

Fig. 9 shows a 2D rectangular SL with muffin tin type SL potential with period and in and directions, and the corresponding energy spectrum. In contrast to 1D SL where Fermi velocity parallel to the SL direction is not affected, Fermi velocity in a rectangular SL is renormalized in every direction. This can be clearly demonstrated by second order perturbation, assuming a weak SL potential strength,^{31}

(22) |

where is the 2D SL potential strength in a circular region of diameter , is the reciprocal lattice vector with and integer, and is the Bessel function. Since can be along any direction, compared to 1D SL where is always along the SL direction, we can see that the Fermi velocity is renormalized in every direction.

The energy spectrum of a 2D rectangular SL also has band crossing points in the middle of an MBZ boundary edge, similar to 1D SL. In addition to these crossing points, at the four corners of the MBZ, energy gap also closes. When similar calculation is carried out for an artificial non-chiral electron, these band crossing points disappear, which truly suggests that they are the consequence of the chirality of low energy excitations in MLG.

Even though band crossing points appear at the MBZ boundary in both 1D and 2D rectangular SLs, the density of states does not vanish at the crossing energy and the newly generated massless Dirac fermions are obscured by other states. However, for triangular SLs, there exists an energy window where the only available states come from the newly generated massless Dirac fermions.^{33} Fig. 10 shows a triangular SL with muffin-tin type SL potential and its corresponding energy spectrum. Again, the Fermi velocity is anisotropically renormalized in every direction. The gap between the first and second conduction bands vanishes in the middle of the MBZ boundary edges, and the density of states also vanishes linearly here. This result will have significant impact on the experiments explained later.^{34, 37, 43, 52}

^{33}. Copyright (2008) American Physical Society.]

When graphene is expitaxially grown on a substrate (i.e., SiC, hexagonal boron nitride (hBN), transition metal surfaces, etc.), the lattice mismatch between graphene and the substrate and also their relative orientation can lead to a 2D SL with large period. Therefore, theoretical results can be tested on these structures. However, previous results are based on an effective Hamiltonian approach which assumes that external potential does not break sublattice symmetry. When such a symmetry breaking effect is taken into account, most of the earlier results will be modified. For example, a gap should open up at Dirac point and minigap should appear at the MBZ where bands are backfolded. Pletikosić et. al.^{34} have observed a minigap in graphene expitaxially grown on Ir(111) surface, which is due to Moiré patterned periodic potential. They could not determine whether the Dirac point is gapped because graphene on Ir(111) is slightly -doped. On the other hand, Rusponi et. al.^{37} showed that, in the presence of sublattice symmetry breaking SL potential, the Dirac point remains intact and, remarkably, the Fermi velocities are anisotropically renormalized and the energy spectrum becomes trigonally warped. This is consistent with the theory of Ortix et. al.,^{52} where it was demonstrated, incommensurate Moiré patterned superstructure preserves the Dirac cone in a renormalized form, with threefold global symmetry due to a substrate-induced trigonal warping. Moreover, additional finite energy Dirac points are also generated, but at different positions of MBZ in contrast to Park et. al.^{33} Since the SL potential also breaks the particle-hole symmetry, the energy spectrum no longer possesses this symmetry, and only in an energy window below the original Dirac point, the newly generated massless fermions are truly Dirac fermions and the density of states can become zero, while for those above, massless fermion states are obscured by the presence of other states. Recently, a scanning tunnelling microscope measurement of graphene on hBN has observed dips in the differential conductance and thus confirmed the existence of finite energy Dirac points.^{43}

## 3 Superlattices in bilayer graphene

We now turn our attention towards 1D electrostatic potential modulations in BLG. In general, the features of the band structure will depend on the details of the superlattice potential^{53, 54}. For the bilayer system, a general modulation can be decomposed into two basic types: i) a chemical potential modulation where both layers sit at same potential and, ii) an electric field modulation where there is a local interlayer bias. If the SL potential is purely of one type, there is a dramatic restructuring of the band structure, particularly at low energy. Notably, the low energy quasiparticles transform from being massive chiral fermions in intrinsic BLG to massless chiral Dirac fermions for certain SL parameters. In both cases, much of the band structure can be understood by appealing to the inherent symmetries and/or to an intuitive effective low energy model. The generation of new zero energy modes has similarly been shown to arise in a periodic array of -function potentials^{35}, twisted BLG ^{55}, and along domain walls in monolayer graphene with broken sublattice symmetry ^{56, 57}.

In this section, we focus on reviewing the band structure of the two rudimentary types of SL in BLG, a chemical potential and electric field superlattice. Both types of SL are of particular interest because each can support the formation of new Dirac points for arbitrarily weak SL strengths, in contrast to SL in the monolayer. In fact, the Dirac points for the electric field SL survive even for strong modulations. A thorough understanding of these two basic SL potentials also provides a firm foundation to understand more generic SL profiles and the formalism reviewed here can readably be applied to more general SLs.

We start here by introducing the low energy Hamiltonian that can be used to study the properties of generic SL of moderate strength. It should be noted that for larger SL potentials, the full tight-binding model is required to correctly describe new features in the band structure. Instances where the full Hamiltonian gives quantitative differences in the band structure will be duly noted. After establishing the formalism, we discuss the band structure generated by a chemical potential and electric field SL in Sec. 3.1 and Sec. 3.2, respectively.

The low energy Hamiltonian for Bernal-stacked BLG can be obtained by expanding its minimal tight binding spectrum near one of the Brillouin zone corners ( points).^{58} When the layer potential (i.e., interlayer potential difference) is not too large, , we find ,^{58} where

(23) |

and , with () being the electron operator on the top (bottom) layer. Here, is the momentum operator, for the Hamiltonian at the valley, m/s is the Fermi velocity, eV is the nearest neighbor hopping integral, Å is the distance between neighboring atoms on the same sublattice (note: where is the nearest neighbor
Carbon-Carbon distance), are the potentials on each layer, and is the interlayer coupling. Unless stated, we set . We will ignore inter-valley scattering assuming the potentials are varying slowly on the scale of , so we only consider the valley (at ). Such an approach has been successfully used to study SLs in monolayer graphene ^{31, 32}.

To diagonalize , we Fourier transform and then make a unitary transformation , , where and . This leads to . Here are energies of electron (hole) states, with an effective mass . This minimal model supports quadratic band touching points at .

When are periodic, we can also Fourier transform the SL potential to obtain , where

(24) |

, and is the angle between momenta and . Our aim is to understand the band structures of SLs described by . We will study 1D SLs with period along , so that the reciprocal lattice vectors, , are integer multiples of , and the mini Brillouin zone (MBZ) boundaries are at .

### 3.1 Band structure of 1D chemical potential superlattices

A chemical potential SL corresponds to the case where . For simplicity, we first consider a step-like potential with with (i) for and (ii) for and use the effective two-band Hamiltonian introduced above. Starting from and increasing the SL strength to moderate values, we observe the following restructuring of the band dispersion (see Fig. 11): i) the zero energy quadratic band touching point splits into two anisotropic Dirac cones located at , ii) further increasing causes the Dirac points to push out towards the MBZ and, iii) upon reach the boundary at , a band gap opens at a critical .^{54, 28, 39} Before considering the band structure for beyond , let us first discuss the formation of the Dirac cones in more detail.

#### 3.1.1 Dirac Cones: Formation

The sequence of semimetal to band insulator with increasing SL strength was shown to not be dependent on any symmetry in SL profile and to be robust even against weak perturbations that vary slowly perpendicular to the principle SL direction.^{28} Given the persistence of the Dirac point, it important to understand why it forms and how it is protected.

The formation of linear band crossing points has been argued to be deeply rooted in the chiral nature of the low energy BLG quasiparticles.^{28, 39} This can be seen from the scattering angle dependence of the matrix elements in Eqn. 24 that arises from the pseudospin structure of the eigenstates. For states with momenta parallel to the modulation direction, or , the off-diagonal matrix elements vanish; the electron and hole states decouple, so that a particle in an electron (hole) state can only forward/back-scattering of the SL potential to another electron (hole) state. Since all such electron (hole) states within the first MBZ in an extended zone scheme only mix with electron (hole) states of higher (lower) energy, the energy of the conduction (valence) band will be globally shifted down (up). This results in two level crossings along the modulation direction, which are protected by the chirality of the low energy BLG quasiparticles. If this electron-hole decoupling was true for all momenta, we would see the two parabolic bands crossing on a full circle in the MBZ, but going to momenta leads to electron-hole mixing that is linear in ; this results in an avoided level crossing and the robust emergence of two Dirac cones in the MBZ.

#### 3.1.2 Dirac Cones: Properties

The location and velocity anisotropy of Dirac cones, as well as the critical modulation amplitude to gap them out, can be estimated using perturbation theory in . The second order energy correction of states with is Since while in the MBZ, this correction is always negative (positive) for electron (hole) states, as expected.

Thus, the two bands will intersect and cross linearly at momenta , where For weak modulations, , and keeping only , we estimate . For a step profile, , and contributions are small.

For small away from the level crossing point, we can estimate the electron-hole mixing term using perturbation theory and we find that the resulting eigenstates have energies . The crossing points at are thus really massless Dirac points in the full MBZ. We find velocities , and for the anisotropic linear dispersion.

Once these Dirac nodes reach the MBZ boundary, Bragg scattering between them opens up a full gap. The critical potential strength, for this is roughly estimated by setting , which yields . For a step profile, with , we find which is close to the numerical result .

#### 3.1.3 Strong Potentials and Electron Screening

For larger SL potentials beyond , it becomes necessary to employ the four band model to capture a the effects of the high energy bands. Increasing above , the band gap continues to grow until a maximum value is reached and then begins to decrease before finally closing by forming two new pairs of Dirac cones. For even greater , two of the four Dirac cones merge at and become gapped, but remaining two Dirac cones retain the semi-metallicity of the system.^{28, 39}

Tan et al.^{28} also performed a self-consistent tight-binding calculation to describe the higher energy effects of the entire band structure and to determine the effects of interactions on the band dispersion. Besides confirming that the simple single particle low energy model correctly describes the qualitative features of the band structure, it showed that it is possible to account for screening at the Hartree level by a dielectric constant . Hence, the main effect of electron interactions is to screen the external SL potential, therefore increasing the critical SL potential required to open a band gap.

### 3.2 Band structure of 1D electric field superlattices

When BLG is subjected to an electric field SL the potentials on the two layers are such that . In contrast to the chemical potential SL discussed above, the band structure is sensitive to the form of the SL profile^{53}. To illustrate this, Killi et al.^{39} considered a more general periodic potential, with for , and for , where we have kept the average potential to zero.

If the parameter the SL potential is symmetric. A numerical calculation of the band structure show a pair of anisotropic massless Dirac cones forming at zero energy at , as seen in Fig. 12 (left panel)^{38, 39}. Here, the zero energy Dirac cones lie along the direction perpendicular to the modulation as opposed to along it for chemical potential SL. Two additional anisotropic Dirac cones are also present at high energy, one in the valence band the other in the conduction band. Irrespective of the strength of the SL potential, these Dirac points are pinned MBZ boundary at (or equivalently ).

When a band gap opens at all of the Dirac points. This suggests that the protection of the Dirac points is governed by a symmetry of the SL profile that is broken when . It is found that the relevant symmetry corresponds to a generalized parity operator that transforms followed by exchanging the two layers of BLG, and the Dirac points persist as long .

In the letter by the present authors,^{39} a simple intuitive picture for understanding all of the features of the Dirac cones was proposed. The idea is to view the SL potential as establishing a periodic array of ‘kink’ and ‘anti-kink’ steps in the potential profile where the parity of interlayer bias reverses. Since it was shown by Martin et al.^{23} that topological zero energy modes are confined along an isolated kink (or anti-kink), it is possible to construct an effective low energy theory evolving these modes. Moreover, the band structure of the SL should be entirely dictated by these modes fore energies below where the bulk states are gapped. As a prerequisite to presenting and analyzing the low energy effective model, it is necessary to understand some of the basic properties of the 1D kink modes. The next we provide a brief overview of these states before returning the construction of the low energy model of the electric field SL.

#### 3.2.1 Kink in the electric field: Soliton modes

For a general potential profile with and , the bulk region far from has a gap , while along this interface, localized ‘topological’ edge modes emerge that are analogous to those in quantum hall systems.^{23} These modes can be thought of as forming chiral 1D quantum wires, since states from opposite valleys are counterpropagating (this follows from the Berry curvature about each valley having opposite sign). With respect to the two band Hamiltonian, the kink interface marks a region where the mass of the quasiparticles, , changes sign.

Solving the full tight binding model for a single kink yields the dispersion depicted in Fig. 13. A single kink interface generates right moving subgap modes in one valley and two left moving modes in the opposite valley (labeled in red). For an anti-kink profile, the dispersion is identical except the modes velocities are reversed in each valley.

In terms of the low energy, the eigenfunctions are then of the form

(25) |

with corresponding eigenvalues of and of the operator , respectively. Solutions with eigenvalues with even symmetry belong to the lower -band while solution with odd symmetry belong to the upper -band.

Increasing the interlayer bias strength results in two important effects that can be seen qualitatively in Fig. 13: i) the Fermi velocity of the two bands is enhanced, and ii) the wavefunctions become more confined to the interface. With respect to the overall width of the wavefunction, a simple scaling analysis suggests that the wavefunction width should go as .

The large width of the wavefunction transverse to the wire direction strongly suppresses the bare backscattering terms due to electron-electron interactions.
A similar effect also seen in wide carbon nanotubes,^{59} and it leads to a dominance
of forward scattering processes, where a simple bosonization analysis predicts a spin-charge
separated gapless Tomonaga-Luttinger liquid.
We thus expect a relatively large energy window where interactions drive the 1D
modes even in these ‘kink’ modes in bilayer graphene into such a Luttinger liquid.
Remarkably, such a bosonization analysis arrives at a novel two-band Luttinger liquid with
tunable mode velocities and tunable Luttinger parameters.^{26}

Naively, it would appear that the localized kink states may not be robust because they are not topologically protected. However, Qiao et al,^{27} performed an extensive study into various potentially detrimental mechanisms. Through numerical conductance calculations and examining the LDOS, they showed the low energy states are remarkably robust to both short and long range disorder, and even to abrupt changes in the interface direction. Although quantized conductance is not unlikely, the mean free path could be as large as a hundred microns for relatively clean samples. The mechanism which quantitatively leads to a strong suppression of
backscattering is again the large wavefunction spread.

It has also been observed that the kink-modes are also quite robust when subjected to a magnetic field.^{46, 45, 44} This is due to the strong magnetic field induced confinement of the wavefunctions.^{46} Interestingly, by coupling a pair of coupled
kink and antikink modes and applying a magnetic field, the current in the kink flows in one direction, opposite to the direction of flow in the antikink. Moreover, all the dispersing modes have the same valley index, making this a potential valley filter.^{44, 60} Recently, it has been proposed that at electron interactions form a charge density pattern in the vicinity of a kink state, which provides a key signature of quantum Hall ferromagnetism.^{47}

Before turning our discussion back over to superlattices, we close this section by emphasizing that similar localized kink states are also expected to form naturally in the presence of charge impurities and in strongly correlated phases where the layer symmetry is spontaneously broken. In the case of the former, charge impurities close to the surface of the sample can generate a local electric field strong enough to induce an interlayer bias. In uniformly biased BLG or where there are multiple charge impurities in close vicinity (on opposite layers or with opposite charge), this can cause the interlayer bias to reverse, generating kink states – a point to be elaborated on in the conclusion. In the latter case, any state with spontaneously broken layer symmetry will naturally form domain walls separating regions with opposite interlayer bias.^{9, 10} Again, kinks states are expected to from percolation networks that permeate throughout the bulk.

#### 3.2.2 Electric Field Superlattice: Effective Model

Equipped with an understanding of the properties of the soliton kink/anti-kink modes, it is now possible to describe how to construct a low energy effective model. ^{39}
To begin, first consider the dilute limit where the period length is much longer than the characteristic spread of the soliton modes, (i.e. ). Each kink supports two (ignoring spin) unidirectional dispersing soliton modes while each anti-kink supports two oppositely moving modes per valley, as shown in Fig. 14 The counterpropagation of the kink and anti-kink modes results in four band crossing points about each K-point, two at zero energy between bands with the same symmetry (- and -) and two at finite energy between modes with opposite symmetry (- and -). As described below, when the wavefunctions of neighbouring soliton modes couple (i.e. ), Dirac cones precipitate precisely at the band crossing points.

At energies and momenta in the the vicinity of any one of the band crossing points the system looks as if it were an array of 1D chiral ‘wires’ lying along the kinks and anti-kinks of the SL. Each wire supports modes that flow in opposite direction to its two neighbors. Now, as the wavefunctions of these modes begin to overlap the electrons can hop between neighbouring wires.

With this in mind, let us consider the - modes at zero energy and at a momentum (away from ). The hopping between neighboring wires along is then between states which have opposite velocities (since it is between a kink and an anti-kink edge state) and it is between a p-wave like state (-odd) and an s-wave like state (-even) (see Sec. 3.2.1).

Careful attention must be made to get the correct form of transfer integral that describes the hopping between the wires. The sign of the hopping can be deduced by taking the wave functions of the -band and -band as having s-wave orbital and p-wave orbital character, respectively, and noting the sign of the overlap between the wires. An illustrative picture for two of the band crossing points is shown on the right in Fig. 13. Hence, the hopping between wires at a zero energy band crossing point is staggered and uniform for the finite energy band crossing points. Further details are provided in our previous work. ^{39}

Using the index to label the wires, the interchain hopping parameter will then alternate as for equally spaced wires and as (with ) if pairs of wires are closer to each other. Linearizing the dispersion at the crossing point, and letting denote the velocity of the linearized modes,

(26) | |||||

where is the location of the crossing point in the single kink or antikink problem, and annihilates an electron on wire with momentum . Let