# The Unusual Milky Way-Local Sheet System:
Implications for Spin Strength and Alignment

M.A. Aragon-Calvo1, Joseph Silk2, Mark Neyrinck3,4,5
1Instituto de Astronomía, UNAM, Apdo. Postal 106, Ensenada 22800, B.C., México
2Institut d’Astrophysique de Paris, UMR7095:CNRS UPMC-Sorbonne University,
F-75014, Paris, France
3 Ikerbasque, Basque Foundation for Science, 48011 Bilbao, Spain
4 Department of Physics, University of the Basque Country UPV/EHU, 48080
Bilbao, Spain
5 DIPC, Basque Country UPV/EHU, 20018 San Sebastián, Spain

###### Abstract

The Milky Way and the Local Sheet form a peculiar galaxy system in terms of
the unusually low velocity dispersion in our neighbourhood and the seemingly
high mass of the Milky Way for such an environment. Using the TNG300
simulation we searched for Milky Way analogues (MWA) located in cosmological
walls with velocity dispersion in their local Hubble flow similar to the one
observed around our galaxy. We find that MWAs in Local-Sheet analogues are
rare, with one per (160-200 Mpc)3 volume. We find that a Sheet-like cold
environment preserves, amplifies, or simplifies environmental effects on the
angular momentum of galaxies. In such sheets, there are particularly strong
alignments between the sheet and galaxy spins; also, these galaxies have low
spin parameters. These both may relate to a lack of mergers since wall
formation. We hope our results will bring awareness of the atypical nature of
the Milky Way-Local Sheet system. Wrongly extrapolating local observations
without a full consideration of the effect of our cosmic environment can lead
to a Copernican bias in understanding the formation and evolution of the Milky
Way and the nearby Universe.

###### keywords:

Cosmology: large-scale structure of Universe; galaxies: kinematics and
dynamics, Local Group; methods: data analysis, N-body simulations

††pagerange: The Unusual Milky Way-Local Sheet System:  Implications for Spin
Strength and Alignment–LABEL:lastpage††pubyear: 2002

##  1 Introduction

The Milky Way (MW) is an SBb/c galaxy with an estimated mass between
$8.3\times 10^{11}$M⊙ and $2\times 10^{12}$ M⊙ (Karukes et al., 2020). It is
located near the center of a flat sheet of galaxies extending up to 4 Mpc
along along the supergalactic plane (Kroupa, Theis & Boily, 2005; Metz, Kroupa
& Jerjen, 2007; Pawlowski, Kroupa & Jerjen, 2013) with a density close to the
mean density of the Universe (Klypin et al., 2015) and delineated by a ring of
luminous galaxies, the “Council of Giants” (McCall, 2014). The “Local Sheet”
in which our galaxy is located belongs to a class of cosmic structures known
as cosmic walls. It is bounded to the supergalactic north and south by the
Local Void and the Southern Void respectively (Tully et al., 2008). There are
no galaxies in the line of sight between the our galaxy and the center of both
north and southern voids. The Milky Way is sitting at the edge of a cosmic
cliff.

An increasing body of evidence points to the Milky Way and its immediate
cosmic environment as somewhat uncommon occurrences in the Universe. The
velocity dispersion measured for galaxies located within 3 Mpc from our own
galaxy is $\sigma_{v}\sim 25-40$ km s-1(Karachentsev et al., 2002, 2003), one
order of magnitude lower than predictions (Peebles, 1980; Governato et al.,
1997) and observed pairwise velocity dispersions (Landy, 2002; Zehavi et al.,
2002; Scoccimarro, 2004). The “coldness” of the local Hubble flow seems to be
a global property of the Local Sheet as indicated by observations around the
nearby Centaurus A / M83 system, showing a dispersion of $\sigma_{v}\sim 32$
km s-1 (Karachentsev et al., 2007). The coldness of the local Hubble flow can
be explained as the result of the characteristic geometry and dynamics of
cosmic walls (Aragon-Calvo, Silk & Szalay, 2011). The central regions of walls
are less dense than their surrounding denser “ring” of filaments, producing a
local expansion on top of the global Hubble flow which damps peculiar
velocities.

The halo mass function in walls also indicates the peculiarity of Milky Way
Analogue (MWA) galaxies as it is roughly one order of magnitude lower than the
global mass function in the MW mass range (Aragón-Calvo (2007); Ganeshaiah
Veena et al. (2018), see also Neuzil, Mansfield & Kravtsov (2020)). The rarity
of MWA galaxies is evident, albeit not explicitly addressed, in studies of MWA
in which constraints on halo mass, local density and geometry are included and
also in constrained realizations of the Local Group. In both cases a large
volume or a large number of realizations are needed in order to find MWA
candidates (see for instance Garrison-Kimmel et al. (2014) and Fattahi et al.
(2016)).

Walls are dynamically young structures, being the result of the gravitational
collapse along one axis (Zel’dovich, 1970). In the tidal torque theory (Hoyle,
1949; Peebles, 1969; White, 1984), the tidal field governs spin. In a
collapsing wall environment, this tidal field squashes protohaloes into a
plane. Any aspherical lumps on the protohalo will get caught in the tidal
field, imparting torque to the halo, with its spin tending to end up
perpendicular to the axis of collapse (see also Neyrinck et al., 2020). This
alignment between galaxy spin axes and their parent wall planes is small but
statistically significant in dark matter computer simulations (Aragón-Calvo et
al., 2007b; Hahn et al., 2007). Locally, the alignment of edge-on spiral
galaxies (including our own) perpendicular to the Local Sheet is so strong
that it can be seen by eye. As shown in Navarro, Abadi & Steinmetz (2004),
gas/stars can retain their primordial orientation better than dark matter, due
to the latter being more affected by mergers. The strength of the alignment
signal in our neighborhood may be the result of disk galaxies keeping their
original spin alignment and a quiet mass growth and low velocity dispersion in
the relatively low density of the Local Sheet. On the other hand, mergers are
a well known mechanism for angular momentum acquisition (Vitvitska et al.,
2002). In the case of galaxy interactions within walls, mergers would tend to
induce a spin aligned with the normal of the wall, not with the plane, as
prescribed by the TTT.

Figure 1: Hubble diagram measured from a MWA in our sample. The dots
correspond to background galaxies inside the search radius, the solid line
indicates the best fit to the model given in sec. 2.1. The peculiar velocity
dispersion is indicated by the dashed lines. The vertical dotted line marks
the radius of zero velocity. We also show the unperturbed Hubble expansion
$v=Hr$ as the dashed gray line. The mass and fitting parameters are shown in
the figure. Gray dots correspond to real galaxies presented in Karachentsev et
al. (2003).

##  2 N-body Simulation and LSS analysis

The analysis presented here is based on Illustris TNG300 (Nelson et al.,
2018), a cosmological hydrodynamical simulation of galaxy formation enclosed
in a box of 205$h^{-1}$Mpc of side containing dark matter and gas particles
with masses of $5.9\times 10^{7}$ and $1.1\times 10^{7}$ $h^{-1}$M⊙
respectively. The cosmological parameters used are given by Planck results
(Planck Collaboration et al., 2016). While TNG300 has a relatively low mass
resolution it has the large volume needed to sample rare cold Local Sheet
environments as we will see in the following sections.

