Moire superlattice effects in graphene/boron-nitride van der Waals heterostructures

Van der Waals heterostructures of graphene and hexagonal boron nitride feature a moir\'e superlattice for graphene's Dirac electrons. Here, we review the effects generated by this superlattice, including a specific miniband structure featuring gaps and secondary Dirac points, and a fractal spectrum of magnetic minibands known as Hofstadter's butterfly.


Introduction
Van der Waals heterostructures formed from stacks of two-dimensional crystals have recently attracted a great deal of interest, as the available materials cover a broad spectrum of physical properties [1][2][3][4]. The first and, by now, the most developed of these heterostructures is that formed by graphene placed on hexagonal boron nitride (hBN). Both materials are made from a single layer of atoms arranged in a honeycomb structure. However, whereas graphene [5,6] is a gapless semiconductor, hBN is an insulator with a band gap of approximately 6 eV [7]. For this reason, hBN was first used as a substrate to preserve graphene's electronic properties [8]. Nevertheless, the small difference between the lattice constants of the two crystals, and any misalignment between their respective crystallographic axes, generates a large-scale quasi-periodic hexagonal pattern, known as the moiré superlattice (mSL) [9][10][11]. In this review, we focus on the qualitatively new phenomena generated by the presence of the mSL, such as the formation of electronic and magnetic miniband structure. Since the perturbing effect of the mSL on graphene's Dirac electrons decreases with increasing misalignment between the hBN and graphene lattices, we concentrate on the well-aligned heterostructures (misalignment angle θ 2 • ).
The mSL provides a long-wavelength periodic potential, which leads to the Bragg scattering of electrons in graphene. This generates a miniband structure with the most pronounced effects found at the edges of the first and second miniband, where the mSL perturbation leads to the generation of secondary Dirac Points (sDPs). The signature of these sDPs has been observed, both in the density of states [11,12] and in transport measurements [13][14][15], as a repetition of the behaviour of the primary Dirac point when the chemical potential is tuned to the energy of the sDP. Notably, the signatures of the sDPs are much more pronounced in the valence band than the conduction band, which is a manifestation of the electron-hole symmetry breaking produced by the mSL perturbation.
Further spectral reconstruction occurs at the neutrality point of graphene. While the first transport measurement [13] reported no detectable band gap, another experiment [14] reported a band gap of approximately 30 meV. The exact details of the spectrum at the neutrality point appear to depend on the architecture of the device, namely whether the graphene is encapsulated from the top with an additional, intentionally misaligned, hBN layer, as was the case in the former experiment, or left unencapsulated, as it was in the latter. A further study [16] showed that the graphene and hBN layers in almost perfectly aligned, yet unencapsulated devices, undergo local displacements subject to their mutual adhesion, so that they assume their preferred stacking order for a large portion of the moiré unit cell. In these devices transport measurements revealed significant band gaps, consistent with the expectation [17][18][19][20] that graphene's sublattice symmetry is strongly broken in this phase. Conversely, the lack of detectable band gap, and absence of strain signatures in the Raman spectroscopy data, was interpreted as evidence that the additional randomised adhesion forces due to encapsulation with a misaligned hBN layer suppress the local displacements of the crystals, and therefore (on average) restore the sublattice symmetry. Also, recent magneto-optical spectroscopy data [21] for unencapsulated epitaxially grown heterostructures has been fitted to a gapped Dirac spectrum with a gap of 38 meV and a transition dependent Fermi velocity, the latter interpreted as a signature of the electron-electron interaction in the graphene/hBN heterostructure [19,20,22].
When the graphene/hBN heterostructure is placed in a magnetic field that is strong enough that the magnetic length is comparable to the periodicity of the mSL, an intricate fractal spectrum of magnetic bands [12][13][14]23] (sometimes referred to as "Hofstadter's butterfly" [24]) is formed. The physics of this phenomena was explored theoretically many years ago [25,26]. However, accessing the fractal electronic spectrum experimentally is more challenging because, for a typical crystalline periodicity, the magnetic fields required to access this regime are unobtainable. Low-field oscillations and internal structure within Landau levels [27][28][29][30][31][32][33][34] have been observed using semiconductor heterostructures, but the technological challenges involved in producing the structures and accessing the two-dimensional electron gas limited further studies. In contrast, graphene/hBN heterostructures are uniquely suited to investigate the "Hofstadter butterfly" spectrum because of their high electronic quality, ease of doping, and convenient range of the moiré wavelength (A 14nm). One striking feature in these heterostructures is the systematic reappearance of tertiary Dirac points at the edges of the magnetic minibands formed at rational values of the magnetic flux [13,35,36].
Graphene/hBN heterostructures have also been studied using other techniques, such as optical absorption [37] and Raman [38] spectroscopy. For the former, the reconstruction of the electronic spectrum at the edge of the first miniband and the modification of optical selection rules, result in a modulation [37,39] of the otherwise flat absorption spectrum of graphene [40]. In the case of the latter, the width of the Raman 2D peak (for unencapsulated samples) was found to increase linearly with moiré wavelength, which was interpreted in terms of increasing strain in the graphene layer [38]. It has also been calculated that the formation of minibands in the graphene/hBN heterostructures should impact graphene's plasmon dispersion [41]. Finally, the moiré minibands in graphene/hBN heterostructures may have some interesting topological properties including nonzero Berry curvature [18,42]. In an applied electric field the Berry curvature generates an 'anomalous' contribution to the electron's velocity [43], which has been observed in recent transport measurements [44].
Geometry of the moiré superlattice. When two isostructural crystals with a small lattice mismatch, δ, are placed on top of each other with a small misalignment angle, θ , between their respective crystal directions, the beating of the two crystalline periodicities produces the large-scale modulation of their local stacking known as the mSL (see Fig. 1 upper panels). For graphene, this geometry-based effect occurs on a variety of different hexagonal substrates/surfaces [71][72][73][74][75][76][77][78][79], and can be observed, for example, using scanning tunnelling microscopy or conducive atomic force microscopy [4, 9-11, 13-16, 23, 69, 70]. In Sec. 5, we discus the Figure 1 Top: examples of the moiré pattern arising from graphene (blue) placed on hBN (grey), for two different misalignment angles, θ = 0 and θ = 1 • . Bottom: the Brillouin zones of graphene (left) and the mSL (right), with various high symmetry points and reciprocal lattice vectors labelled.

