Unconventional satellite resistance peaks in moir\'e superlattice of h-BN/ AB-stacked tetralayer-graphene heterostructure

Most studies on moir\'e superlattices formed from a stack of h-BN and graphene have focused on single layer graphene; graphene with multiple layers is less understood. Here, we show that a moir\'e superlattice of multilayer graphene shows new features arising from the anisotropic Fermi surface affected by the superlattice structure. The moir\'e superlattice of a h-BN/AB-stacked tetralayer graphene heterostructure exhibited resistivity peaks showing a complicated dependence on the perpendicular electric field. The peaks were not due to secondary Dirac cones forming, but rather opening of the energy gap due to folding of the anisotropic Fermi surface. In addition, superlattice peaks resulted from mixing of light- and heavy-mass bilayer-like bands via the superlattice potential. The gaps did not open on the boundary of the superlattice Brillouin zone, but rather opened inside it, which reflected the anisotropy of the Fermi surface of multilayer graphene.


Introduction
Moiré is a kind of geometrical interference that appears when two nearly identical lattices are aligned with a small angle mismatch ( Fig. 1(a)). The discovery of moiré in a stack of h-BN and graphene [1][2] has led to numerous studies associated with superstructures of graphene formed by incommensurate substrates [3][4][5][6][7][8][9][10][11][12]. The moiré structure virtually gives graphene another periodic structure, and therefore, these systems can be regarded as superlattices [13][14]. The conventional understanding of superlattices is that an energy gap forms at the boundary of the Brillouin zone of the superlattice. However, recent findings in monolayer graphene have shown that Dirac cones are formed by locally closing the energy gap, and this results in additional peaks in the carrier density dependence of resistance [3][4][5][6]10]. Further studies have indicated that the secondary Dirac cones are formed at the K or K' point of the superlattice Brillouin zone [15][16][17][18]. The moiré potential forms a band gap of about 20 meV at the secondary Dirac cone in the hole branch [7,12]. In addition, the superlattice potential forms a band gap of about∼ 35 meV at the primary Dirac cone in the case of perfect alignment, which appears at the charge neutrality point [4][5][6][7]12,[19][20]. These band-gap values are highly dependent on the alignment angle and sample structure [4,12].
Although there have been a number of studies on moiré superlattices of h-BN/single layer graphene heterostructures, there have only few studies on moiré superlattices of graphene with multiple layers. As far as the authors know, bilayer graphene has been studied in terms of Hofstadter butterfly's diagram in the energy spectrum obtained in a magnetic field [9,18,21]. The superlattice potential due to the moiré would primarily influence the surface layer of graphene in contact with it. Therefore, it is not clear how the electronic band structure is modified by the moiré potential in graphene with multiple layers. Moreover, trigonal warping of the band structure becomes increasingly dominant as the number of layers increases [22][23]. This effect would be important in forming an energy gap due to the moiré structure, and it would also be important in reconstructing the Fermi surface. In this work, we studied the electronic band structure of a moiré superlattice of AB-stacked tetralayer graphene, as an example of multilayer graphene, by conducting low-temperature transport experiments. We found new features intrinsic to the moiré superlattice of multilayer graphene.
The band structure of AB-stacked tetralayer graphene consists of two bilayer-like bands (light and heavy mass) as shown in Fig. 1(b). A simplified band model indicates that the tetralayer graphene should have a semi-metallic band in which there is an overlap of the electron and hole bands near = 0; however, mixing between the heavy and light-mass bilayer-like bands creates local energy gaps, and thus, the band overlap of the conduction and valence bands is expected to be considerably smaller than predicted by the simplified band model [24][25][26][27]. The band overlap vanishes by applying a perpendicular electric field; an energy gap opens at the bottoms of the bilayer-like bands [25][26] as in the case of bilayer graphene [28][29][30][31][32][33][34][35][36]. The bottoms of the heavy-mass and light-mass bilayer-like bands become nearly flat [25][26][27]. However, trigonal warping locally closes the energy gap, resulting in mini-Dirac cones at the bottom of the heavy-mass bilayer-like band.
These characteristic structures in the dispersion relation were revealed by studying the intrinsic resistance peaks that as a function of carrier density and perpendicular electric flux density [25][26][37][38].