We defined two set of haloes: background haloes used to sample the velocity
field and Milky Way Analogues. The background halo sample consisted of subfind
haloes in the mass range $1\times 10^{9}<$M$<5\times 10^{12}$ M⊙ and
$V_{\textrm{\tiny max}}>20$ km s-1. MWAs were defined as central haloes with
$V_{\textrm{\tiny max}}>20$ km s-1 in the mass range $10^{12}<$ M $<2\times
10^{12}$ M⊙. We extracted mass accretion histories for the MWA sample and
identified merger events as increases in $30\%$ of the total mass of the halo
between snapthots. The MWA sample was classified into spiral and ellipticals
using a custom convolutional neural network feed with synthetic images created
from the stellar component of the halos observed in the $r,g,i$ SDSS filters
(Aragon-Calvo, in preparation, see for instance Rodriguez-Gomez et al.
(2019)). For every MWA halo we computed the spin parameter following the
definition of Bullock et al. (2001):
$\lambda=|J|/(\sqrt{2}M_{200}R_{200}V_{200})$, where $J$ is the angular
momentum, $M_{200}$ and $R_{200}$ are the virial mass and radius respectively
and $V_{200}$ the circular velocity.

The density field used for cosmic web analysis was computed by running a
resimulation of TNG300 sampled on a regular grid with $512^{3}$ particles. The
initial conditions were sharp-$k$-smoothed in order to target structures
larger than 4 h-1Mpc as described in Aragon-Calvo et al. (2010) and evolved to
$z=0$ using the GADGET-2 N-body code (Springel et al., 2005). From the final
snapshot we labeled wall regions using the Hessian-based MMF method (Aragón-
Calvo et al., 2007a; Aragon-Calvo & Yang, 2014).

###  2.1 Measuring the Hubble flow

The Hubble flow around MWAs was measured by first selecting background haloes
inside a shell of $1<r<4$ Mpc centered in the target MWA galaxy. This distance
range was chosen to roughly follow the range used in real measurements in the
local Universe (Karachentsev et al., 2009) although extending to a larger
radius in order to increase the number of available galaxies. Only MWAs with
more than 10 neighbours were included in the analysis. We computed velocity
dispersions $\sigma_{v}$ from the model-subtracted radial velocities centered
on the MWA using the model (Peirani & de Freitas Pacheco, 2008):
$v(r)=H_{0}r-H_{0}R_{0}(R_{0}/r)^{1/2}$, where $r$ is the distance from the
MWA and $H_{0}$ and $R_{0}$ are the fitted Hubble constant and radius of the
zero-velocity surface respectively (see Fig. 1).

Figure 2: Hubble flow excess ($H_{0}-H$) as a function of $\sigma_{v}$ around
the MWA sample. The vertical dotted lines indicate the range of observed
velocity dispersions around our galaxy (see text for details). The 2D
distribution is shown as a gray background with inverted color scale. The
black cross shows the median of $\sigma_{v}$ and $H-H_{0}$ The contours
enclose the areas containing 86,95 and 99.7 percent of galaxies. The top
subplot shows the cumulative distribution of $\sigma_{v}$.

##  3 Results

Figure 2 shows the distribution of $\sigma_{v}$ and Hubble flow excess (H0-H)
measured from the MWA sample. The distribution of both $\sigma_{v}$ and H0-H
have a roughly log-normal character. The peak of the velocity dispersion is
surprisingly low, around $\sigma_{v}\sim 50$km s-1. However our MWA sample is
restricted to relatively massive galaxies and does not include galaxies in
pairs, groups or clusters which can have a much higher $\sigma_{v}$. The tail
of the distribution, being $\sim$log-normal, extends to much higher values in
$\sigma_{v}$ and the mean of the distribution is $178$ km s-1. The peak in the
Hubble excess distribution is located at H0-H $\sim$ 10 km s-1. The tail of
the distribution extends into the negative hundreds (not shown in Fig. 2 for
clarity). The median value of H0-H is -2.2 km s-1, close to the global
expansion. The reported range in the velocity dispersion $25<\sigma_{v}<40$
corresponds to $3\sigma-1\sigma$ events for H0-H=0. Interestingly, from Fig. 2
we should expect to measure a Hubble flow excess around our galaxy of $\sim
H+(10-15)$ km s-1 Mpc-1.

Figure 3: Mass function in different cosmic environments defined by their
$\sigma_{v}$ and LSS morphology. We show the mass function of haloes in cold
($\sigma_{v}<25$ km s-1) and warm ($\sigma_{v}<40$ km s-1) environments. For
each $\sigma_{v}$ we also show the mass function of haloes in cosmic walls.
The vertical light shaded area shows the mass range of MWA and the dark grey
area highlights the most likely region where the MW would be located in the
diagram. For comparison we show the global mass function. Figure 4:
Distribution of halo spin parameter $\lambda$ for different MWA sub-samples.
Left panel: Distribution of $\lambda$ for spirals and ellipticals, the solid
line corresponds to the full MWA sample and the dashed line to MWA galaxies
with 3 or more major merger events since $z=2$. We also show synthetic images
(three orthogonal orientations) of one spiral and one elliptical galaxy.
Center panel: The distribution of $\lambda$ for spiral galaxies in cold
($\sigma_{v}<40$ km s-1) and increasingly hot environments ($\sigma_{v}>100$
km s-1 and $\sigma_{v}>200$ km s-1). The vertical arrow and error bar indicate
the most recent estimate of the spin parameter for the Milky Way
($\lambda=0.061^{+0.022}_{-0.016}$, Obreja, Buck & Macciò (2022) (OBM22)).
Right panel: Same as the center panel but for elliptical galaxies.

###  3.1 How common are MWAs in cold walls?