Review Article
Ann. Phys. (Berlin) 527, No. 5-6 (2015) magnifying effect in which defects in either of the two crystal layers (such as dislocations or wrinkles) generate a large partner defect in the mSL [80][81][82]. Here, we assume that graphene and its substrate are perfect crystals that are placed rigidly on top of each other. Then, the six (m = 0, . . . 5) shortest reciprocal vectors of the mSL, are obtained from the graphene reciprocal lattice vectors g m =R 2π m/6 (0, 4π √ 3a ) by the subtraction of the substrate reciprocal lattice vector g m = (1 + δ) −1R θ g m (whereR θ is the matrix for anticlockwise rotation, a is the graphene lattice constant,ẑ is the unit vector in the z-direction, and θ, δ 1). The corresponding real space lattice vectors of the mSL are given by, where a m =R 2π m/6 (a, 0) are the graphene lattice vectors. It should be noticed that the mSL described by vectors A m is distinct from the crystallographic lattice of the two layer-system, the latter requiring a perfect translational symmetry that is only present when the two layers are commensurate. As shown in the top row of Fig. 1, as θ increases, the primitive lattice vectors of the mSL are rotated by ≈ − tan −1 (θ/δ) and their length decreases as A ≡ |A m | ≈ (δ 2 + θ 2 ) −1/2 a. For the graphene/hBN heterostructure, δ ≈ 1.8% [83] so that mSLs with lattice constants as large as A ≈ 14 nm are possible for θ = 0. The energy of the first miniband edge is approximately vb/2 ≈ v(δ 2 + θ 2 ) 1/2 |g m |/2, which is as low as 170 meV for θ = 0 (where b ≡ |b m |, = 1, and v is the Dirac velocity). Importantly, this energy corresponds to a carrier density of n ∼ 2 × 10 12 cm −2 in the graphene flake, so that electrostatic doping can be used to shift the Fermi level to the miniband edge, and the strong spectral reconstructions found at the edge of the superlattice Brillouin zone (sBZ) can be probed in transport measurements. To compare, the characteristic energy vb/2 increases rapidly with the misalignment angle, so that for θ = 2 • almost five times greater doping is required to fill the first miniband.
Superlattice perturbation for Dirac electrons. To model the low-energy band structure (| | vb) of the graphene/hBN heterostructure, we use the fact that hBN has a large band gap ( hBN ∼ 6 eV [7]), to project the full two-layer Hamiltonian onto the Hilbert space of the graphene π orbitals only. Then, the dominant effect of the substrate on graphene Dirac electrons can be modelled by scattering processes using the six shortest non-zero reciprocal lattice vectors of the hBN substrate (sometimes referred to as the "first star" of reciprocal lattice vectors). One such process, in which an electron (red point) at graphene's K Brillouin zone corner is scattered by −g 4 , is displayed in the left lower panel of Fig. 1. Hence, upon addition of a graphene reciprocal vector, the mSL perturbation provides intravalley scattering by the simplest harmonics of the mSL (b 4 = g 4 − g 4 in this case). Scattering processes involving higher reciprocal lattice vectors of the substrate (or, equivalently, of the mSL) are suppressed, as they necessarily depend on the overlap between higher Fourier components of the graphene 2p z orbital (see section 2.3). However, at interlayer distances, z, comparable to the graphene-hBN interlayer spacing, d, these Fourier components decay rapidly as a function of the in-plane momentum q, whereψ(q, z) is the two-dimensional in-plane Fourier transform of the graphene 2p z orbital, and a 0 is the effective carbon Bohr radius. In addition to the dominance of the simplest mSL harmonics in the perturbation, the hBN substrate can only affect the graphene electrons via three distinct mechanisms [84]; (i) an electrostatic potential, which does not distinguish between the two carbon sublattices, (ii) a sublattice-asymmetric part of the potential, and (iii) spatial modulation of the nearest neighbour carbon-carbon hopping amplitude. Each of these can be thought of as made of two contributions, either symmetric or antisymmetric under the in-plane spatial inversion. It can be argued [84] that in the two limits where either, both boron and nitrogen sublattices perturb the Dirac electrons with almost the same strength, or, the dominant perturbation arises from one sublattice only, the inversion symmetry of the system would only be weakly broken. Then, the mSL potential can be modelled as a combination of a dominant inversionsymmetric part with the addition of a small inversionasymmetric perturbation. Accordingly, the mSL perturbation can be parametrised by three phenomenological parameters, U + 0 , U + 3 and U + 1 , each controlling the strength of the inversion symmetric component of modulation mechanisms (i), (ii) and (iii) described above. For such perturbations, the electronic miniband spectra in either graphene's valence or conduction bands can Panel (d) also shows the fractal spectrum of magnetic minibands, with an enlargement of the indicated region. The full dispersion of the magnetic minibands at φ = φ 0 /2 is also presented with a colour map displaying the Berry curvature of the lower band.
be classified into three groups depending on the mutual arrangement of the first and second miniband: either (a) they do not overlap and are connected by a single isotropic sDP, (b) they do not overlap and are connected by a triplet of anisotropic sDPs, or (c) the minibands overlap. The regions in the parameter space of the three inversion-symmetric parameters, corresponding to the valence band spectrum of a particular type, are shown in different colours (orange, green and white, respectively) in Fig. 2(a). Representative examples of spectra from each region are shown in Fig. 2(d,c,b). Note, that while a purely inversion-symmetric mSL perturbation generates both gapless primary DP and sDPs, the addition of a finite inversion-asymmetric component will open gaps at both. The experimental evidence [11][12][13][14][15] clearly points to the existence of sDPs in the valence band, with recent temperature-dependent measurements attributed to the spectrum featuring only a single isotropic sDP [12].
In a perpendicular magnetic field, a hierarchy of magnetic minibands and gaps is generated for each magnetic field providing a rational fraction of the flux quantum φ = φ 0 p/q (here φ 0 = h/e is the flux quantum) per moiré unit cell of the graphene/hBN heterostructure. For a weak magnetic field (φ 0.2φ 0 ) these can be traced to the weakly broadened Landau levels of the zeromagnetic field bandstructure (see Fig. 2d). In contrast, in stronger magnetic fields (φ 0.2φ 0 ), the edges of many of the consecutive magnetic bands can be described in terms of a weakly gapped Dirac spectrum [13,35]. One example is shown as an inset in Fig. 2(d) with the Berry curvature of the lower magnetic miniband shown as a colour map. The fractal spectrum surrounding each magnetic miniband can be interpreted [85] in terms of the weakly broadened Landau levels of these next generation Dirac electrons (see right panel of inset). This can be seen best in the cases when the gap between consecutive minibands is small, and the largest gaps are around
Finally, we note that superlattice perturbations for Dirac electrons have been discussed in a variety of different contexts, some of which pre-date the realisation of the graphene/hBN heterostructure. Many theoretical works have investigated the influence of one or two-dimensional electrostatic potentials on graphene electrons , with the former situation realisable using patterned gates [109,110]. Magnetic and pseudomagnetic field superlattices (the latter arising from periodically strained graphene) have also been extensively studied [95,[111][112][113][114], with steps towards experimental realisation [115,116]. There has also been significant work on the aligned heterostructures of bilayer graphene with hBN, including the observation of Hofstadter's butterfly in transport measurements [23]. Theoretically, the mSL perturbation of this heterostructure can be modelled in a similar manner to the monolayergraphene/hBN heterostructure, except that the perturbation is felt much more strongly by the graphene layer which is closest to the hBN [117,118]. Because of this, the inversion symmetry of this heterostructure is broken, typically leading to gaps both at zero energy and at the edge of the first miniband [117]. In a magnetic field, the broken inversion symmetry manifests itself as a strongly broken valley symmetry in the magnetic miniband structure [118]. Other systems in which miniband structure is generated by a mSL include twisted bilayer graphene [119][120][121][122][123][124][125][126][127][128], graphene with almost commensurate √ 3 × √ 3 hexagonal crystals [129], and graphene on metal catalysts such as Ir(111) [130].

Phenomenological Hamiltonian
To describe the effect of the moiré perturbation on the graphene electrons, we assume that the local form of the Hamiltonian is determined by the local relative displacement, of atoms in the underlay with respect to the graphene layer. The first two terms in u(r) describe the effect of a finite lattice mismatch, δ, and misalignment angle, θ , between the two lattices. The final two terms account for additional displacements in either the hBN layer (u hBN ) or the graphene layer (u Gr ) which may be caused, e.g., by a defect in one of the crystal layers, or spontaneous lattice distortion produced by their mutual adhesion. We also recognize that the interlayer separation between graphene and hBN may vary locally, by a small amount z, for regions of the mSL with different local values of u(r). Below, this is included as an overall prefactor to the mSL potential h(u(r)) ∼ 1 + O( z). Then, consistent with the dominance of the first star of moiré harmonics in the mSL potential, a general [82] phenomenological Hamiltonian of the graphene/hBN heterostructure reads, Above, the Pauli matrices σ i=1,2,3 and σ = (σ 1 , σ 2 ) act on the electron amplitudes on the graphene A and B sublattices and ζ = +1 (ζ = −1) is used for graphene's K (K ) valley. The first term inĤ is the Dirac Hamiltonian, with momentum p = −i∇ + e A taking into account the magnetic vector potential A related to the external magnetic field using B = ∇ × A. In this term we also include one of the two dominant strain-induced effects on graphene electrons [5,[131][132][133][134] using the pseudo-vector potential A = − √ 3η0 2ea (∂ y u Gr y − ∂ x u Gr x , ∂ x u Gr y + ∂ y u Gr x ) T (with η 0 = ∂lnγ 0 /∂lna used to express the change in the nearest neighbour hoping γ 0 with the change of the lattice constant, and u Gr = (u Gr x , u Gr y )). The second term describes the other strain-induced effect: a scalar potential accounting for modification of the on-site energy 2 p of the graphene 2p z orbitals due to the local changes of bond lengths. The remaining terms inĤ describe the mSL perturbation, where δĤ + (δĤ − ) contains terms symmetric (anti-symmetric) under the in-plane spatial inversion. These two contributions are written using two sets of phenomenological parameters, U + i and U − i , which characterise the various perturbation mechanisms described in the introduction.
When the atomic rearrangements are suppressed, u Gr = u hBN = z = 0, e.g. by encapsulation of the graphene layer, the Hamiltonian is simplified to [84], Indeed, since the higher harmonics of the mSL, induced by either u Gr , u hBN , or z = 0 in Eq (4), only lead to second order corrections to the energy of the lowest energy minibands, they can often be neglected. Then Hamiltonian (4) can be re-parametrised as Hamiltonian (5), albeit with two additional non-spatially dependent terms [18]: (i) a mass term ζ mσ 3 which opens a gap in the primary Dirac point and (ii) a trivial term which causes a shift of the energy scale. For this reason, unless stated otherwise, we limit ourselves to the discussion of Eq. (5). Also, for convenience [84] we usually choose θ = 0. The bandstructure of Hamiltonian (5) calculated in either valley obeys the c 3v symmetry. Moreover, if either the time-reversal (i.e. B = 0) or spatial inversion symmetry (i.e. U − i = 0) is satisfied, the spectrum of Hamiltonian (5) obeys the symmetry K+ p = K − p so that we can limit the discussion of minibands to the K valley. Moreover, using the commutation properties of σ i , the following symmetry relations for the spectrum can be obtained, The first equality above allows one to obtain the phase diagram for the behaviour of the first miniband edge in the conduction band from that of the valence band ( Fig. 2(a)). Finally, a simultaneous transformation of all the perturbation parameters [135], is equivalent to a coordinate shift H (r) → H (r − 4π 3b 2 b 0 ) and therefore leaves the band structure unchanged. Figure 2 shows representative spectra calculated perturbatively [84] for inversion-symmetric mSL perturbations in a basis of plane wave states from graphene's K valley. All of these spectra show a gapless Dirac spectrum persisting near the conduction-valence band edge with an almost unchanged Dirac velocity. In contrast, the inversion-asymmetric terms U − i are able [136] to open a minigap, 0 , at the primary Dirac point [39,117]. To second order in the perturbation parameters this is given by,

Moiré minibands for undeformed heterostructures
More sophisticated calculation shows that this gap may be enhanced by distortions of the graphene and hBN lattices, and also the electron-electron interaction [17][18][19][20]22]. For the point μ = b 0 /2 on the edge of the first sBZ, zone folding brings together two degenerate plane wave states, |μ + q and |μ + b 3 + q . Splitting of these degenerate states by the moiré potential in Eq. (5) can be studied using degenerate perturbation theory. The corresponding 2 × 2 matrix, expanded in small deviation q of the electron momentum from the μ-points reads, For the inversion-symmetric perturbation, the dispersion calculated using Eq. (8) contains an anisotropic sDP [11,84,90,101] with Dirac velocity component ≈ 2U + 0 /b in the direction of the sBZ edge and ≈ v in the perpendicular direction. A similar Hamiltonian to Eq. (8) is obtained for the other two non-equivalent edges of the first sBZ, so that three of these sDPs are visible between the first and second valence miniband in Fig. 2(c). The gap at the sDP, opened by the inversion-asymmetric terms is, For the points κ = (b 4 + b 5 )/3 and −κ, zone folding brings together three degenerate plane wave states, |ξ (κ + q) , |ξ (κ + b 1 + q) , and |ξ (κ + b 2 + q) (where ξ = ± for the two inequivalent sBZ corners), whose splitting is determined bŷ The inversion-symmetric terms inĤ ξ (κ+q) partially lift the ξ κ-point degeneracy into a singlet with energy ( svb √ 3 − 2w ξ ) and a doublet with energies ( svb √ 3 + w ξ ), so that