Results and Discussion
Figure 1 (c) shows an optical micrograph of an AB-stacked tetralayer graphene sample that consisted of encapsulated graphene with a top and a bottom gate electrode. The lowtemperature electrical mobility of this device was above 4 × 10 4 cm −2 /Vs. As shown in The resistivity traces ( vs ) showed significant variation with top gate voltage ( Fig.   2(a)), and detailed measurements of the resistivity map revealed a fine structure of resistivity peaks ( Fig. 2(b)). The straight ridge from the top left to bottom right of Fig.   2(b) is the condition of charge neutrality, which is given by = 0, where is the carrier density induced by the top and bottom gate voltage, Here, and are the specific capacitances for the top and bottom gate electrodes.
The perpendicular electric flux density ( ⊥ ), which is roughly proportional to the strength of the perpendicular electric field, can be calculated as Figure 3(a) is a replot of the map. One can see ridges of resistivity (bb+, bb-, MDP) that were recently found in AB-stacked tetralayer graphene [25][26][27][37][38]. Ridges bb+ and bb-reflect the bottoms of the bilayer-like bands, whose dispersions are nearly flat in a perpendicular electric field. MDP is due to mini-Dirac cones. However, the other ridges (A+, B+, C+, etc.) are absent from pristine AB-stacked tetralayer graphene.
The Landau levels were also influenced by the moiré structure. Figure 1 (e) shows a map of longitudinal resistivity ( ), which was measured as a function of and magnetic field (B ) for ⊥ = 0 cm −2 As. Each streak indicated by bright and dark lines is a Landau level and the energy gap between them. The Landau level structure between = −2 × 10 12 and 2 × 10 12 cm -2 are specific to pristine AB-stacked tetralayer graphene [25][26][27]39].
In addition, at ∼ 3.5 × 10 12 cm −2 , one can see a Landau fan reminiscent of that for the secondary Dirac cones in the case of moiré superlattice in bilayer graphene. Around this carrier density, ridge A+ appears (Fig. 2(c)). Such Landau fans could not be clearly seen in the hole regime, which implies that the moiré superlattice for multilayer graphene is different from those in the mono-or bilayer case.
To study the origin of the ridges, we compared the results of the experiment with a theoretical calculation. Figure 2 (d) shows a map of resistivity that was numerically calculated from the dispersion relations by using Boltzmann transport theory with the constant relaxation time approximation. The band structure calculation took account of an effective moiré potential based on the h-BN-graphene hopping model presented by Moon and Koshino [18]. The amplitude of the effective potential and mismatch angle of the crystal axes between h-BN and graphene were adjusted so that the calculation  The detailed structure of the numerically calculated dispersion relation is highly dependent on the models of the effective moiré potential and its amplitude. However, the essential feature of the Fermi surface topology is unchanged by the choice of model. We calculated the dispersion relations and resistivity maps for h-BN-graphene hopping models [16][17][18]42], the 2D charge modulation model [43], and the potential modulation model [3,16,42] (see the Supplementary material). The h-BN-graphene hopping models approximately reproduce the experimental resistivity map if the potential amplitude is less than about half of those given in Refs. [16][17][18]42]. Using any of these models, the resistivity map was approximately reproduced for a sufficiently small potential amplitude, which indicates that the topological transitions of the Fermi surface due to nesting occur regardless of the model used.
If the resistivity peaks (ridges) originate from nesting of the Fermi surface, the resistivity maps should show a significant dependence of the trigonal warping because the nesting is significantly influenced by the shape of the Fermi surface. Figure  Finally, we comment on the period of the superlattice potential and mismatch angle between the crystallographic axes of h-BN and that of graphene. In the case of the moiré superlattice of h-BN/monolayer graphene, one can calculate the period of the moiré potential [3][4][5][6] from the size of the superlattice Brillouin zone, which can be estimated from the carrier density of the resistivity peak due to the superlattice ( ) [3][4][5][6]. The mismatch angle can be estimated as [3][4][5][6], Here, is the wave length (or lattice constant) of the moiré potential, is the difference between the lattice constants of h-BN and graphene, and is the lattice constant of graphene [3]. In the multiband system, however, does not reflect the size of the superlattice       Namiki, Tsukuba 305-0044, Japan, 3 International Center for Materials Nanoarchitectonics, National Institute for Materials Science, 1-1 Namiki, Tsukuba 305-0044, Japan.