Figure 3 shows the effect of cosmic environment, defined in terms of geometry
(walls) and $\sigma_{v}$, on the halo mass function (see also Aragón-Calvo
(2007); Ganeshaiah Veena et al. (2018)). The mass function decreases in height
and slope with decreasing Hubble velocity dispersion. The effect of
environment is even stronger when we consider only haloes in walls. At low
masses the mass function in all environments is a scaled version of the global
mass function. The slope in the mass function begins to decrease with
increasing mass around $10^{12}$ M⊙. The “knee” of the curve associated with a
characteristic mass in cold environments is more than one order of magnitude
lower than the “knee” of the global population. The slope in the mass function
in cold wall environments is significantly steeper than the global population,
the effect being even more pronounced after the characteristic mass range. The
Milky Way is right at the border of the knee in the mass function in walls in
the two regimes of $\sigma_{v}$ considered here. The probability of drawing a
MWA in cold wall environments from the general population
$P(>M)=(N(>M)_{\textit{\tiny env}}/N(>M)_{\textit{\tiny global}}$ lies
somewhere between 0.2% and 0.001% for the ranges in mass and $\sigma_{v}$
considered here. This is equivalent to finding one MWA in a box of side $\sim
160$ and $\sim 200$ Mpc respectively. If we take the lowest estimate of
$\sigma_{v}\sim 25$ km s-1Milky Way analogues in cold Local Sheet environments
can be considered extremely rare events.

###  3.2 Galaxy/Halo Spin in Local Sheet

####  3.2.1 Spin parameter

Before addressing environmental effects, we study the spin parameter
distribution of galaxies separated only by their morphology. Figure 4 (a)
shows that spirals have in general a higher $\lambda$ compared to ellipticals
with the peaks in their $\lambda$ distribution, $P(z)$, at $\lambda\sim 0.03$
and $\lambda\sim 0.015$ respectively (see also Du et al. (2021); Rodriguez-
Gomez et al. (2022)). Note that we use a strict definition of morphology,
rejecting galaxies with $<90\%$ probability of being a spiral/elliptical given
by the neural network. Both spirals and ellipticals with 3 or more major-
merger events since z = 1 have a higher spin parameter, supporting the idea of
angular momentum acquisition via mergers (Vitvitska et al., 2002).

The effect of environment on $\lambda$ is shown in Fig. 4 (b). Both spirals
and ellipticals in walls have lower $\lambda$ in cold environments. The peak
in $P(\lambda)$ increases with $\sigma_{v}$ and the trend saturates for values
higher than $\sigma_{v}\sim 100$ km s-1. In all cases spirals have a higher
$\lambda$ compared to ellipticals in their respective environment.
Importantly, $P(\lambda)$ for galaxies in cold and hot environments not only
differ in the position of their peaks but also in their character. The
$P(\lambda)$ of galaxies in cold walls has a more narrow peak centered at a
lower $\lambda$ value and has a lower high-$\lambda$ tail indicating a
systematic lack of high-$\lambda$ galaxies in this environment. This effect is
more pronounced in the case of spirals.

####  3.2.2 Spin alignment

Finally, we show the alignment of galaxy spins with walls they inhabit. Figure
6 (right panel) shows that the alignment signal increases with decreasing
$\sigma_{v}$ and increasing halo mass. The alignment signal of haloes with
mass $<2\times 10^{12}$ is weak, consistent with previous findings (Aragón-
Calvo et al., 2007b; Hahn et al., 2007) and showing a weak dependence with
environment. In comparison, haloes above this mass present a significant
alignment signal and dependence with environment. The star/disk component of
galaxies is more strongly aligned than the dark matter component. Regarding
cosmic environment, galaxies in cold environments are more strongly aligned
than galaxies in hot environments. The strength in the alignment of the
stellar component of galaxies in cold walls, similar to our Local Sheet, is
very close to the observed value reported in Navarro, Abadi & Steinmetz
(2004). Interestingly, the Milky Way has a somewhat high $\lambda$ for its
environment (see arrow in Fig. 6 center panel). We note that the high
$\lambda$ galaxies are also the strongest aligned. However the sample size is
too small to draw a firm conclusion and report here. We will study this in a
future work.

Figure 5: Mean mass accretion histories for spiral and elliptical MWA galaxies
in increasingly hotter environments.

##  4 Conclusions and discussion

In this work we present evidence that the Milky Way, a common galaxy by mass
with respect to the global population, is an outlier when we consider the
environment in which it is located. Improving on previous works, we describe
the Local Sheet environment in terms not only of its geometry, but also its
dynamics. We find that a Local Sheet and Milky Way environment, consisting of
an unusually massive galaxy in a cold wall, is rare, occurring only once in a
160-200 Mpc box. We note that given the uncertainties on the mass of our
galaxy, the measured $\sigma_{v}$ and the general size and shape of the Local
Sheet, our results should only be taken as an approximate indication of the
rarity of such systems. Adding conditions to a situation will always
eventually give a vanishingly small probability so we adopt somewhat
conservative constraints and define the MW-Local Sheet system in terms of only
three parameters: LSS geometry, LSS dynamics and halo mass. Our results apply
only to Milky Way analogues, if we consider the more massive Milky Way -
Andromeda system, with a combined mass of $\sim 3-4\times 10^{12}$ M⊙, such
event is even rarer, so we do not analyze the joint situation of two large
galaxies in a cold wall. But note that after the two galaxies merge, the
situation will be the same as we analyze here, except with a larger combined
mass. A cold wall with a galaxy of this larger mass would be extremely rare; a
much larger box than the TNG300 volume we analyze would be necessary to obtain
an example.

We show that both spiral and elliptical galaxies located in cosmological walls
with low $\sigma_{v}$ have a lower spin parameter, characterized by a narrower
peak in the distribution located at lower $\lambda$ and a lack of
high-$\lambda$ galaxies. This could indicate a lack of mergers in the quiet
wall environments (Hellwing et al., 2021) (see also Fig. 4) emphasized by the
low $\sigma_{v}$. At the same time, we found a strong alignment signal between
the spin of spiral galaxies (selected using a neural network) and the plane of
their host wall. The value of the alignment strength is close to the observed
value around our galaxy. To our knowledge this is the strongest spin alignment
in walls reported in $N$-body simulations. This is significant since we
observe an unusually strong spin alignment of spiral galaxies around our
galaxy, yet no previous study to date has been able to reproduce it.

Our results highlight the importance of carefully characterizing the
environment around our galaxy. The effect of the geometry and coldness of
Local Sheet environment on angular momentum processes may help us better
understand current problems in galaxy formation such as the co-rotating plane
of satellites around M31 (Ibata et al., 2013), the missing satellite problem
(Klypin et al., 1999) among others.

Figure 6: Spin alignment of spiral galaxies in walls. We include galaxies with
95% probability of being spiral. The alignment $\cos\theta$ corresponds to the
angle between the spin vector of the galaxy/halo and the plane of its host
wall. We compute the alignment for the stars and dark matter components and
cold and hot environments. The gray area indicates the dispersion inside each
bin for the spin alignment of the star component of galaxies in cold
environments. The other cases have similar dispersion.

##  5 Data Availability

The data underlying this article were accessed from the Illustris TNG
datababse (https://www.tng-project.org). The derived data generated in this
research will be shared on reasonable request to the corresponding author.

##  6 Acknowledgements

This research was partially funded by the UNAM-PAPIIT grant IA102020. M.A
Aragon-Calvo would like to thank Dylan Nelson for kindly providing the TNG300
parameter file used for generating ICs.

## References

  * Aragón-Calvo (2007) Aragón-Calvo M. A., 2007, PhD thesis, University of Groningen, Netherlands
  * Aragón-Calvo et al. (2007a) Aragón-Calvo M. A., Jones B. J. T., van de Weygaert R., van der Hulst J. M., 2007a, A&A, 474, 315
  * Aragon-Calvo, Silk & Szalay (2011) Aragon-Calvo M. A., Silk J., Szalay A. S., 2011, MNRAS, 415, L16
  * Aragon-Calvo et al. (2010) Aragon-Calvo M. A., van de Weygaert R., Araya-Melo P. A., Platen E., Szalay A. S., 2010, MNRAS, 404, L89
  * Aragón-Calvo et al. (2007b) Aragón-Calvo M. A., van de Weygaert R., Jones B. J. T., van der Hulst J. M., 2007b, ApJLett, 655, L5
  * Aragon-Calvo & Yang (2014) Aragon-Calvo M. A., Yang L. F., 2014, MNRAS, 440, L46
  * Bullock et al. (2001) Bullock J. S., Dekel A., Kolatt T. S., Kravtsov A. V., Klypin A. A., Porciani C., Primack J. R., 2001, ApJ, 555, 240
  * Du et al. (2021) Du M., Ho L. C., Debattista V. P., Pillepich A., Nelson D., Hernquist L., Weinberger R., 2021, ApJ, 919, 135
  * Fattahi et al. (2016) Fattahi A. et al., 2016, MNRAS, 457, 844
  * Ganeshaiah Veena et al. (2018) Ganeshaiah Veena P., Cautun M., van de Weygaert R., Tempel E., Jones B. J. T., Rieder S., Frenk C. S., 2018, MNRAS, 481, 414
  * Garrison-Kimmel et al. (2014) Garrison-Kimmel S., Boylan-Kolchin M., Bullock J. S., Lee K., 2014, MNRAS, 438, 2578
  * Governato et al. (1997) Governato F., Moore B., Cen R., Stadel J., Lake G., Quinn T., 1997, New Astronomy, 2, 91
  * Hahn et al. (2007) Hahn O., Porciani C., Carollo C. M., Dekel A., 2007, MNRAS, 375, 489
  * Hellwing et al. (2021) Hellwing W. A., Cautun M., van de Weygaert R., Jones B. T., 2021, Phys. Rev. D, 103, 063517
  * Hoyle (1949) Hoyle F., 1949, Problems of Cosmical Aerodynamics, Proceedings of the Symposium on the Motion of Gaseous Masses of Cosmical Dimensions, Paris, 195?197.
  * Ibata et al. (2013) Ibata R. A. et al., 2013, Nature, 493, 62
  * Karachentsev et al. (2009) Karachentsev I. D., Kashibadze O. G., Makarov D. I., Tully R. B., 2009, MNRAS, 393, 1265
  * Karachentsev et al. (2003) Karachentsev I. D. et al., 2003, A&A, 398, 479
  * Karachentsev et al. (2002) Karachentsev I. D. et al., 2002, A&A, 389, 812
  * Karachentsev et al. (2007) Karachentsev I. D. et al., 2007, A&A, 133, 504
  * Karukes et al. (2020) Karukes E. V., Benito M., Iocco F., Trotta R., Geringer-Sameth A., 2020, JCAP, 2020, 033
  * Klypin et al. (2015) Klypin A., Karachentsev I., Makarov D., Nasonova O., 2015, MNRAS, 454, 1798
  * Klypin et al. (1999) Klypin A., Kravtsov A. V., Valenzuela O., Prada F., 1999, ApJ, 522, 82
  * Kroupa, Theis & Boily (2005) Kroupa P., Theis C., Boily C. M., 2005, A&A, 431, 517
  * Landy (2002) Landy S. D., 2002, ApJLett, 567, L1
  * McCall (2014) McCall M. L., 2014, MNRAS, 440, 405
  * Metz, Kroupa & Jerjen (2007) Metz M., Kroupa P., Jerjen H., 2007, MNRAS, 374, 1125
  * Navarro, Abadi & Steinmetz (2004) Navarro J. F., Abadi M. G., Steinmetz M., 2004, ApJLett, 613, L41
  * Nelson et al. (2018) Nelson D. et al., 2018, MNRAS, 475, 624
  * Neuzil, Mansfield & Kravtsov (2020) Neuzil M. K., Mansfield P., Kravtsov A. V., 2020, MNRAS, 494, 2600
  * Neyrinck et al. (2020) Neyrinck M., Aragon-Calvo M. A., Falck B., Szalay A. S., Wang J., 2020, The Open Journal of Astrophysics, 3, 3
  * Obreja, Buck & Macciò (2022) Obreja A., Buck T., Macciò A. V., 2022, A&A, 657, A15
  * Pawlowski, Kroupa & Jerjen (2013) Pawlowski M. S., Kroupa P., Jerjen H., 2013, MNRAS, 435, 1928
  * Peebles (1969) Peebles P. J. E., 1969, ApJ, 155, 393
  * Peebles (1980) Peebles P. J. E., 1980, The large-scale structure of the universe
  * Peirani & de Freitas Pacheco (2008) Peirani S., de Freitas Pacheco J. A., 2008, A&A, 488, 845
  * Planck Collaboration et al. (2016) Planck Collaboration et al., 2016, A&A, 594, A13
  * Rodriguez-Gomez et al. (2022) Rodriguez-Gomez V. et al., 2022, MNRAS, 512, 5978
  * Rodriguez-Gomez et al. (2019) Rodriguez-Gomez V. et al., 2019, MNRAS, 483, 4140
  * Scoccimarro (2004) Scoccimarro R., 2004, Phys. Rev. D, 70, 083007
  * Springel et al. (2005) Springel V. et al., 2005, Nature, 435, 629
  * Tully et al. (2008) Tully R. B., Shaya E. J., Karachentsev I. D., Courtois H. M., Kocevski D. D., Rizzi L., Peel A., 2008, ApJ, 676, 184
  * Vitvitska et al. (2002) Vitvitska M., Klypin A. A., Kravtsov A. V., Wechsler R. H., Primack J. R., Bullock J. S., 2002, ApJ, 581, 799
  * White (1984) White S. D. M., 1984, ApJ, 286, 38
  * Zehavi et al. (2002) Zehavi I. et al., 2002, ApJ, 571, 172
  * Zel’dovich (1970) Zel’dovich Y. B., 1970, A&A, 5, 84


# Magnetoconductance Oscillations in Electron-hole Hybridization Gaps and
Valley Splittings in Tetralayer Graphene

Illias Klanurak  Department of Physics, Faculty of Science, Chulalongkorn
University, Patumwan, Bangkok 10330, Thailand  Kenji Watanabe  Research Center
for Functional Materials, National Institute for Materials Science, 1-1
Namiki, Tsukuba 305-0044, Japan.  Takashi Taniguchi  International Center for
Materials Nanoarchitectonics, National Institute for Materials Science, 1-1
Namiki, Tsukuba 305-0044, Japan.  Sojiphong Chatraphorn  Department of
Physics, Faculty of Science, Chulalongkorn University, Patumwan, Bangkok
10330, Thailand  Thiti Taychatanapat<EMAIL_ADDRESS>Department of
Physics, Faculty of Science, Chulalongkorn University, Patumwan, Bangkok
10330, Thailand

###### Abstract

We investigate magnetotransport on Bernal-stacked tetralayer graphene whose
band structure consists of two massive subbands with different effective
masses. Under a finite displacement field, we observe valley splitting of
Landau levels (LLs) only in the light-mass subband, consistent with a tight-
binding model. At low density, we find unexpected magnetoconductance
oscillations in bulk gaps which originate from a series of hybridizations
between electron-like and hole-like LLs due to band inversion in tetralayer
graphene. In contrast to a trivial LL quantization gap, these inverted
hybridization gaps can lead to a change in number of edge states which
explains the observed oscillations.

The electronic properties of Bernal-stacked tetralayer graphene (4LG) have
been shown to exhibit many intriguing phenomena such as insulating state [1],
interlayer interaction [2], tunable Lifshitz transitions [3], helical edge
states [4], and unconventional satellite peaks [5]. Its band structure
features band inversion as a result of band overlap and hybridization between
the two bilayer-graphene-like (BLG-like) subbands [6, 7, 8]. Bernal-stacked
trilayer graphene (TLG) also possesses a band overlap but it occurs at a
narrow density range which is hard to resolve experimentally. In addition, the
subbands in TLG only hybridize in a presence of perpendicular electric field,
while those in 4LG always hybridize due to the next-nearest interlayer hopping
parameters [8].

The band inversion in 4LG provides a unique platform to study quantum Hall
effect (QHE) in a regime in which both electron-like and hole-like states
coexist. Previously, the energy inversion of electron-like and hole-like
Landau levels has been engineered by electric field or both in-plane and out-
of-plane magnetic fields to form helical edge states [9, 10, 4]. Some of two-
dimensional quantum wells also host inverted electron-hole system which, under
a magnetic field, reveals novel aspects of QHE [11, 12] and quantum spin Hall
effect [13, 14, 15, 16, 17, 18].

Figure 1: (a) Crystal structure of 4LG. The purple dashed lines indicate two
hopping sites for $\gamma_{i}$ parameters. (b) Band structure of 4LG from TB
model. The heavy- and light-mass bands are labeled as B (red) and b (blue),
respectively. The lattice constant $a$ is $0.246$ nm. (c) Left: an optical
image of an encapsulated 4LG device. Scale bar, $2$ $\mu$m. Right: schematic
diagram illustrating a cross-sectional view of the 4LG device. (d) Two-
terminal conductance as a function of $V_{\rm BG}$ at $0$ (black), $3$ (pink),
and $7$ T (green), respectively.

In this work, we investigate magnetotransport properties of 4LG encapsulated
by hexagonal boron nitride (hBN). The electronic band structure of 4LG can be
characterized by a set of hopping parameters $\gamma_{0}$-$\gamma_{5}$ [Fig.
$1$(a)] and energy imbalance between dimmer and non-dimmer sites $\delta$ [8].
The low energy band structure of 4LG, shown in Fig. $1$(b), comprises two BLG-
like subbands with different effective masses [19, 8, 20]. We denote the
light-mass subband by b and the heavy-mass subband by B. The presence of
skewed lattice sites hopping $\gamma_{3}$ induces trigonal warping. The next-
nearest interlayer hopping parameters $\gamma_{2}$ and $\gamma_{5}$ cause the
b and B subbands to overlap and hybridize at low energy as shown in Fig.
$1$(b) [8, 3]. The remaining hopping parameters $\gamma_{4}$ and $\delta$
generate electron-hole asymmetry in the band structure.

Figure 2: (a) conductance as a function of $n$ and $B$. (b) ${\rm d}g/{\rm
d}n$ numerically calculated from data in (a). The numbers at the right border
indicate filling factors. (c) LL spectra in 4LG calculated from TB model. The
filling factors associated with some energy gaps are labeled. The labels at
the right border denote quantum states for each LL. At high field, LLs from
the b and B subbands are in blue and red, respectively. At low field, LLs from
both subbands are in black due to LL hybridization. (d) Calculated density of
states of 4LG as a function of $n$ and $B$.

In our devices, we use a dry transfer technique to assemble the
heterostructures and etch them in CHF3/O2 plasma. Edge contacts are defined by
electron beam lithography and formed by sputtering 80-nm Mo [21]. The standard
lock-in technique is used to investigate magnetotransport properties. All
measurements are performed at 2.4 K unless stated otherwise. Figure $1$(c)
shows an optical image and a schematic diagram of a device.

We first examine two-terminal conductance $g$ as a function of back-gate
voltage $V_{\rm BG}$ at magnetic field $B=0$ [Fig. $1$(d), black line]. We
find that the conductance curve exhibits multiple local minimums, which are
associated with band edges and Lifshitz transitions near zero energy in the
band structure of 4LG as shown in Fig. $1$(b) [3, 22]. At finite $B$, we
observe conductance oscillations in $V_{\rm BG}$ due to the QHE [Fig. $1$(d)].
Our conductance does not develop into well-defined plateaus even at $7$ T,
likely due to relatively high contact resistance ($\sim$1 k$\Omega$) and
geometry effect on two-terminal conductance which leads to the distortion of
the QH plateaus [23, 24]. At high field, the conductance exhibits a single
minimum at which we associate to the charge neutrality point (CNP) of the
sample. Near-zero gate voltage for the CNP indicates pristine quality of our
samples.

To further investigate the magnetotransport properties of the system, we
measure $g$ as a function of $V_{\rm BG}$ and $B$ [Fig. $2$(a)]. Charge
carrier density $n$ is determined from the period of conductance oscillations
at high $B$. The Landau fan diagram displays rich features associated with LLs
from the b and B subbands. To see the oscillations more clearly, we calculate
${\rm d}g/{\rm d}n$ [Fig. $2$(b)]. The dark lines correspond to energy gaps
with associated filling factors $\nu$ shown in the figure. Multiple LL
crossings are evident. For instance, the dark line for $\nu=12$, seen clearly
at $7$ T, disappears at about $5$ T and reemerges again at $3$ T. The absence
of the dark line for $\nu=12$ between $3$ and $5$ T is the result of the LL
crossing.

To gain a better understanding of the Landau fan diagram, we calculate
energies of LLs as a function of $B$ using the tight-binding (TB) model [Fig.
$2$(c)]. In this plot, we assume that potential difference between layers is
zero and use $\gamma_{0}=3.1$, $\gamma_{1}=0.39$, $\gamma_{2}=-0.022$,
$\gamma_{3}=0.315$, $\gamma_{4}=0.12$, $\gamma_{5}=0.018$ and $\delta=0.020$
eV. These TB parameters are determined by matching the LL crossing positions
from the experiment at low density with those from the calculation [25]. We
note that our device has a single gate. As we induce higher carrier density
via back gate, 4LG is inevitably subject to a stronger displacement field
which induces larger potential difference between layers. Therefore, our
simulation gives a good agreement with data from low density at which
potential difference is still small.

From the spectra, LLs of 4LG at high $B$ can be viewed as a combination of two
sets of BLG-like LLs from the b and B subbands [8, 26]. At high energy, a LL
energy is approximately linear in $B$, as expected from bilayer nature of the
subbands. The energy spacing of the LLs from the light-mass band b is larger
than that of the heavy-mass band B because cyclotron frequency is inversely
proportional to effective mass. The mixing between LLs due to $\gamma_{2}$ and
$\gamma_{5}$ parameters and trigonal warping effect from $\gamma_{3}$
parameter leads to hybridization gaps, more visible at low energy [8]. We
label each LL in Fig. $2$(c) with two indices, indicating the subband (B or b)
and LL index $n$ ($n+$ for electron-like and $n-$ for hole-like LLs). For the
zeroth index LL, we label them as (b/B, 0) and (b/B, $-1$). These two LLs are
degenerate in bilayer graphene but the degeneracy is slightly lifted in 4LG
with a small energy gap of $\sim$2 meV at $7$ T. The numbers inside LL energy
gaps in Fig. $2$(c) indicate values of $\nu$ associated with the gaps. Each LL
has four-fold spin and valley degeneracy [8].

To compare our data with the calculation, we simulate density of states as a
function of $B$ and $n$ from the LL spectra in Fig. $2$(c) to obtain the plot
in Fig. $2$(d). Here, we assume a Lorentzian line-shape for each LL with
broadening of 1.5 meV, estimated from a LL gap from the b subband at the onset
of the oscillations of $0.5$ T [Fig. $2$(b)]. The simulation captures main
features of the experimental data in Fig. $2$(b) such as the positions of LL
crossing at low density and a position of the horizontal line in the electron
side which originates from the zeroth LL of the light-mass band.

However, some discrepancies exist between our data and simulation. In our
data, some LLs are two-fold degenerate. For example, in Fig. $2$(b) at $6$ T,
we observe the filling factor sequences of {$8$, $10$, $12$} and {$24$, $26$,
$32$} which imply LLs with degeneracy of $2$ and $6$. Comparing the positions
of these LLs with the LL spectra in Fig. $2$(c), we conclude that each of the
4-fold degenerate (b, $0$) and (b, $1+$) split into $2$ LLs with $2$-fold
degeneracy while the LLs from the B subband remain 4-fold degenerate. The
$6$-fold degeneracy observed when filling factors change from $26$ to $32$ is
the result of LL crossing between the $2$-fold degenerate LL from the
splitting of (b, $1+$) and the $4$-fold degenerate (B, $5+$).

Figure 3: (a) LLs energy of the bulk 4LG as a function of $B$ with potential
difference between adjacent layers of $16$ meV. The black and red lines are
LLs for $K$ and $K^{\prime}$, respectively. (b) Calculated density of states
of the LLs in (a) as a function of $n$ and $B$.

In 4LG, a 2-fold degenerate LL can occur by applying a displacement field to
generate a potential difference between layers which breaks inversion symmetry
and lifts valley degeneracy. To capture the effect of the potential difference
on LLs, we simulate LL spectra using a constant value of the potential
difference between adjacent layers of $16$ meV and the same set of the hopping
parameters used in Fig. $2$(c). We note that, in our measurement, a value of
potential difference will vary as we change density (see more details in
Supplemental Material [27]). The black and red lines in Fig. $3$(a) represent
LLs from $K$ and $K^{\prime}$ valleys, respectively. The valley splittings of
the LLs from the B subband are much smaller than LL broadening and therefore
they continue to appear $4$-fold degenerate in our measurement. However, the
splittings are much more pronounced for LLs from the b subband, consistent
with the data in which the LL splittings are observed in the b subband only.
Although we can explain most observed features within a single-electron
picture, we cannot completely rule out interaction-induced LL splitting. For
example, the (b, $0$) and (b, $-1$) are so close in energy that they should
experience similar value of the potential difference. However, we observe the
splitting of (b, $0$) but not (b, $-1$) which may suggest that other symmetry-
breaking mechanisms are involved.

Figure 4: (a) Color map of ${\rm d}g/{\rm d}n$ at low $n$ and $B$. The integer
numbers indicate filling factors of the dark diagonal lines. (b) Conductance
oscillations at $n=0$ as a function of $B$ at various temperatures. The
temperature at which the oscillations disappear ($\sim$20 K) is consistent
with the size of the $\nu=0$ gap. (c) The bulk LL spectra of 4LG (black
lines). The blue and red lines are LL spectra of b and B subbands,
respectively, calculated without band hybridization. (d) Diagrams of LL edge
states at $B_{c1}^{<}$ and $B_{c1}^{>}$, indicated by pink and orange lines in
(c). We omit the zeroth LLs (B, $0$) and (B, $-1$) (two red horizontal lines).

A more surprising discrepancy between our data and the simulation occurs at
zero filling factor. Figure 4(a) shows a plot of ${\rm d}g/{\rm d}n$ at low
$n$ and $B$ while Figure 4(b) displays $g$ at $n=0$ as a function of $B$ at
various temperatures. We observe oscillations in magnetoconductance clearly
along the zero density in all three devices we have measured (see data in
Supplemental Material [27]). Typically, the conductance oscillations in QHE
occur when Fermi energy passes through a LL. However, from the calculation of
the LL spectra from the TB model [Fig. 4(c), black lines], there is no LL
crossing inside $\nu=0$. As a result, the conductance should exhibit no
oscillation in $B$ at $\nu=0$, contradicting our results.

To resolve this discrepancy, we examine the complex nature of low-energy LL
spectra. Due to the band inversion in 4LG at zero magnetic field [see Fig.
1(b)], hole states from the b subband reside at higher energy than electron
states from the B subband for wave vectors around zero. As $B$ increases,
energies of hole-like LLs decrease while those of electron-like LLs increase.
The opposite magnetic dependence of LL energies leads to a series of crossings
and anti-crossings which manifests as three adjacent LLs braided together at
low field [see Fig. 4(c), black lines]. These anti-crossings, whose energy
gaps depend on $\gamma_{2}$ and $\gamma_{5}$ parameters, are the result of
hybridization between hole-like LLs and electron-like LLs from the b and B
subbands, respectively. A braiding of three LLs is the consequence of the
trigonal warping effect from $\gamma_{3}$ parameter. The effect causes the
anti-crossings to occur when the LL indices of the unperturbed LLs are the
same or differ by multiples of $3$ [28, 29, 30, 31].

To identify the underlying LLs that hybridize into the braided LLs, we
calculate LL spectra by setting the mixing terms between the subbands to zero
[8]. The result is displayed in Fig. 4(c). The blue and red lines represent
LLs of the b and B subbands, respectively. Comparing the LLs with and without
the mixing terms, we find that the energy gap at $\nu=0$ originates from two
different mechanisms with a crossover at the critical field $B_{c1}$ of
$\sim$2.1 T. For $B>B_{c1}$, the $\nu=0$ gap is a trivial LL gap which is
always present regardless of LL hybridization. For $B<B_{c1}$, the $\nu=0$ gap
emerges from a series of hybridizations between hole-like (b, $i-$) and
electron-like (B, $i+$) for $i=1,2,\ldots$.

Let us focus on a $\nu=0$ gap at $B_{c1}$ arising from the hybridization
between (b, $1-$) and (B, $1+$). At $B_{c1}^{<}<B_{c1}$, the hybridized LLs at
higher energy (solid blue star) and lower energy (solid red triangle) in Fig.
4(c) are mostly dominated by hole-like (b, $1-$) and electron-like (B, $1+$),
respectively. As we increase $B$ beyond $B_{c1}$, the admixture of each
hybridized LL gradually changes and the situation becomes reversed. Now, the
higher-energy LL (hollow red triangle) evolves into the electron-like (B,
$1+$) while the lower-energy LL (hollow blue star) turns into the hole-like
(b, $1-$). Therefore, as we increase $B$, the characteristic of the higher-
energy LL switches from electron-like to hole-like LL while that of the lower-
energy LL changes from hole-like to electron-like LL. As a result, the edge
states change their behavior significantly across a hybridization gap.

Figure 4(d) illustrates edge state diagrams at $B_{c1}^{<}$ and $B_{c1}^{>}$,
indicated by pink and orange lines in Fig. 4(c). Due to a confining potential,
the energy of an electron-like LL will bend up near the edge while that of a
hole-like LL will bend down. At $B_{c1}^{>}$ [Fig. $4$(d), right], an energy
ordering of the LLs is a conventional one in which the electron-like (B, $1+$)
has higher energy than the hole-like (b, $1-$). In this case, the energies of
both LLs will bend away from each other near the edge. Therefore, these two
LLs do not contribute any edge state to the system at $\nu=0$. However, at
$B_{c1}^{<}$, we have an inverted energy ordering of the LLs in which the
energy of hole-like (b, $1-$) is higher than that of the electron-like (B,
$1+$) in the bulk. Near the edge, their energies will bend toward each other
(down for (b, $1-$) and up for (B, $1+$)) [Fig. 4(d), left]. We therefore
obtain two counter-propagating edge states even though the filling factor is
zero in the bulk. We note that these counter-propagating edge states are
likely not helical since they can interact via mixing terms. As a result, we
expect an energy gap to open at the crossing.

As we lower $B$ further, we encounter another hybridization gap between (b,
$2-$) and (B, $2+$) at $B_{c2}\sim 1.2$ T [Fig. $4$(a-c)]. With the same
argument as the $B_{c1}$ case, the number of edge states will increase from
$2$ to $4$ when $B$ drops below $B_{c2}$ because two hole-like (b, $1-$) and
(b, $2-$) now sit at higher energy than two electron-like (B, $1+$) and (B,
$2+$). In general, as we move through the hybridization gap between (b, $i-$)
and (B, $i+$), the number of edge states changes from $2(i-1)$ to $2i$. We
find that the positions of the conductance peaks at $\nu=0$ in Fig. 4(b) are
in excellent agreement with the theoretical positions of $B_{ci}$ in Fig. 4(c)
which are magnetic fields at which a slope is zero for a hybridized LL
separating $\nu=0$ and $4$. For $B>B_{c1}$, the energy gap at $\nu=0$ turns
into a trivial LL gap and there is no further inversion of electron-like and
hole-like LLs for LL indices $|n|\geq 1$. Therefore, the number of edge states
stays constant and we no longer observe any oscillation [Fig. 4(b)].
Similarly, following the $\nu=12$ line in Fig. 4(a), we observe
magnetoconductance oscillations when $B\lesssim 1.2$ T even though the
$\nu=12$ gap remains finite in the bulk. We find that this $\nu=12$ gap below
$1.2$ T arises from a series of hybridization gaps between (b, $i-$) and (B,
$(i+3)+$) while the gap above $1.2$ T is a trivial LL gap.

We therefore conclude that the oscillations arise from the change in the
number of edge states in a hybridization gap between electron-like and hole-
like LLs. We emphasize that a hybridization gap between LLs of the same type
will not lead to an oscillation because the number of edge states remains
unchanged [30, 32]. Our results show that it is not sufficient to predict
Shubnikov-de Haas oscillations from Landau level spectra only. One needs to
consider if a gap is a trivial Landau level gap or a hybridization gap between
electron-like and hole-like LLs to obtain a complete picture of the
oscillations. We expect our result to be useful for understanding Landau level
spectra of other few layer graphene systems because of the band inversion in
their band structures.

For helical edge modes, two-terminal conductance in the hybridization gap
regime along $\nu=0$ should appear as steps with quantized values of
$4Ne^{2}/h$ where $N$ is the number of edge states and a factor of $4$ comes
from spin and valley degeneracy [9, 10, 4]. For instance, when $\nu=0$ and
$B_{c2}<B<B_{c1}$, we expect conductance of $8e^{2}/h$ from two counter-
propagating edge states from (b, $1-$) and (B, $1+$) but the measured value is
less than $2e^{2}/h$. A few mechanisms could contribute to the observed low
value of conductance. Since the edge states in our system are not helical and
they counterpropagate on the same edge, these two states could mix and tunnel
to each other. As a result, they form 1D localized states and conductance is
no longer quantized at $4Ne^{2}/h$ because the edge states do not have a
perfect transmission [33]. An interaction between the edge states could induce
a small gap at the Fermi energy, reducing conductance further. In addition,
our conductance appears oscillatory which is likely due to the geometry effect
observed in a long sample for two-terminal measurement [23, 24].

In summary, we study magnetotransport properties of 4LG. We observe LL
crossings between the b and B subbands. At finite displacement field, we find
that the LLs in the b subband become valley polarized while those in the B
subband remain valley degenerate, in agreement with the TB calculation with
finite potential difference. At low $n$ and $B$, the band inversion gives rise
to a series of bulk hybridization gaps between electron- and hole-like LLs. As
a result, alternating characteristic of the hybridized LLs between electron
and hole states leads to the change in the number of edge states and manifests
as magnetoconductance oscillations in our measurement. Finally, our proposed
mechanism for magnetoconductance oscillations should be applicable to other
Bernal-stacked multilayer graphenes since they also host similar band
inversion [22, 34].

###### Acknowledgements.

We thank K. Jaruwongrungsee for experimental help and S. Hodak for useful
discussion. This research has been primarily supported by the Research Fund
for DPST graduate with First Placement (Grant no. 002/2015), the NSRF via the
Program Management Unit for Human Resources & Institutional Development,
Research and Innovation (Grant no. B05F640152), and National Research Council
of Thailand (NRCT) and Chulalongkorn University (Grant no. N42A650266). K.W.
and T.T. acknowledge support from JSPS KAKENHI (Grant Numbers 19H05790,
20H00354 and 21H05233).

## References

  * Grushina  _et al._ [2015] A. L. Grushina, D.-K. Ki, M. Koshino, A. A. L. Nicolet, C. Faugeras, E. McCann, M. Potemski, and A. F. Morpurgo, Insulating state in tetralayers reveals an even–odd interaction effect in multilayer graphene, Nat. Commun. 6, 6419 (2015).
  * Wu  _et al._ [2015] Z. Wu, Y. Han, J. Lin, W. Zhu, M. He, S. Xu, X. Chen, H. Lu, W. Ye, T. Han, Y. Wu, G. Long, J. Shen, R. Huang, L. Wang, Y. He, Y. Cai, R. Lortz, D. Su, and N. Wang, Detection of interlayer interaction in few-layer graphene, Phys. Rev. B 92, 075408 (2015).
  * Shi  _et al._ [2018] Y. Shi, S. Che, K. Zhou, S. Ge, Z. Pi, T. Espiritu, T. Taniguchi, K. Watanabe, Y. Barlas, R. Lake, _et al._ , Tunable Lifshitz transitions and multiband transport in tetralayer graphene, Phys. Rev. Lett. 120, 096802 (2018).
  * Che  _et al._ [2020] S. Che, Y. Shi, J. Yang, H. Tian, R. Chen, T. Taniguchi, K. Watanabe, D. Smirnov, C. N. Lau, E. Shimshoni, _et al._ , Helical edge states and quantum phase transitions in tetralayer graphene, Phys. Rev. Lett. 125, 036803 (2020).
  * Mukai  _et al._ [2021] F. Mukai, K. Horii, R. Ebisuoka, K. Watanabe, T. Taniguchi, and R. Yagi, Unconventional satellite resistance peaks in moiré superlattice of h-BN/ AB-stacked tetralayer-graphene heterostructures, Commun. Phys. 4, 109 (2021).
  * Latil and Henrard [2006] S. Latil and L. Henrard, Charge carriers in few-layer graphene films, Phys. Rev. Lett. 97, 036803 (2006).
  * Aoki and Amawashi [2007] M. Aoki and H. Amawashi, Dependence of band structures on stacking and field in layered graphene, Solid State Commun. 142, 123 (2007).
  * Koshino and McCann [2011] M. Koshino and E. McCann, Landau level spectra and the quantum Hall effect of multilayer graphene, Phys. Rev. B 83, 165443 (2011).
  * Young  _et al._ [2014] A. Young, J. Sanchez-Yamagishi, B. Hunt, S. Choi, K. Watanabe, T. Taniguchi, R. Ashoori, and P. Jarillo-Herrero, Tunable symmetry breaking and helical edge transport in a graphene quantum spin Hall state, Nature 505, 528 (2014).
  * Sanchez-Yamagishi  _et al._ [2017] J. D. Sanchez-Yamagishi, J. Y. Luo, A. F. Young, B. M. Hunt, K. Watanabe, T. Taniguchi, R. C. Ashoori, and P. Jarillo-Herrero, Helical edge states and fractional quantum Hall effect in a graphene electron–hole bilayer, Nat. Nanotechnol. 12, 118 (2017).
  * Suzuki  _et al._ [2004] K. Suzuki, K. Takashina, S. Miyashita, and Y. Hirayama, Landau-level hybridization and the quantum Hall effect in $\mathrm{I}\mathrm{n}\mathrm{A}\mathrm{s}/(\mathrm{A}\mathrm{l}\mathrm{S}\mathrm{b})/\mathrm{G}\mathrm{a}\mathrm{S}\mathrm{b}$ electron-hole systems, Phys. Rev. Lett. 93, 016803 (2004).
  * Gusev  _et al._ [2010] G. M. Gusev, E. B. Olshanetsky, Z. D. Kvon, N. N. Mikhailov, S. A. Dvoretsky, and J. C. Portal, Quantum Hall effect near the charge neutrality point in a two-dimensional electron-hole system, Phys. Rev. Lett. 104, 166401 (2010).
  * König  _et al._ [2007] M. König, S. Wiedmann, C. Brüne, A. Roth, H. Buhmann, L. W. Molenkamp, X.-L. Qi, and S.-C. Zhang, Quantum spin Hall insulator state in HgTe quantum wells, Science 318, 766 (2007).
  * Du  _et al._ [2015] L. Du, I. Knez, G. Sullivan, and R.-R. Du, Robust helical edge transport in gated $\mathrm{InAs}/\mathrm{GaSb}$ bilayers, Phys. Rev. Lett. 114, 096802 (2015).
  * Ma  _et al._ [2015] E. Y. Ma, M. R. Calvo, J. Wang, B. Lian, M. Mühlbauer, C. Brüne, Y.-T. Cui, K. Lai, W. Kundhikanjana, Y. Yang, M. Baenninger, M. König, C. Ames, H. Buhmann, P. Leubner, L. W. Molenkamp, S.-C. Zhang, D. Goldhaber-Gordon, M. A. Kelly, and Z.-X. Shen, Unexpected edge conduction in mercury telluride quantum wells under broken time-reversal symmetry, Nat. Commun. 6, 7252 (2015).
  * Karalic  _et al._ [2016] M. Karalic, S. Mueller, C. Mittag, K. Pakrouski, Q. Wu, A. A. Soluyanov, M. Troyer, T. Tschirky, W. Wegscheider, K. Ensslin, and T. Ihn, Experimental signatures of the inverted phase in InAs/GaSb coupled quantum wells, Phys. Rev. B 94, 241402(R) (2016).
  * Yahniuk  _et al._ [2019] I. Yahniuk, S. S. Krishtopenko, G. Grabecki, B. Jouault, C. Consejo, W. Desrat, M. Majewicz, A. M. Kadykov, K. E. Spirin, V. I. Gavrilenko, N. N. Mikhailov, S. A. Dvoretsky, D. B. But, F. Teppe, J. Wróbel, G. Cywiński, S. Kret, T. Dietl, and W. Knap, Magneto-transport in inverted HgTe quantum wells, npj Quantum Mater. 4, 13 (2019).
  * Chen and Ng [2019] Z. Chen and T. K. Ng, Interaction-induced edge states in HgTe/CdTe quantum wells under a magnetic field, Phys. Rev. B 99, 235157 (2019).
  * Koshino [2010] M. Koshino, Interlayer screening effect in graphene multilayers with $ABA$ and $ABC$ stacking, Phys. Rev. B 81, 125304 (2010).
  * Yagi  _et al._ [2018] R. Yagi, T. Hirahara, R. Ebisuoka, T. Nakasuga, S. Tajima, K. Watanabe, and T. Taniguchi, Low-energy band structure and even-odd layer number effect in AB-stacked multilayer graphene, Sci. Rep. 8, 13018 (2018).
  * Wang  _et al._ [2013] L. Wang, I. Meric, P. Huang, Q. Gao, Y. Gao, H. Tran, T. Taniguchi, K. Watanabe, L. Campos, D. Muller, _et al._ , One-dimensional electrical contact to a two-dimensional material, Science 342, 614 (2013).
  * Hirahara  _et al._ [2018] T. Hirahara, R. Ebisuoka, T. Oka, T. Nakasuga, S. Tajima, K. Watanabe, T. Taniguchi, and R. Yagi, Multilayer graphene shows intrinsic resistance peaks in the carrier density dependence, Sci. Rep. 8, 13992 (2018).
  * Abanin and Levitov [2008] D. A. Abanin and L. S. Levitov, Conformal invariance and shape-dependent conductance of graphene samples, Phys. Rev. B 78, 035416 (2008).
  * Williams  _et al._ [2009] J. R. Williams, D. A. Abanin, L. DiCarlo, L. S. Levitov, and C. M. Marcus, Quantum Hall conductance of two-terminal graphene devices, Phys. Rev. B 80, 045408 (2009).
  * Taychatanapat  _et al._ [2011] T. Taychatanapat, K. Watanabe, T. Taniguchi, and P. Jarillo-Herrero, Quantum Hall effect and Landau-level crossing of Dirac fermions in trilayer graphene, Nat. Phys. 7, 621 (2011).
  * Yin  _et al._ [2017] L.-J. Yin, K.-K. Bai, W.-X. Wang, S.-Y. Li, Y. Zhang, and L. He, Landau quantization of Dirac fermions in graphene and its multilayers, Front. Phys. 12, 127208 (2017).
  * [27] See Supplemental Material at [url] for effect of potential difference and additional data from other devices .
  * Serbyn and Abanin [2013] M. Serbyn and D. A. Abanin, New Dirac points and multiple Landau level crossings in biased trilayer graphene, Phys. Rev. B 87, 115422 (2013).
  * Campos  _et al._ [2016] L. C. Campos, T. Taychatanapat, M. Serbyn, K. Surakitbovorn, K. Watanabe, T. Taniguchi, D. A. Abanin, and P. Jarillo-Herrero, Landau level splittings, phase transitions, and nonuniform charge distribution in trilayer graphene, Phys. Rev. Lett. 117, 066601 (2016).
  * Shimazaki  _et al._ [2016] Y. Shimazaki, T. Yoshizawa, I. V. Borzenets, K. Wang, X. Liu, K. Watanabe, T. Taniguchi, P. Kim, M. Yamamoto, and S. Tarucha, Landau level evolution driven by band hybridization in mirror symmetry broken ABA-stacked trilayer graphene, arXiv:1611.02395 [cond-mat.mes-hall] (2016).
  * Zibrov  _et al._ [2018] A. A. Zibrov, P. Rao, C. Kometter, E. M. Spanton, J. I. A. Li, C. R. Dean, T. Taniguchi, K. Watanabe, M. Serbyn, and A. F. Young, Emergent Dirac gullies and gully-symmetry-breaking quantum Hall states in $ABA$ trilayer graphene, Phys. Rev. Lett. 121, 167601 (2018).
  * Zhang  _et al._ [2006] X. C. Zhang, I. Martin, and H. W. Jiang, Landau level anticrossing manifestations in the phase-diagram topology of a two-subband system, Phys. Rev. B 74, 073301 (2006).
  * Lee and Ramakrishnan [1985] P. A. Lee and T. V. Ramakrishnan, Disordered electronic systems, Rev. Mod. Phys. 57, 287 (1985).
  * Horii  _et al._ [2019] K. Horii, T. Nakasuga, T. Hirahara, S. Tajima, R. Ebisuoka, K. Watanabe, T. Taniguchi, and R. Yagi, Magnetotransport study of the mini-Dirac cone in AB-stacked four- to six-layer graphene under perpendicular electric field, Phys. Rev. B 100, 245420 (2019).

Supplemental Material

Effect of Potential Difference on Landau Levels

In our experiment, the potential difference between adjacent layers $U$ varies
with carrier density but, in the calculation of Landau levels, we use a
constant value of the potential difference. Therefore, the calculated Landau
levels may not agree with the experimental data for a whole range of density
values. One of the discrepancies is an absence of $\nu=28$ gap at high field
in our data. To account for it, we consider the effect of $U$ on the Landau
level spectra.

Figures S1(a, b) and S1(c, d) show Landau level spectra for $U=16$ meV and
$U=22$ meV, respectively. As $U$ becomes larger, the LL (b, 1+) shifts up in
energy relatively to the LL (B, 5+). As a result,the energy at which LL (b,
1+) and LL (B, 5+) cross increases. With varying $U$, the LL (b, $1+$) could
disperse along the LL (B, 5+) which leads to the absence of nu = 28 gap over a
considerable range of magnetic field, similar to the data in Fig. 2(b).

Figure S1: (a, b) Landau level spectra and the calculated density of states of
the bulk 4LG as a function of n and B, using the potential difference between
adjacent layers $U$ of $16$ meV. (c, d) Same with (a, b) but with $U=22$ meV.
The position of the crossing between (B, 5+) and (b, 1+) can be seen to shift
up in energy which causes the absence of $\nu=28$ gap at $7$ T.

Additional Data on Unexpected Oscillations at Zero Density

We have measured conductance on 3 devices [see Fig. S2], two of which are two-
terminal devices (Device A and B) and the other is a four-terminal device
(Device C). Figures 2 and 4 from the main text use data from Device A and
Device B, respectively. All three devices display very similar oscillations in
which positions of the peaks occur at almost identical values of magnetic
field.

Figure S2: (a, b) Color map of d$g$/d$n$ from device A at 0 – 7 T and at 0 – 3
T, respectively. (c, d) Color map of d$g$/d$n$ from device B at 0 – 7 T and at
0 – 3 T, respectively. (e, f) Color map of $R_{xx}$ measured in device C with
4-terminal geometry at 0 – 7 T and at 0 – 3 T, respectively. (g, h)
Conductance as a function of magnetic field at $n=0$ from devices A (blue
line), B (red line), and $1/R_{xx}$ as a function of magnetic field at $n=0$
from device C (black line) at 0 – 7 T and at 0 – 3 T, respectively. The
oscillations in all 3 devices show consistent behaviors.

Periodicity in $1/B$ of magnetoconductance oscillations at $n=0$

We observe that the magnetoconductance oscillations are periodic in inverse
magnetic field, suggesting the origin related to Landau levels. We have
included a plot of conductance oscillations with a smooth background
subtracted as a function of inverse magnetic field (Fig. S3).

Figure S3: (a) Conductance as a function of magnetic field at $n=0$ and
$T=2.4$ K. (b) Conductance from (a) with a smooth background subtracted as a
function of inverse magnetic field.