Review Article
Ann. Phys. (Berlin) 527, No. 5-6 (2015) a distinctive sDP [84,101,137] characterized by Dirac ve- [135] is always present at ±κ somewhere in the spectrum. This single isotropic sDP is visible in Fig. 2(d) at the first valence miniband edge. Once again, the inversion symmetric terms open a gap at the sDP,

Microscopic models for moiré perturbation parameters
Interlayer hopping models. To describe the interlayer coupling in graphene/hBN heterostructure, several studies [17,18,118,136] made use of the similarity between this system and twisted bilayer graphene. In the latter [119][120][121][122], the electronic structure can be described by a bilayer-like Hamiltonian, in which the intralayer blocks are given by the Dirac Hamiltonian and the interlayer blocks describe the modulation of the interlayer coupling as a function of the position within the mSL [121]. Applied to the K valley of graphene placed rigidly on top of hBN, such a Hamiltonian takes the form Above, the Dirac Hamiltonian is used for the graphene layer,Ĥ hBN describes the hBN layer with B and N characterising the on-site energy of the boron and nitrogen sublattices, andT (r) describes the spatially varying interlayer coupling, with γ B and γ N the hopping integrals from graphene to the boron and nitrogen sites respectively. For the energies of interest, | | | N |, | B |, the Hamiltonian can be projected on the Hilbert space of the graphene layer, so that the perturbation is parametrised by Eq. (4) with, where, Setting V + = 0.064vb and V − = 0 results in the spectrum displayed in panel (d) of Fig. 2. Interestingly, even if the inversion asymmetric perturbation is included in this model using V − = 0, the primary Dirac point remains gapless ( 0 = 0 in Eq. (7)), although a gap will be opened in the sDP as expected.
Intralayer potential models. A complementary description of the mSL perturbation for electrons in the graphene/hBN heterostructure concentrates on the significant [138] electrostatic potentials generated by the presence of two distinct atomic species in the underlay. In contrast to graphene, in which every carbon atom (atomic number 6) provides one π-electron, in hBN the boron atoms (atomic number 5) provide no π-electrons whereas the nitrogen atoms (atomic number 7) provide two π-electrons each. Then, the electrostatic potential in the underlay can be modelled as a trigonal lattice of +2|e| point charges mimicking the core charges of the nitrogen atoms, compensated by a homogeneous background charge density mimicking the dispersed cloud of the π-electrons [84]. The matrix elements of the resulting potential, taken between sublattice Bloch states i and j (i, j = A or B), acting on the low energy Dirac spinors of the graphene K valley, are given by the long wavelength components of Here, R N are the positions of the nitrogen sites, L 2 is the total area of the graphene sheet, 3.22Å ≤ d ≤ 3.5Å [139] is the graphene-hBN separation and K ,i (r, z) are the Bloch wavefunctions of graphene 2p z orbitals exactly at the K point. The long wavelength components of Eq. (10) are extracted [84] by re-writing the Bloch wave functions using the in-plane 2D Fourier transform of the graphene 2p z orbitalsψ(q, z), and the dominance of the simplest moiré harmonics in the mSL potential is set by the rapid decay ofψ(q, z) for |q| > |K| at z = d (see Eq. (2)). Interestingly, this model is parametrised by Eq. (9) except with, The relation V − = 0 is prescribed by the inversion symmetry resulting from the assumption that only the nitrogen atoms are responsible for the mSL potential. Nonzero values for V − are obtained when a charge on the boron sites is included, which can be obtained from Eq. (9) using the coordinate shift, Eq. (6). Ref. [137], also models the effect of potentials generated by the hBN layer on the graphene π-electrons, and obtainins a similar mSL perturbation but without the modulation of the nearest neighbour hopping between carbon atoms. Numerical models. Many studies of the graphene/hBN heterostructure have employed numerical models based on ab initio and/or numerical tight-binding approaches. One obstacle to the application of such techniques arises from the fact that no true (perfectly translationally invariant) unit cell exists for this system, due to the incommensurability between the graphene and hBN lattices. Several early ab initio studies avoided this problem by assuming the two crystals to be lattice-matched, for example by contracting the hBN lattice and assuming perfect rotational alignment. The graphene in this lattice-matched system typically has a strongly broken sublattice symmetry, resulting in predictions of large band gaps ∼ 50 meV opened at the primary Dirac point [139][140][141][142][143]. However, already calculations by Sachs and co-authors [143] showed that the van der Waals adhesion of the graphene to the hBN underlay is insufficient to cause the system to become fully lattice-matched.
Other ab initio models of the graphene/hBN heterostructure have used a large commensurate unit cell (e.g. 55 carbon atoms on 56 hBN atoms) to approximate the incommensurate moiré unit cell [142,144]. Also, numerical tight-binding calculations using a large commensurate unit was employed to explain some of the experimental work [14,23], and more recently [118] this approach was found to agree with a continuum model similar to Eq. (5).
A further development is detailed by Jung and coauthors in Ref. [138], with similar approaches employed in Refs. [19,20,143]. Using this methodology, the numerical simulation is performed using unit cells containing only four atoms in total (the two atoms from each of the graphene and hBN unit cells) and a forced commensuration. However, the incommensurability is then included by simulating many such unit cells with a graphene-hBN coordination specific to each local area in the mSL supercell. Then the mSL perturbation can be parametrised using the dominance of the simplest mSL harmonics in the mSL perturbation.
Models incorporating spontaneous lattice deformations. As mentioned in the introduction, a relaxation of the atomic positions of the crystal layers in unencapsulated highly aligned graphene/hBN heterostructure has been observed experimentally [16]. Here, the local displacements of the two layers results from a competition between the graphene-hBN adhesion landscape and the elasticity of the graphene and hBN layers. The former contribution varies on the atomic scale and is found to be minimised for a commensurate AB stacking in which the boron atom is situated directly beneath a carbon site, and the nitrogen directly beneath the hexagon centre [139,141,143,145]. In contrast, the elastic contribution prefers the incommensurate arrangement (described by Eq. (1)) in which both layers maintain the rigid arrangement expected of individual layers. Similar to experimental observations [16], the analytical [17,18] and ab initio [19,20] models show an expansion of the regions of the mSL in which the preferred graphene-hBN stacking is obtained. This is particularly pronounced at small misalignment angles (θ 1 • ). Similar effects have also been observed in molecular dynamics simulations [146], with comparison to the experiment [16] suggesting that the carbon-nitrogen interaction is two or three times stronger that the carbon-boron interaction. The breaking of the sublattice symmetry in the commensurate phase opens a band gap at zero-energy [18], which in ab initio models is enhanced significantly by the electronelectron interaction to a value of 0 ∼ 20 − 30 meV [19,20]. Also, Refs. [18][19][20] report that far lower values of 0 are obtained if the lattice deformations are artificially suppressed (note that, even in the total absence of deformations, a small but non-zero band gap can still be generated similar to Eq. (7)). Finally, it should be mentioned that the strain-induced vector and scalar potentials in Hamiltonian (4) produce a significant non-trivial contribution to the mSL potential [18,114,147]. Figure 3 provides a simple example of how a small strain, w ∼ √ δ 2 + θ 2 , (either in the graphene layer or the substrate) can induce qualitative changes in both the geometry of the mSL and its corresponding miniband structure. The top left panel shows the usual hexagonal mSL of an unstrained and perfectly aligned (θ = 0) heterostructure, while in the lower left panel a uniaxial strain w = δ has been applied to compensate the lattice mismatch along the x-axis so that the mSL becomes onedimensional. This strain is incorporated into the Hamiltonian (4) using u Gr