A Sample fabrication and method of measurement
The schematic structure of our sample is depicted in Figure S1 (a). The sample was h-BN-encapsulated AB-stacked tetralayer graphene formed on Si substrate covered with SiO2. The h-BN and graphene flakes were prepared by mechanical cleaving using adhesive tape. The stack was formed using dry transfer process with polypropylene carbonate (PPC) film [1]. The sample was equipped with top and bottom gate electrodes.
The top one was made of a few-layer graphene and the bottom one was the Si substrate which was heavily doped and remained conducting down to liquid helium temperatures.
The sample was patterned into a Hall bar via reactive ion etching using a mixture of low-pressure CF4 and O2 gas. The electric contacts to the graphene were formed by using the edge-contact-technique [1].
The number of graphene layers was roughly determined by making optical contrast measurements [2][3][4]. Then, the number of layers and stacking were verified via a Landau fan diagram, which is a map of resistivity measured as a function of carrier density and magnetic field. The structure of the Landau levels characterizes the electronic band structures of graphene. It is different for different numbers of layers and stackings, and accordingly, the fan diagrams exhibit specific patterns. By comparing the diagram with those in the past studies [5], the present sample was identified to be AB-stacked tetralayer graphene.
Resistivity measurements were performed by using the standard four-probe method and lock-in technique. The magnetic field was applied using superconducting solenoid.
B Effect of perpendicular electric field on moiré superlattice of h-BN/bilayer graphene heterostructure.
The dependence of resistivity on the top and bottom gate voltage was studied in a moiré superlattice of a bilayer h-BN graphene hetero structure to compare with the case of ABstacked tetralayer graphene. Pristine bilayer graphene has one massive band at low energy, which contrasts with the AB-stacked tetralayer case with two massive bands.  (17)- (18) in Ref. [7] was added to the diagonal element of the first layer of graphene that contacts h-BN.
The main text presents the results calculated using the effective moiré Hamiltonian of Moon and Koshino [7], which is given by (S1) Here, = ( On the other hand, Wallbank has derived an effective moiré potential for h-BN-graphene heterostructure [8,9]; Here, = ( 1 , 2 , 3 ), and ( = 1, 2, 3) are the Pauli matrices. is a primitive vector of graphene, which is defined by rotating ( , 0) by (2 /6). is a reciprocal vector of the moiré superlattice. In ref. [9] values of ( 0 + , 1 + , 3 + ) are given for different models of effective moiré potential. In case of the potential modulation model, with 0 ∼ 60 meV [9][10]. In the 2D charge modulation model [8,9,11], the parameters are given by 0 is a phenomenological parameter. We used 0 = 60 meV in the calculation. In the one-site version of the G-hBN hopping model [12], it is approximately given by The numerical calculation used 0 = 30 meV for a small mismatch angle .
To compare with the numerically calculated band structures with the experimental results, we used a potential amplitude scaled by a factor s, which served as an adjustable parameter. We also adjusted the mismatch angle between the crystal axes of h-BN and graphene.
The dispersion relations were calculated by using the Hamiltonian of the continuum model [6], Here, 2ℏ . Here, is and ( = 1 ∼ 4) is the electrostatic potential due to the top and bottom gate voltage.
reflects the perpendicular electric field so that one must consider screening of the external electric field properly in multilayer graphene [13][14][15][16]. We calculated by assuming that charges induced by each gate voltage attenuate with a relaxation length of ∼ 0.45 nm. [16][17][18][19].
Resistivity was numerically calculated by using the constant relaxation time approximation of the Boltzmann transport theory. In this calculation the derivative of the Fermi function with respect to energy was approximated using a Gaussian function with a width of 1 meV. D Maps of resistivity for different effective potential models and Fermi surface topologies Figures S3-S6 show maps of resistivity numerically calculated for different models of effective moiré potential. Figure S3 is for the h-BN-graphene hopping model of Moon and Koshino [7], Fig. S4 is for the one-site version of the G-hBN hopping model [9,12], S5 is for the 2D charge modulation model [9,11], and S6 is for the potential modulation model [10]. In each figure, the factor s was varied from left to right as 0.25, 0.5 and 1, with the right-most panel showing the experimental results. It can be seen that the resistance ridge structures are strongly dependent on s, which reflects the variation in the detailed structure in the dispersion relations. The maps are slightly asymmetric in ⊥ = 0. This is because the effective moiré potential is present only in the first layer of graphene in the Hamiltonian, which results in an effective offset in ⊥ . In all the potential models, the map for = 1 shows a significant difference from the experiment.     Numerically calculated maps of resistivity by using the effective moiré potential of the 2D charge modulation model [8,9,11]. s is the scaling factor of the potential. = 0.35 ∘ .
The SWMcC parameters of graphite were used in the calculation.