Insulating state in tetralayers reveals an even–odd interaction effect in multilayer graphene

Close to charge neutrality, the electronic properties of graphene and its multilayers are sensitive to electron–electron interactions. In bilayers, for instance, interactions are predicted to open a gap between valence and conduction bands, turning the system into an insulator. In mono and (Bernal-stacked) trilayers, which remain conducting at low temperature, interactions do not have equally drastic consequences. It is expected that interaction effects become weaker for thicker multilayers, whose behaviour should converge to that of graphite. Here we show that this expectation does not correspond to reality by revealing the occurrence of an insulating state close to charge neutrality in Bernal-stacked tetralayer graphene. The phenomenology—incompatible with the behaviour expected from the single-particle band structure—resembles that observed in bilayers, but the insulating state in tetralayers is visible at higher temperature. We explain our findings, and the systematic even–odd effect of interactions in Bernal-stacked layers of different thickness that emerges from experiments, in terms of a generalization of the interaction-driven, symmetry-broken states proposed for bilayers.

C lose to charge neutrality and at zero magnetic field, B ¼ 0, electron-electron interactions in graphene and its multilayers have a strong effect in a very narrow range of electron density 1 . In monolayer graphene, Coulomb interactions renormalize the Fermi velocity of the Dirac fermions, which-in agreement with theoretical predictions 1 -increases as the Fermi level approaches the charge neutrality point (CNP) 2 . In Bernal bilayers, interactions are considered to be responsible for the opening of a gap 1,3-7 , and they turn the system into an insulator at a low temperature [8][9][10][11] . Thicker Bernal-stacked multilayers have received less attention, and it is only known that Bernal trilayers remain conducting at low temperature (rhombohedral trilayerson the contrary-have again been found to open a large gap due to interactions) 12 . In all cases-that is, for mono and bilayersthe effects of interactions are visible experimentally only at charge densities below 2-3 Â 10 10 cm À 2 , and their observation requires the study of suspended devices of the highest quality 13,14 , since otherwise the magnitude of charge inhomogeneity is too large.
Here we realize very high-quality multiterminal suspended devices based on Bernal-stacked tetralayer graphene (4LG), in which charge inhomogeneity can be reduced to the level of the best mono and bilayer structures realized in the past. Thanks to their very high quality enabled by employing a recently developed fabrication process (see ref. 15 and Methods for details), these devices allow us to reveal experimentally the occurrence of a very pronounced and robust insulating state at the CNP in 4LG, clearly visible already at relatively high temperature (at T ¼ 4.2 K the resistance at the CNP is several hundreds of kiloOhms, 50 times larger than in analogous bilayer devices). We find that the size of the energy gap in 4LG coincides with the values reported previously for bilayers, strongly suggesting a common origin for the insulating state in the two system. We discuss the phenomenon theoretically, and suggest that a symmetry-broken state in which a staggered potential is present in 4LGcorresponding to the simplest generalization of the symmetrybroken states that are considered to explain the behaviour of bilayers-naturally explains our observations, as well as the coincidence between the values of energy gap observed in 4LG and bilayers. These same theoretical considerations point to an even-odd interaction effect in Bernal-stacked graphene multilayers at the CNP (odd multilayers remain conducting at low temperature, whereas even multilayers become insulating) that is consistent with all experimental observations so far, and has implication for the crossover from graphene to graphite.