Homogeneously strained graphene in heterostructure with hBN
where we have taken w = δ and used σ = 0.165 [148] for . The corresponding band structure is then characterised by a wavevector k = (k x , k y ) with |k y | < |G|/2 and k x unbounded. This is shown, for a finite range of k x in the right panel of Fig. 3, and does not contain any of the sDPs characteristic of the bandstructure of the unstrained heterostructure with a similar mSL perturbation (compare with Fig. 2(d)). Interestingly, the translational symmetry of the strained one-dimensional mSL depends very sensitively on the direction of the applied strain, so that the corresponding spectral support develops a fractal structure as a function of the angle at which the strain is applied [82].

Topological aspects of the graphene/hBN heterostructure
One fundamental property of graphene is the Berry phase [149] of ±π acquired by an electron after traversing one orbit around the Dirac point [150][151][152]. For a general graphene/hBN heterostructure with broken inversion symmetry and gapped Dirac spectrum (at either the primary or secondary Dirac point), this phase is smeared out into a non-zero Berry curvature (k) over a range of energies of the order of the band gap [18,42]. When an external electric field is applied, this Berry curvature generates an 'anomalous' contribution to the electron's velocity in the direction perpendicular to the applied field [43], with u n (k) = e −ik·r ψ n,k (r) the cell-periodic part of the eigenstate ψ n,k (r), and n (k) energy of that eigenstate. In transport experiments [44], this is manifested as a transverse current, yielding a Hall-like conductivity for B = 0 [43,153], with f (k) the occupancy factor and g s = 2 the spin degeneracy. Due to the time-reversal symmetry, the transverse current has opposite signs for the two graphene valleys. However, intervalley scattering of Dirac electrons in graphene/hBN heterostructures is weak. In the transport experiment [44], the effect was observed as a nonlocal voltage in a narrow energy range near primary and secondary Dirac points at finite distances from the direct current path.
An alternative line of enquiry is motivated by the fact that, the Dirac Hamiltonian subject to a spatial varying mass term has topologically protected zero energy modes wherever the mass changes sign [154,155]. For the graphene/hBN heterostructure this has been applied to models [136,156] in which the sublattice-asymmetric mSL potential is treated as the dominant perturbation. In this model the conducting/insulating behaviour at zeroenergy is controlled by the classical percolation [157] of the zero energy modes (which is achieved when the spatial average of the mass vanishes). It has also been predicted [158] that the electron-electron interaction leads to a spontaneous spin-valley ordering for these zeroenergy modes.

Optical spectroscopy
In contrast to transport measurements, optical absorption spectroscopy can probe the strongly reconstructed spectrum at the edge of the first miniband of the graphene/hBN heterostructure without (in principle) the necessity of first shifting the chemical potential to such energies. The optical absorption coefficient is affected by both the changes to the energy spectrum, and the changes to the wavefunction of graphene electrons induced by the hBN underlay [37,39]. Similar to the studies on twisted bilayer graphene [123,128], this leads to a modulation of graphene's otherwise universal optical absorption coefficient of g 1 = πe 2 / c ≈ 2.3% [40], which will be most pronounced in the spectral range around ω ∼ vb. For much lower photon frequencies, the electron states are almost the same as in the unperturbed Dirac spectrum, whereas photons of much higher energies induce transitions between numerous overlapping minibands so that individual spectral features will be indistinguishable. Several examples of absorption spectra of an unstrained undoped graphene/hBN heterostructure, calculated [39] for realisations of the moiré perturbation in which only one of the three inversion-symmetric terms is present, are displayed in Fig. 4. Each absorption spectrum contains a characteristic pattern of peaks originating from various transitions between the minibands indicated with arrows on the corresponding bandstructure ( Fig. 4(b-d)).
The experimental measurement of the optical absorption spectra of graphene/hBN heterostructures is compounded by a large absorption background produced by the hBN underlay. To overcome this difficulty, Shi and co-workers [37] produced an electrostatically gated heterostructure and, by comparing absorption spectra for different graphene chemical potentials, could subtract the background contribution. In-terestingly, they found the measured absorption spectra consistent with that calculated using the mSL perturbation set by Eq. (9) which is based on the simple microscopic models. The right panel of Fig. 5 shows the theoretical absorption spectra, calculated using Eq. (9) with V + = 0.63vb for several values of hole doping. The left panel shows the miniband structure with coloured lines used to indicate the chemical potentials used in the right panel. By comparing the absorption spectra calculated for the undoped heterostructure (red line) to that in a slightly doped heterostructure (green line), it is seen that the effect of a small shift in chemical potential (|μ| < vb/2) consist of the usual Pauli blocking of transitions with ω < 2μ. This is the same as the effect on unperturbed graphene. In contrast, for higher hole doping (blue line, chemical potential tuned to the valence band sDP), new absorption channels can be opened. This results in the additional absorption peak marked with an asterisk in the right panel.