Results
Highly resistive state in charge neutral Bernal-stacked 4LG. Before discussing the detailed characterization of the thickness and the stacking of the layers in our devices, we present the experimental indications of the high device quality and of the unusual highly resistive state at the CNP, which can be appreciated already in the most basic transport measurements. Figure 1a shows the gate voltage (V G ) dependence of the twoterminal resistance (R 2-4 ) of two different devices measured at 4.2 K. In both cases, a very high and narrow peak is observed around the CNP, B50 times higher than that observed in suspended monolayer and Bernal-stacked bi/trilayer graphene at the same temperature 2,8-16 (B0.3 MO as compared with 4-8 kO; see Supplementary Fig. 1; the two devices exhibit identical behaviour and in the rest of the manuscript we show the data from one of them). The peak width is extracted from the log-log plot of the conductance G 2-4 ¼ 1/R 2-4 versus carrier density n, as illustrated in Fig. 1b (n is calculated from the applied V G , with the gate capacitance obtained from the analysis of the quantum Hall effect; see discussion of Fig. 2d). The range in which the conductance stays constant upon increasing n-a measure of the carrier density fluctuations 13,14 -extends only up to n * B2-3 Â 10 9 cm À 2 , comparable to the best values reported in the literature for suspended graphene of any thickness 2,[8][9][10][11][12][13][14][15][16] . Finally, the multiterminal configuration allows magnetotransport to be measured in both longitudinal and transverse configurations 15,16 . Already at moderately low B, a fully developed integer quantum Hall effect (IQHE) is observed (see Fig. 1c), with integer conductance plateaus (in units of e 2 /h) in the transverse resistance (R xy ) and concomitant vanishing of the longitudinal resistance (R xx ).
The thickness and stacking order of our samples are determined by magneto-Raman spectroscopy. The 2D Raman peak, occurring around 2,710 cm À 1 at B ¼ 0 when measured with an excitation at l ¼ 514 nm, exhibits the contribution of two components (see inset of Fig. 2a), as expected for 4LG with Bernal stacking [17][18][19] . This observation, however, is not sufficient to n* S2 S1 S2 Figure 1 | Very high resistance at charge neutrality in high-quality suspended tetralayer graphene devices. (a) V G dependence of the twoterminal resistance R 2-4 measured in two devices (S1 and S2) at 4.2 K, showing a peak resistance B0.3 MO at the CNP, much larger than the value normally observed in mono or BLG at the same temperature (R a-b represents the two-terminal resistance measured between contacts a and b, with the indexes labelling the contacts shown in the left inset). Left inset: optical microscope image of a device covered by a PMMA mask before etching (see Methods; scale bar, 2 mm long). Right inset: crystal structure of Bernal-stacked 4LG. (b) Conductance G 2-4 ¼ 1/R 2-4 as a function of n in a double-logarithmic scale, showing a very low level of charge inhomogeneity, n * ¼ 2-3 Â 10 9 cm À 2 (pointed to by the arrow). The density is given by n ¼ a(V G -V CNP ), with a ¼ 4.66 Â 10 9 cm À 2 V À 1 determined from the analysis of QHE (see Fig. 2d). (c), Fully developed QHE measured at B ¼ 0.45 T, with vanishing longitudinal resistance R xx ¼ R 1-4,2-3 (black curve) and quantized transverse resistance R xy ¼ R 1-3,4-2 ¼ 1/n Â h/e 2 (magenta curve) at n ¼ 8, 12 and 16 (data taken at 250 mK; R a À b,g À d corresponds to the ratio of the voltage difference measured between contacts g and d, and the current flowing from contact a to b).
exclude tri-and pentalayer graphene unambiguously. To gain conclusive evidence, Raman spectra were collected in the regime of Landau quantization, by applying a finite B to split the singleparticle energy spectrum into discrete Landau levels (LLs). In this regime, the G Raman peak exhibits characteristic anti-crossings upon increasing B-an oscillatory shift of the peak position with a concomitant broadening (see Fig. 2a,b)-when the energy of the E 2g phonon matches that of specific inter-LL transitions (see ref. 20 for details). In few layer graphene, whose electronic spectrum consists of several bands [21][22][23][24][25][26][27] , the effect of coupling is enhanced in the presence of nearly degenerate transitions (that is, when inter-LL transitions from different bands occur at nearly the same energy) 28 , and become measurable even at room temperature (for non-degenerate transitions, the effect of coupling is too weak to be detected). In the simplest approximation-that we will refer to as the effective bilayer model-the single-particle band structure of 4LG consists of two independent bilayer-like bands with different effective masses [21][22][23][24][25][26][27] , and predicts pairs of nearly degenerate inter-LL transitions to match the phonon energy at BB5.8, 7.8 and 12 T as marked by the blue circles between Fig. 2a-b 28 . As indicated by the blue dashed lines in Fig. 2a-b, the observed features are in virtually perfect agreement with these predictions, whereas they do not match the expectations for tri-and pentalayer graphene (respectively, magenta and cyan circles between Fig. 2a,b) nor those for non-Bernal-stacked multilayers (note, in this regard, that since our devices are realized by exfoliating natural graphite, the probability to find two-out-of-two identically non-Bernalstacked 4LG is entirely negligible). Magneto-Raman spectroscopy, therefore, unambiguously identifies our devices as made of Bernal-stacked 4LG. In addition, by showing the persistence of Landau quantization at room temperature, these observation further provide another clear indication of the high quality of our devices.
To characterize the device quality in more detail, and to gain additional understanding of the low-energy electronic properties, we have investigated transport in the quantum Hall regime by taking advantage of the multiterminal geometry that enables both R xx and R xy to be measured 15,16 . Figure 2c shows a colour plot of R xx (V G , B), exhibiting clear features originating from the QHE, with Shubnikov-de Haas oscillations becoming visible from BB0.1 T. Through the condition for their visibility (mBc1; refs 13,14), we estimate the carrier mobility m to be larger than 100,000 cm 2 V À 1 s À 1 . The oscillations evolve into fully developed IQHE states-with plateaus in R xy and vanishing R xx -starting from BB0.3 T at filling factor nnh/eB ¼ 8 (see the white dashed line in Fig. 2c), and subsequently at n ¼ 12 and 16 (starting from BB0.4-0.5 T). These states are clearly apparent in Fig. 2d, which shows the longitudinal (s xx ) and transverse (s xy ) conductivity measured as a function of V G for B in the range 0.475-0.8 T, and . At low n, below 2.5 Â 10 10 cm À 2 , the device becomes highly insulating (see Fig. 3), preventing multiterminal measurements. (d) Traces of the longitudinal and transverse conductivities, s xx ¼ r xx /(r xx 2 þ r xy 2 ) and s xy ¼ r xy /(r xx 2 þ r xy 2 ), as a function of n, measured for B between 0.475 T and 0.8 T (r xx ¼ R xx Â W/L and r xy ¼ R xy with the device aspect ratio W/L ¼ 0.87). The filling factor n ¼ nh/eB is calculated from B and n ¼ a(V G -V CNP ) with a ¼ 4.66 Â 10 9 cm À 2 V À 1 determined by imposing scaling of all the measured curves in this magnetic field-range. The good scaling enables averaging of the magnetoconductivity curves to suppress noise, as shown in e.
NATURE COMMUNICATIONS | DOI: 10.1038/ncomms7419 ARTICLE plotted as a function of n, with the gate capacitance relating V G and n determined by enforcing the scaling of all curves (note that the good scaling allows us to reduce noise by averaging 16 as shown in Fig. 2e; it also indicates the absence of significant inhomogeneity in the carrier density close to the edges of the 4LG, as it may be expected from basic electrostatic considerations 29 ). These fully developed states at low B follow the sequence expected from the simplest description of the electronic structure of 4LG in terms of two decoupled bilayer-like bands, for which a 16-fold degenerate LL at zero energy needs to be completely filled before electrons occupy (fourfold degenerate) LLs at higher energy 25 . The appearance of the n ¼ 4 state, which starts from B ¼ 0.8 T, indicates that, at higher fields, the 16-fold degeneracy of the zero energy LLs starts to break. Possible mechanisms responsible for breaking this degeneracy are the effect of next-nearest-layer hopping 30 (that is not considered in the effective bilayer model) or electron-electron interactions 31 , for which an analysis of the behaviour at larger B provides additional evidence (see discussions below). Irrespective of these details, the observation of the 16-fold degenerate zero-energy LL provides useful information about the low-energy electronic properties of suspended 4LG.
Insulating state and electron-electron interaction in 4LG. With the device structure and quality characterized well, we look at the high resistance peak found close to the CNP (Fig. 1a). At base temperature (T ¼ 250 mK; see Methods), the conductance G becomes unmeasurably small (the resistance is 4100 MO), as shown in Fig. 3a. The different curves in the figure correspond to two-terminal measurements done using different pairs of contacts (as indicated in the legend; for these high resistance values, fourterminal measurements cannot be done reliably), and show that the strong suppression always occurs in a same range of carrier density |n|o2-3 Â 10 10 cm À 2 , irrespective of the measurement configuration. This finding is an indication of the very high device homogeneity, which can be checked directly when working with multiterminal devices 15,16 . Figure 3b shows the temperature dependence of G(V G ) measured in one of the two-terminal configurations (other configurations give identical results), showing a pronounced insulating behaviour upon lowering T. Above 1 K, the minimum conductance is thermally activated with an activation energy E A B15 K (see Fig. 3c), indicating the presence of an energy gap at the CNP. The gap closes very rapidly when the device is gate-biased away from the CNP, as it can be seen in the measurement of the differential conductance (dI/dV) as a function of different V G and applied bias voltage V B (Fig. 3d). Figure 3e further shows that at low bias, the dI/dV vanishes, within the accuracy of the measurements, for all bias voltages V B below a threshold (B1.5 mV, corresponding approximately to the activation energy extracted from the Tdependent measurements), above which it increases sharply. In the case of bilayer graphene (BLG), a qualitatively similar behaviour has been observed 8,9 (see also Supplementary Note 1) and interpreted in spectroscopic terms, with the applied bias taken to be a direct measure of the energy of the injected electrons 8,9 . It cannot be excluded, however, that the dominant effect of the applied bias is the electrostatic accumulation of charges on the suspended 4LG flake. Indeed, since the gap closes very rapidly upon increasing the charge density away from the CNP, the accumulation of even a very small amount of charges can result in a drastic suppression of the gap, leading to a large and sharp conductance increase.
The occurrence of an insulating state at the CNP is unexpected. It cannot be accounted for in terms of the gapless single-particle band structure of 4LG [21][22][23][24][25][26][27] , and its explanation calls for the effects of electron-electron interactions (a scenario based on the presence of macroscopic defects such as AB-BA stacking faultswhich has been proposed for bilayers 32,33 to explain why some of high-quality suspended devices remain conducting at low T (ref. 34)-is incompatible with our data, and with different observations reported in the literature; Supplementary Note 2). Indeed, the appearance of fully developed QHE at B ¼ 2.3 T with plateaus in s xy at 1, 2 and 3 e 2 /h and vanishing s xx (Fig. 4a) demonstrates that the interactions are strong enough to completely lift the degeneracy of the zero-energy LL already at low B. Plotting dG/dB as a function of B and n shows that these symmetry-broken states start developing from BB1 T (Fig. 4b).
This value corresponds to the magnetic field at which the n ¼ 4 state becomes visible, suggesting that the appearance of this state is also due to interactions. Note that, in the fractional quantum Hall effect (likely to be within experimental reach at higher B with devices of the quality shown here), the 16-fold degeneracy of the zero-energy LL is expected to lead to new unexplored regimes, occurring when states belonging to different LLs are mixed by interactions [35][36][37] . The first manifestations of these new regimesthe appearance of an even-denominator fractional state at n ¼ À 1/2 (ref. 16) and of electron-hole asymmetry 38,39 -have just been reported in BLG: with its larger degeneracy, 4LG should lead to an even richer spectrum of new phenomena. Finally, additional evidence for the relevance of the interactions in 4LG is provided by the enhancement of the insulating state around charge neutrality, which at higher magnetic field extends throughout a larger density range (Fig. 4c). This behaviour is identical to that observed in mono and BLG, where it has been shown to originate from a canted anti-ferromagnetic state due to interactions [40][41][42] (a similar explanation may be valid for 4LG. Indeed, if we consider the magnitude of the gap B1.5 meV, the data in Fig. 4d, which show that the insulating state is not affected by a parallel magnetic field B up to 15 T, indicates that the electron spin is not a good quantum number in the insulating state). The finding that the enhancement starts already at low B (Fig. 4c) and that the insulating state at finite field evolves continuously into the B ¼ 0 insulating state upon reducing B, suggests a related origin of both the states.
Gap origin and even-odd effect in graphene multilayers. A comprehensive analysis of the effect of electron-electron interactions in Bernal-stacked 4LG is complex and goes beyond the scope of this paper. However, our observation that the magneto-Raman measurements and the IQHE at low B are in surprisingly good agreement with the behaviour expected from the effective bilayer model suggests that treating 4LG by analogy with the case of graphene bilayers is a good starting point. At a mean-field level 3,5,7,8 , the low-energy properties of the symmetry-broken states in BLG are commonly described within a two-band model, in terms of a 2 Â 2 matrix Hamiltonian of the type for each valley and spin (the '±' sign changes going from the K to the K' valley). The first term is the kinetic energy of the electrons (m * is the effective mass) and the second one represents the order parameter of the broken symmetry state (s x,y,z are Pauli matrices acting on the layer degree of freedom). Different ground states emerge, depending on how the sign of the order parameter D changes upon switching valley or spin 7,8 , but, in all cases, the 'bulk' of BLG is gapped (see Fig. 5a) as the term Ds z in the Hamiltonian has an opposite sign in opposite layers for any fixed spin and valley configuration. The analogy with BLG suggests that this mean-field treatment can be extended to the case of 4LG in a natural way, by generalizing the interlayer potential asymmetry in BLG as a staggered potential V i ¼ ( À 1) i þ 1 D (i ¼ 1,y,4 is the layer index) which changes sign between the neighbouring layers. More specifically, the Hamiltonian of 4LG is decomposed into blocks describing two effective graphene bilayers using a unitary transformation, and the staggered potential then works as the Ds z term for each bilayer (see Supplementary Note 3). It is straightforward to show by direct diagonalization of the 4LG Hamiltonian that, within the effective bilayer model, the staggered potential D opens a gap at the CNP in 4LG (see Fig. 5a; again, in analogy to the BLG, ground states with different properties shown in Supplementary Fig. 2  above discussion based on the analogy with BLG certainly cannot substitute a complete theoretical analysis-and more complex possibilities cannot be ruled out at this stage-it has the benefit of showing explicitly a realistic scenario based on electron-electron interactions, which leads to the opening of a gap at the CNP and naturally explains why the gap size in 4LG is essentially identical to that observed in BLG. Interestingly, the scenario that we propose also captures the systematic behaviour of charge neutral Bernal multilayers, which emerges experimentally with our finding of an interaction-driven insulating state in 4LG. Specifically, as it has been reported in previous experiments, 'odd' (that is, mono and tri) layers remain conducting at low temperature 2,12-14 , whereas, as we can conclude from this work, 'even' (that is, bi [8][9][10][11] and tetra) layers exhibit a gap and become insulating. This phenomenon cannot simply be explained in terms of an enhanced susceptibility to interactions due to a larger low-energy density of states, since Bernal-stacked trilayers, for instance, have a larger density of states than bilayers, and yet they remain conducting at low temperature. On the contrary, describing the effect of interactions in terms of a staggered mean-field potential at the level of approximation used in this study (see Equation (1) and the Supplementary Note 3) leads precisely to the even-odd effect that is observed experimentally. In fact, at the same level of approximation used to discuss 4LG, the electronic structure of even multilayers generally decouples into bilayer-like bands [21][22][23][24][25][26][27] , and since each band opens a gap in the presence of a staggered potential, the electronic state of a generic even multilayer is gapped (see Fig. 5a). When performing the same analysis for odd multilayers, next to the bilayer-like bands, a linearly dispersing monolayer-like band emerges and remains ungapped in the presence of a staggered potential (see Fig. 5b). As a result, odd multilayers remain metallic without opening a gap (see Supplementary Note 3 for more details). It remains to be understood why the simple approximation used in this study to describe the Bernal multilayers (which neglects, for instance, next-layer hopping processes, responsible for the g 2 and g 5 parameters in the tight-binding description of graphite) correctly captures the basic aspects of the electron-electron interactions. The robustness of the even-odd effect may be rooted in the symmetry of these systems, which is different in the two cases: even multilayers are inversion symmetric, whereas odd ones possess reflection symmetry 27 (for thickness larger than the monolayer). It therefore seems plausible that the appearance of a staggered potential-which breaks inversion and preserves reflection symmetry-can account for a drastic modification of the behaviour of even multilayers, and not of that of odd ones.

Discussion
The experimental observation of the even-odd effect just discussed is particularly interesting because it shows how the interactions induce robust, systematic behaviour in graphene multilayers, irrespective of the specific values of the tight-binding hopping parameters. It is worth noting the similarity between the effect discussed here and the phenomenology observed in a class of very different materials, namely ladders of antiferromagnetically coupled spin ½ chains 43 . There, the spectrum of spin excitations also exhibits an even-odd effect, being gapped for ladders composed by an even number of chains and gapless when the number of chains is odd. Theoretical concepts that have been considered to explain the even-odd effect observed in spin ladders include topological terms related to the non-trivial Berry phase of spin ½ and the presence of a Dirac dispersion (ubiquitous in truly one-dimensional systems), which clearly have analogies with the physics of graphene 44 . Nevertheless, there are important differences, for example, in spin ladders the evenodd effect is seen when going from 1D to 2D systems, whereas with graphene multilayers we go from 2D to 3D.
An important implication of the robustness of the even-odd interaction effect, and of the experimental findings reported here, is that electron-electron interactions can be expected to determine also the properties of multilayers that are so thick to be normally considered as bulk graphite. Indeed, finding an insulating state in 4LG already clearly visible well above 4.2 K, that is, more pronounced than in BLG (where it becomes visible only below 1 K; see Supplementary Note 1), shows that the role of interactions does not become weaker, upon increasing layer thickness. This finding and the theoretical considerations made above, therefore, indicate that the commonly held notion that bulk graphite can be well understood within a single-particle model is unlikely to be correct, and that at a sufficiently lowenergy the effects of electron-electron interactions become important. These effects can easily be eclipsed in experiments on samples of microscopic size, which unavoidably contain structural defects, but when probing well-characterized, small multilayers-which can be made defect free-they become fully apparent.

Methods
Device fabrication and measurement. Samples in this study were fabricated by using polydimethylglutarimide (PMGI)-based lift-off resist (LOR, MicroChem) as a sacrificial layer 15,45 . Graphene flakes were mechanically exfoliated from natural graphite onto an B1-mm-thick LOR layer covering the highly doped Si/SiO 2 substrate (acting as a back gate), and the metallic electrodes (10 nm Ti/70 nm Au) contacting the selected flakes were defined by using a standard microfabrication process: electron-beam lithography, deposition and lift-off. Before removing the sacrificial LOR resist locally under the flake to achieve suspension, the flake shape was defined by the exposure to an Oxygen plasma through a patterned polymethylmethacrylate resist mask (see the inset of Fig. 1a for an optical image of etching pattern). After suspension, the devices were mounted in the vacuum chamber of a Heliox 3 He system (Oxford Instruments; base temperature of 250 mK) and current annealed at low temperature (T ¼ 4. configuration, by using shorted adjacent pair of contacts (1 and 2) and (3 and 4) as a source and drain, respectively 15 . As we adopt the strategy to proceed with the current annealing until all obvious manifestations of inhomogeneity in the transport properties of the device disappear, multiple annealing steps are required and the number of devices that survive the process in the end is small (in other words, we do not start with in-depth investigations of transport, which can easily take several months for a single device, as long as clear non-idealities are present in the transport characteristics of the devices).