Fractal spectrum of magnetic minibands of the graphene/hBN heterostructure in a strong magnetic field
The fractal spectrum of electron waves constrained to two dimensions and exposed to both an in-plane periodic potential and a strong perpendicular magnetic field [24][25][26] is one of the most spectacular results in the quantum theory of solids [159]. When the magnetic flux, φ ≡ BS = p q φ 0 threading each unit cell (area S) of the system becomes a rational fraction ( p, q integers) of the flux quantum, φ 0 = h e , each Landau level is fractured into p magnetic bands. Since the width of the magnetic bands only becomes appreciable when φ 0.2φ 0 , their experimental observation at sustainable magnetic fields requires the use of systems with a large lattice constant (e.g. A ∼ 10 nm).
Early attempts to observe the fractal spectra focused on two-dimensional electrons in periodically patterned GaAs/AlGaAS heterostructures [27][28][29][30][31][32][33][34], where the superimposed superlattice period was made sufficiently large. However, the technological difficulties in producing such devices without introducing substantial disorder, and accessing the electron gas buried relatively deep within the sample, resulted in only partial success. In contrast, the graphene/hBN heterostructure, with the high quality of the exfoliated flakes, its naturally occurring mSL, and the ease of doping of graphene by electrostatic means, provides the perfect system to observe the fractal magnetic bands.
Already, several observations of the fractal bands have been reported in transport measurements, both for encapsulated [13] and non-encapsulated [14] heterostructures, as well as for devices made of bilayer graphene on aligned hBN [23]. Here, the signatures of the magnetic bands are the strongest when φ/φ 0 = 1/q [13,14,23]. Upon doping with 4 holes per moiré unit cell (corresponding to emptying the first B = 0 miniband), the Hall resistivity is observed to change sign for every φ/φ 0 = 1/q, indicating the recurrent generation of electron-like orbits in graphene's valence band due to the magnetic bands [13]. Also, transport measurements [14], and capacitance measurements [12] have revealed the lifting of spin and valley degeneracies of the magnetic minibands as a result of the electron-electron interaction. There have also been several theoretical works on the magnetic bands of graphene/hBN heterostructure [35,36,118,160,161], in addition to a plethora of works on the magnetic bands in the related twisted bilayer system [124][125][126][127][128].

Magnetic translational symmetry
To include a magnetic field in Hamiltonian (5) it is convenient [35] to use a Landau gauge for the magnetic vector potential, , written using a coordinate system adapted to the hexagonal symmetry of the mSL, r = x 1

A2
A . Then Hamiltonian (5) commutes with the magnetic translations [25,26,159] in the group, where T X is a geometrical translations on the mSL, and X X = e i2π p q m 1 m2 X +X , The subgroup of G M made of translations R = m 1 q A 1 + m 2 q A 2 on a (q × q)-enlarged superlattice is isomorphic to the simple group of translations, T R , so that its eigenstates, R | n, j t,k = e ik·R | n, j t,k , form a plane wave basis with a magnetic Brillouin zone with area q 2 -times smaller area than the mSL Brillouin zone. Moreover, since the non-abelian group G M has q-dimensional irreducible representations [26], the spectrum of such planewave states is q-fold degenerate.
To calculate the magnetic minibands, we project the Hamiltonian (5) on to the basis of Bloch functions, | n, j t,k , built from the wave functions of the unperturbed Landau levels, ψ k2 n (r), Here, β = B/|B|, Since the above basis has the convenient property, it is sufficient to fix k and consider t = 0 only. Also, the | n, j t,k diagonalize the Dirac part of Eq. (5) giving the energy of an unperturbed Landau level, E n = sign(n)vλ −1 B √ 2|n|, so that the mSL potentials can be treated perturbatively and the resulting magnetic bandstructure converges when only a finite range of n is included in the basis.

Generic features of the magnetic miniband spectra
Magnetic miniband spectra calculated for various mSL perturbations are displayed in Fig. 6. In each case the magnetic miniband spectra for weak magnetic fields, φ 0.2φ 0 , can be traced to a weakly broadened sequence of Landau levels associated with both the primary Dirac point as well as the sDPs. In higher magnetic fields, each Landau level is fractured to reveal the fractal magnetic miniband spectrum.
One feature of the magnetic minibands is their tendency to form a weakly gapped Dirac-like spectra. An example of this for φ/φ 0 = 1/2 is shown in the right panel of the inset in Fig. 2(d). The lower portion of this inset shows a colour map of the Berry curvature [43,149] for the lower magnetic miniband. This is concentrated at the tip of the gapped Dirac point, manifesting the topological character of such features. The left panel of the inset shows the magnetic miniband spectra (black points), at small deviations of the magnetic field from φ = φ 0 /2. This is grouped around [85] the Landau levels (blue lines) of the Dirac-like magnetic bandstructure calculated for φ = φ 0 /2 (albeit with the additional q-fold degeneracy due to the magnetic translational symmetry, and a finite gradient of the zeroth Landau level due to the finite magnetic momentum [43] of the magnetic bands). Similar features, identified in all panels of Fig. 6, are particularly pronounced in Fig. 6(c). Here the electron-hole symmetry H (−r) = −H (r), present in Hamiltonian (5) for the perturbation U + 3 = 0.15vb, prescribes a gapless Dirac point at zero energy whenever p is even [162] (clearly visible for e.g. φ/φ 0 = 2/1). A separate study by Diez and co-workers [36] also reveals that the generically gapped Dirac-like features, can become gapless upon varying a single mSL parameter (e.g. the misalignment angle θ ).
The remaining panels of Fig. 6 exemplify other interesting features in the magnetic miniband spectra. The principle novelty in Fig. 6(b), calculated using the mSL perturbation U + 1 = 0.15vb, is that it contains a completely unperturbed "n = 0" Landau level traced to the main Dirac point. This is a consequence of the electronhole symmetry, σ z H (r)σ z = −H (r), which is present for this particular mSL perturbation [162]. Fig. 6(a) is calculated for the mSL perturbation U + 0 = 0.15vb. Here the zero-field spectra contains a triplet of anisotropic sDPs on the edge of the first valence miniband, in contrast to Figs. 6(b-d) which display only a single isotropic sDP. Consequently, the three "n = 0" Landau levels traced to these three sDPs are hybridised by the mSL perturbation so that they undergo a magnetic breakdown at a lower magnetic field (φ ∼ 0.1φ 0 ) than similar features in the other three plots (φ ∼ 0.25φ 0 ). Finally, the spectra displayed in Fig. 6(d) was calculated for a mSL perturbation which contains a large inversion asymmetric component. In this case, the combination of spatial and timeinversion asymmetry lifts the valley degeneracy [35,118], so that the spectral supports for the two valleys must be shown separately (black for K and red for K ).

Interplay between magnetic bands and the electron-electron interaction
It is known from graphene devices manufactured on non-aligned hBN [8,163,164], suspended devices [165,166], and SiO 2 [167][168][169], that the electron-electron interaction can lift the four-fold spin-valley degeneracy of each Landau levels. This effect, known as quantum Hall ferromagnetism (QHFM) [170], is usually manifested through the appearance of gaps at all integer filling factors, ν = ρφ 0 /B (ρ the carrier density), and is particularly pronounced for ν = 0, ±1, corresponding to lifting the spin-valley degeneracy within the n = 0 Landau level. Related effects have also been observed for the magnetic minibands of the aligned graphene/hBN heterostructure, both in transport measurements [14] and capacitance measurements [12]. In the latter, incompressible states were traced to the band gaps created by the Landau level sequences of both the primary and secondary Dirac point. Other gaps were associated with the Landau quantization of magnetic minibands nearrational values of the flux. Further sequences of gaps could not be explained using the single-particle picture, and appeared due to the interplay [171] between the electron-electron interaction and the fractal magnetic miniband. Of these, the filling factors related to the QHFM expected in the absence of a mSL were suppressed for all rational fluxes φ = pφ 0 /q for which the band width of the magnetic minibands is large. This effect was found to be most pronounced around φ = φ 0 where the gaps at filling factors ν = ±1 (counted from charge neutrality) disappeared. Moreover, the Landau level sequences traced to small deviations in the magnetic field δ B from φ = φ 0 exhibited their own QHFM,

Moiré magnification of defects in graphene heterostructures
Section 2.4 discussed the effects of small homogeneous strains on the graphene/hBN heterostructures and its corresponding miniband structure. Here we use hBN as an example of a hexagonal substrate for graphene to discuss the geometrical magnification effect in which general extended or topological defects in either crystal layer generate a similar defect in the mSL, magnified by a factor M ∼ [δ 2 + θ 2 ] −1/2 [82].
To track the effect of these defects on the geometry of the mSL, we follow the coordinate dependent phase factor in Hamiltonian (4), so that we identify the formal "sites", of the mSL as points with a constant phase, g m · u(R) = 2π N. Then, comparing the location of these sites with (R d ) and without (R) the deformation due to a defect, we define the displacement field in the mSL by whereM is the "magnifying matrix" given in Eq. (1).
To demonstrate the use of Eq. (12), the left panel of Equation (13) can also be applied if the discontinuity in the mSL is caused by a step edge in the underlay [80]. In this case d would include both the effect of curving the graphene flake over the step edge, as well as the shift in stacking characterised by the different atomic layers in the underlay. Similarly, in the right panel in Fig 7,

Conclusion
The development of van der Waals heterostructures, created from stacks of two dimensional materials, is motivated by the possibility of creating devices in which the stacking, alignment, or interaction of the individual crystal layers adds new functionality to the device or opens a window on novel physical phenomena. For the graphene/hBN heterostructure the hBN underlay is not just an excellent substrate for graphene, contributing to the high electronic quality of these devices. Rather, it generates a large quasi-periodic moiré superlattice and provides a source of Bragg scattering for the Dirac electrons in graphene. This, in turn creates a specific miniband structure featuring gaps [14] and secondary Dirac points [11], which depends sensitively on the misalignment angle between the two crystal layers, and also appears to depend on whether the graphene is encapsulated or not with additional misaligned hBN layers [16]. The high quality of the two materials and the nanometerscale period of the moiré superlattice also enabled the first observations of the fractal magnetic miniband structure [13,23], and the subsequent lifting of the spin and valley degeneracies of the magnetic minibands by the electron-electron interaction [12,14]. Further studies have observed the modulation of graphene's otherwise flat optical absorption spectra [37]; the topological currents generated by regions of finite Berry curvature with in the minibands [44]; the generation of spontaneous strains in the graphene lattice [38]; and of course images of the moiré superlattice itself [9]. However all this is perhaps no more than just the tip of the iceberg, and we expect many more discoveries waiting for the graphene/hBN heterostructure. Moreover, there is now a growing library of two dimensional materials [3], and we hope the techniques developed in the study of this particular heterostructure will aid and inspire the development of many further devices created using this new diverse range of available materials.