# Non-identical moiré twins in bilayer graphene

E. Arrighi† Université Paris-Saclay, CNRS, Centre de Nanosciences et de
Nanotechnologies (C2N), 91120 Palaiseau, France  V.-H. Nguyen† Institute of
Condensed Matter and Nanosciences, Université catholique de Louvain
(UCLouvain), 1348 Louvain-la-Neuve, Belgium  M. Di Luca  Université Paris-
Saclay, CNRS, Centre de Nanosciences et de Nanotechnologies (C2N), 91120
Palaiseau, France  G. Maffione  Université Paris-Saclay, CNRS, Centre de
Nanosciences et de Nanotechnologies (C2N), 91120 Palaiseau, France  K.
Watanabe  National Institute for Materials Science, 1-1 Namiki, Tsukuba, Japan
T. Taniguchi  National Institute for Materials Science, 1-1 Namiki, Tsukuba,
Japan  D. Mailly  Université Paris-Saclay, CNRS, Centre de Nanosciences et de
Nanotechnologies (C2N), 91120 Palaiseau, France  J.-C. Charlier  Institute of
Condensed Matter and Nanosciences, Université catholique de Louvain
(UCLouvain), 1348 Louvain-la-Neuve, Belgium  R. Ribeiro-Palau  Université
Paris-Saclay, CNRS, Centre de Nanosciences et de Nanotechnologies (C2N), 91120
Palaiseau, France

The superlattice obtained by aligning a monolayer graphene and boron nitride
(BN) inherits from the hexagonal lattice a $60^{\circ}$ periodicity with the
layer alignment. This implies that the mechanical, optical and electronic
properties of the heterostructure must be identical for $0^{\circ}$ and
$60^{\circ}$ of layer alignment. Surprisingly, this periodicity apply only for
the monolayer case. Here, we demonstrate that the moiré superlattice formed by
a bilayer graphene aligned with BN, is present every $60^{\circ}$, but the
symmetry is broken between the $0^{\circ}$ and $60^{\circ}$ alignments,
creating non-identical moiré twins with different electronic properties. In
particular, electron transport measurements display a fully developed valley
Hall effect at $0^{\circ}$ while it is completely absent at $60^{\circ}$. The
origin of the non-identical moiré twins is explained by studying the
relaxation of the atomic structures using state-ot-the-art numerical
techniques, which highlight the central role of the Bernal stacking
configuration with respect to the BN layer. Our results challenge the current
understanding of the valley Hall effect, since a simple Berry curvature
argument do not hold to explain this $120^{\circ}$ periodicity in bilayer
graphene/BN. This unique interplay between mechanical and electronic
properties, demonstrated by the in situ control of the rotational order,
increases the possibilities for band-structure engineering on van der Waals
heterostructures.

When the crystallographic alignment of monolayer graphene and BN is almost
perfect (close to zero degrees between layers), the electronic, mechanical and
optical properties of graphene are strongly modified Hunt2013 ; Wang2015 ;
Ponomarenko2013 ; Ribeiro-Palau2018 . This is caused by the combination of two
effects: i) a long-wavelength geometric interference pattern, called a moiré
pattern, which effectively acts as a periodic superlattice, and ii) a local
enlargement of the lattice constant of graphene to match the one of BN at the
inner part of the moiré pattern, leading to a local commensurate state.
Outside of the commensurate areas the accumulated stress, due to the
stretching of the lattice constant, is released in the form of out-of-plane
corrugations where the stacking order changes rapidly in space Woods2014 .
These corrugations have the same periodicity as the moiré pattern. For
monolayer graphene, both the long-wavelength pattern and the commensurate
state are observed every time one of the layers is rotated by sixty degrees.

The commensurate state creates an imbalance of the interaction that the carbon
atoms have with the BN substrate. This imbalance of carbon-on-boron and
carbon-on-nitrogen stacking breaks the sublattice symmetry Woods2014 . In
monolayer graphene/BN the breaking of inversion symmetry has been proposed as
the origin of the opening of an energy gap at the charge neutrality point
(CNP) Jung2015  and non-trivial quantum geometry characteristics of the
electronic band structure Song2015 . However, very little is known about how
the graphene/BN alignment affects systems with more than one layer, such as
bilayer graphene.

Here, we demonstrate that the commensurate state is not directly transmitted
to consecutive layers of graphene. Instead it is strongly dependent on the
atomic details of the crystallographic alignment and therefore give rise to a
new $120^{\circ}$ periodicity. We present experimental results from local and
non-local electron transport measurements in dynamically rotatable van der
Waals heterostructures Ribeiro-Palau2018  made of Bernal stacked bilayer
graphene and BN, Fig. 1a. Our measurements focus on three angular alignments
between the graphene and BN layers: $0^{\circ}$, where we observe the fully
developed valley Hall effect, $30^{\circ}$ where the electron transport is
perturbed by localized states generated by disorder, and $60^{\circ}$ which
show a standard ohmic response. This surprising difference between $0^{\circ}$
and $60^{\circ}$ cannot be explained by arguments of Berry curvature and
therefore challenge the current understanding of the valley Hall effect in
aligned graphene/BN heterostructures. However, comparing our experimental
results with atomistic numerical simulations, we explain the existence of a
different transport response at $0^{\circ}$ and $60^{\circ}$ by an
inhomogeneous atomic configuration in the commensurate state. This highlights
the critical importance of the interplay between structural and electronic
properties when engineering new van der Waals heterostructures.

Figure 1: Dynamically rotatable heterostructure. a, Schematic representation
(left) of a dynamically rotatable van der Waals heterostructure, the BN
between the graphite and graphene layers has been omitted for clarity. The
central circular shape on the graphene represents the range of action of the
graphite gate and the moiré superlattice. Cross section of the same
heterostructure (right). b, AFM images at $0^{\circ}$, $30^{\circ}$ and
$60^{\circ}$ of alignment. Dashed lines highlight the Hall bar shape of the
graphene layer, $W=1.7~{}\mu$m and $L=2.3~{}\mu$m. Scale bar 3 $\mu$m. c-d,
Four probes resistance measurement as a function of carrier density for
several angular alignment between the bilayer graphene and the BN handle, from
completely misaligned (brown) to fully aligned (light blue) at room
temperature (c) and at 20 K (d) for a selection of angular alignments. The
misalignment in (d) is calculated from the position of the satellite peaks.

Figure 1a shows a cartoon schematic of our device and its cross section. The
dynamically rotatable van der Waals heterostructures are realized as described
in Ribeiro-Palau2018 , with the improvement of having a pre-shaped local
graphite gate. The latter controls the carrier density only in the central
area of our device, and has the same dimensions as the BN structure used to
create the moiré pattern. It is important to mention that the bottom BN and
graphene layers are intentionally misaligned, to more than $10^{\circ}$, to
avoid a double moiré formation Sun2021 ; Wang2019 ; Finney2019 . The carrier
density of the external parts of graphene is tuned by the global Si gate,
acting effectively as a tunable contact resistance Ribeiro-Palau2019 . The
angular alignment of the bilayer graphene/BN heterostructure is controlled in
situ by means of a pre-shaped BN handle deposited on top of graphene. This
handle can be rotated by applying a lateral force with the tip of an atomic
force microscope (AFM), Fig. 1a. The AFM images on Fig. 1b show the three main
positions described in this report: $0^{\circ}$, $30^{\circ}$ and
$60^{\circ}$. An initial alignment, which we refer to in the following as
$0^{\circ}$, is fixed at room temperature inside the AFM using electron
transport measurements (see below), after which the sample is moved to a
cryostat for low temperature experiments. The carrier mobility of the main
sample presented here ranged from 150.000 to 200.000 cm2V-1s-1 for
intermediate densities $\pm 0.65$x$10^{12}$ cm-2 at $T<10$ K. The mean free
path was calculated to be between 1.2 $\mu$m and 2 $\mu$m for the same carrier
density and temperature range (see supplementary information for details). All
measurements were taken using lock-in amplifiers at $f\approx 33.37$ Hz and
applied currents of 10 nA. The non-local voltages were measured using a high
input impedance voltage amplifiers to ensure the measurement had no leaking
current effects (for more details see supplementary information).

In contrast with monolayer graphene Ribeiro-Palau2018 , for bilayer aligned
with BN the presence of the satellite peaks in the electron transport
measurements - a clear signature of the moiré superlattice - becomes evident
only at low temperatures. At room temperature, these satellite peaks are not
visible in a four probe measurements, instead the signature of
crystallographic alignment is given by a broadening of the resistance peak
around the charge neutrality point (CNP), Fig. 1c. The combination of room and
low temperature measurements, Fig. 1c and 1d respectively, allows us to have a
calibration of the angular alignment at room temperature.

Figure 2: Local and non-local transport measurements. a Artistic
representation of the local measurement configurations. b, Local measurements
as a function of the applied gate voltage at $T=$3 K. c, Energy gaps for
different angular alignments. Extracted from the local measurements at the CNP
for the sample described in the main text, a second sample and numerical
simulations (see text). d Artistic representation of the non-local measurement
configurations in the absence (left) and presence of the valley Hall effect
(right). e, Non-local resistance $R_{\mathrm{NL}}=V_{\mathrm{NL}}/I$ as a
function of the gate voltage around the CNP at $T=3$ K.

For both $0^{\circ}$ and $60^{\circ}$ of alignment, the local electron
transport measurement configuration, represented in Fig.2a, shows at low
temperatures the presence of very well pronounced satellite peaks at both
sides of the CNP, Fig. 2b. These satellite peaks are accompanied by a sign
reversal of the Hall resistance $R_{\mathrm{xy}}$ when a low magnetic field is
applied (see supplementary information). In the case of $30^{\circ}$, as
expected, these satellite peaks are not present since graphene and BN are
completely misaligned and the moiré pattern is absent. From magneto-transport
measurements we extract a moiré wavelength of $\lambda$=14.1$\pm$0.4 nm and
$\lambda$=14.5$\pm$0.3 nm for $0^{\circ}$ and $60^{\circ}$, respectively (see
supplementary information). The good coincidence between the positions of the
satellite peaks and very close values of the extracted moiré wavelength for
the measurements at $0^{\circ}$ and $60^{\circ}$ allow us to say that a good
alignment is reached in both cases. Notice that so far, except for a
difference in the magnitude of the local resistance at the satellite peak and
CNP, $0^{\circ}$ and $60^{\circ}$ look indistinguishable.

By following the local measurements at the CNP as a function of the
temperature, at all alignments, we can extract thermally activated energy
gaps, Fig. 2c. At lower temperatures $T<10$ K, the transport at the CNP enters
a hopping regime, characterized by a very slow variation of the conductance
versus $1/T$. The thermally activated behaviour is not observed at the
satellite peaks for $0^{\circ}$ nor $60^{\circ}$, for more detail see
supplementary information.

The measurement configurations is now changed to a non-local one, represented
in Fig. 2d. The non-local electrical signal refers to the appearance of a
voltage across contacts that are well outside the expect path of the current.
This technique is largely used to detect spin/pseudospin signals Wang2015Jan ;
Gorbachev2014 ; Sui2015 ; Shimazaki2015 ; Brune2010 . Non-local signals,
attributed to the existence of valley topological currents, have been measured
previously in aligned monolayer graphene/BN Gorbachev2014 ; Komatsu2018 ;
Li2020 , in bilayer graphene aligned with BN Endo2019  and in the presence of
a strong displacement field Shimazaki2015 ; Sui2015  as well as in other 2D
materials Wu2019 . In this report we focus only on non-local signals at the
CNP, since the non-local signal around the satellite peaks are too weak to be
studied systematically with our current experimental setup. As in previous
reports Gorbachev2014 ; Sui2015 , the non-local resistance, $R_{\mathrm{NL}}$,
appears in a narrow range of gate voltages around the CNP, Fig. 2e, and decays
rapidly with carrier density to values lower than our experimental measurement
noise. Furthermore, we observe a very strong dependence of the non-local
signal with the angular alignment. Indeed, the maximum value of the non-local
resistance at the CNP decreases by a factor four between the measurements at
$0^{\circ}$ and the one at $60^{\circ}$, Fig. 2e.

Figure 3:  Alignment effects in the valley Hall effect. a, Non-local
resistance as a function of the local resistivity $\rho=(W/L)R_{\mathrm{4P}}$,
both at the CNP, for $0^{\circ}$. The color scale indicates different
temperatures. Dashed line indicates the expected value of the non-local
resistance for a valley angle of $\pi/2$, see text. Solid black lines are fit
to the experimental data. b, Non-local resistance as a function of the local
resistivity, both at the CNP, for different angular alignments. Solid lines
are linear fits to the experimental data meant to extract the power law
dependence. Both a and b are measurements between $T=1.4$ K and $T=110$ K.

Plotting the non-local resistance versus local resistivity, $\rho$, at the CNP
for different temperatures, for $0^{\circ}$ of alignment (Fig. 3a), we observe
three regimes: for $T>40$ K an approximately quadratic dependence
$R_{NL}\propto\rho^{2.2}$; for 40 K $\geq T\geq 12$ K a near-cubic relation,
$R_{NL}\propto\rho^{3.2}$, and a saturation regime for $<12$ K. Let us start
the discussion with what happens for $T\leq 40$ K. In analogy with the spin
Hall effect Abanin2009 , this cubic relation is expected in graphene when the
valley Hall and inverse valley Hall effect are in operation. In particular,
the non-local resistance and local resistivity are related by Beconcini2016 :

$R_{\mathrm{NL}}=\frac{W}{2l_{\mathrm{v}}}(\sigma^{\mathrm{v}}_{\mathrm{xy}})^{2}\rho^{3}e^{-L/l_{\mathrm{v}}},$
(1)

when $\sigma^{\mathrm{v}}_{\mathrm{xy}}\ll\sigma$. Here, $l_{\mathrm{v}}$ is
the inter-valley scattering length, $\sigma^{\mathrm{v}}_{\mathrm{xy}}$ is the
valley Hall conductivity, $\rho=1/\sigma$ is the local resistivity and $W$ and
$L$ are the width and length of the sample, respectively. The regime
$\sigma^{\mathrm{v}}_{\mathrm{xy}}\ll\sigma$ is represented in Fig. 2d (left).
The range of temperatures $T\leq 40$ K, for which this cubic relation is
observed, is consistent with the opening of the energy gap
($E_{\mathrm{g}}^{0^{\circ}}\approx 51\pm 11$ K), Fig. 2c.

For $T\leq 12$ K, Fig. 3a, there is a clear change of behavior, characterized
by a saturation of the non-local resistance. This is consistent with the
regime where $\sigma^{\mathrm{v}}_{\mathrm{xy}}\gg\sigma$, and the non-local
response becomes independent of the local resistivity Beconcini2016 :

$R_{\mathrm{NL}}=\frac{W}{2l_{\mathrm{v}}}\frac{1}{\sigma^{\mathrm{v}}_{\mathrm{xy}}}.$
(2)

Since the Fermi energy is placed at the CNP, and the temperature of the system
is much lower than the energy gap, the valley Hall conductivity is maximal
having a value of $\sigma^{\mathrm{v}}_{\mathrm{xy}}=4e^{2}/h$ ($e$ is the
electron’s charge and $h$ is Planck’s constant). In this case all occupied
states in the valence band contribute to the valley Hall effect Yamamoto2015 ;
Beconcini2016 ; Gorbachev2014 ; Sui2015 ; Shimazaki2015 ; Li2020 . Using the
saturation value of the non-local resistance, dashed horizontal line of Fig.
3a, and maximum valley Hall conductivity we obtain an inter-valley scattering
length $l_{\mathrm{v}}=1.5~{}\mu$m, in good agreement with previous reports
Gorbachev2014 ; Shimazaki2015 ; Sui2015 . This behavior is characteristic of a
system where the valley conductivity is larger than the local conductivity,
which implies a fully developed valley Hall effect, where the Hall angle
becomes $\pi/2$, represented in Fig.2d (right).

If we now change the crystallographic alignment of the layer, by in situ
rotation of the BN handle with the AFM tip at room temperature, we can see
that the $60^{\circ}$ case is very different, Fig. 3b. In this case the
relation between the local and non-local resistance is close to linear,
$R_{\mathrm{NL}}\propto\rho^{0.9}$. This behavior is consistent with an ohmic
contribution $R_{\mathrm{NL}}=4\rho e^{-\pi L/W}/\pi$, and can be adjusted by
using only the geometry of our sample, as is expected given its van der Pauw
geometry. The linear behavior can be observed over the full temperature range,
and it does not change when the temperature of the system becomes comparable
with the energy gap, ($E_{\mathrm{g}}^{60^{\circ}}\approx 59\pm 16$ K) Fig.
2c, in contrast to the case for the $0^{\circ}$ of alignment. This striking
difference of the development of the valley Hall effect at $0^{\circ}$ but not
at $60^{\circ}$ reveals that the consequences of the moiré patterns for the
two alignments are non-identical. Changing the alignment further to
$120^{\circ}$ restores the signatures of the valley Hall effect (for more
details see supplementary information).

Figure 4: Atomic relaxation for bilayer graphene/BN. a-b, In-plane
displacement, Dxy, of the carbon atoms for both graphene layers at $0^{\circ}$
and $60^{\circ}$ of alignment. Pink dashed lines in a and b are guides for the
eye to follow the moiré superlattice. c, Atomistic sketch of a BN/bilayer
graphene for $0^{\circ}$ (left) and $60^{\circ}$ (right) alignments. d, AFM
image in PeakForce mode, height sensor, of an aligned bilayer graphene BN.
Scale bar 7 nm. The angular alignment ($0^{\circ}$ or $60^{\circ}$) is
unknown.

Previously the observation of the valley Hall effect has been attributed to
the existence of a Berry curvature, $\Omega(\mathbf{k})$ Gorbachev2014 ;
Komatsu2018 ; Shimazaki2015 ; Sui2015 ; Li2020 , which acts as a pseudo-
magnetic field, and gives rise to an “anomalous” velocity, perpendicular to
the external longitudinal electric field $\mathbf{E}$. The larger the Berry
curvature, the stronger will be this anomalous velocity. The Berry curvature
will spike to a maximum at small values of the energy gap and then reduce
rapidly as the energy gap increases (see supplementary information). This has
recently been demonstrated for bilayer graphene with a strong displacement
field using circularly polarized light Yin2022 . However, our experimental
results show the same energy gap, within error bars, for $0^{\circ}$ and
$60^{\circ}$. This makes the origin of the valley Hall effect at $0^{\circ}$,
and its absence at $60^{\circ}$, more puzzling, requiring further theoretical
investigation.

Interestingly, neither the valley Hall effect, observed at $0^{\circ}$, or the
ohmic behavior, observed for $60^{\circ}$, are reproduced in the fully
misaligned case, $30^{\circ}$, Fig. 3b. In this case, where no signatures of
alignment are observed in the local electron transport or magneto-transport
measurements, we observe a $R_{\mathrm{NL}}\propto\rho^{1.6}$ relation. This
behavior, unrelated to the valley Hall effect, could be explained by the
existence of localized edge states resulting from disorder. It leads to a
superposition of the local ohmic contribution and the conductivity of these
localized, non-topological edge states Aharon-Steinberg2021 . The existence of
these localized states contrasts with the ohmic response observed for the
$60^{\circ}$ of alignment. Our hypothesis is that the periodic potential,
created by the moiré superlattice, may prevent the formation of the localized
states, in the same way as artificial disorder generated by a scanning gate
does, as reported in Aharon-Steinberg2021 . At this point it is important to
point out that we have reproduced these experiments in a “not-optimized
sample” (with a global graphite gate) in which the exponent for $30^{\circ}$
is not the same as the one observed here. This supports our hypothesis that
disorder plays an important role in the misaligned case (see supplementary
information).

To clarify the difference between the two angular alignments we investigate
the in-plane structural relaxation of each layer, Fig. 4a-b. We use atomistic
numerical simulations to obtain the atomic relaxation of each layer and its
consequence in the electronic band structure of both alignments.

Classical molecular dynamics is used to relax the system i.e., minimize forces
and energy. In particular, intralayer forces are computed using the optimized
Tersoff and Brenner potentials Lindsay2010 , while interlayer interactions are
modeled using the Kolmogorov-Crespi potentials Leven2016 ; Kolmogorov2005 .
The lattice mismatch ($\sim$1.8%) between BN and graphene is taken into
account. The atomic structure is optimized until all force components are
smaller than 1 meV/atom. As illustrated in Fig. 4a and b, the in-plane atomic
displacement, $D_{\mathrm{xy}}$, clearly shows that, for the layer that is
closer to the BN (layer 1), there is an almost circular symmetry around the
centre of each moiré superlattice (marked by the pink dashed lines). On the
other hand, the second layer shows a breaking of this symmetry into a $2\pi/3$
rotational symmetry. The in-plane atomic displacement of the second layer is
also at least one order of magnitude smaller than for the first layer.
Additionally, we can see that the in-plane atomic displacement on the second
layer is larger in the case of $0^{\circ}$. The differences in the stretching
of each layer results in a the spatial variation of stacking structure
(initially, perfect AB stacking) as illustrated in the supplementary
information.

The difference of the in-plane relaxations for $0^{\circ}$ and $60^{\circ}$
can be traced back to the Bernal stacked configuration. For this stacking the
distribution of boron and nitrogen atoms for the second layer is not the same
in the $0^{\circ}$ and $60^{\circ}$ cases, Fig. 4c. The direct consequence of
the Bernal stack is that, in a commensurate state, while half of the atoms in
the second layer do not have a direct neighbouring atom, the other half has a
boron neighbouring atom in the $0^{\circ}$ and a nitrogen atom in the
$60^{\circ}$ orientation. Since the interaction of carbon-on-boron and carbon-
on-nitrogen are not the sameJung2015  this creates an inhomogeneous stretch of
the second layer. This inhomogeneous atomic configuration is at the heart of
our non-identical moiré twins.

As we mentioned before, the stress generated by the commensurated state is
released in the form of corrugations, as in the case of monolayer graphene
Woods2014 . These corrugations are transmitted to the second layer and can be
observed in the height sensor of our AFM measurements (PeakForce mode) of a
different bilayer graphene/BN aligned sample, Fig. 4d (for more details see
supplementary information). In this AFM image, a moiré pattern of $\lambda\sim
14$ nm is clearly observed. Note that the asymmetry on the AFM image and large
value on the height of the deformation is an artefact given by the width of
the AFM tip being comparable to the size of the features we want to measure
($\sim 5$ nm). This qualitatively confirms the existence of a commensurate
state and the transmission of the corrugations to the second layer in aligned
bilayer graphene/BN heterostructures, supporting our numerical simulations.

The existence of the commensurate state reveals the challenge of explaining
the existence of the valley Hall effect in our non-identical moiré twins with
the current theoretical understanding. It also highlights that simple
arguments of relation of the energy gap size and Berry curvature, currently
used to explain the valley Hall effect in bilayer graphene in a displacement
field, do not hold to explain the existence of the valley Hall effect in our
case. This difference has its origin in the fact that bilayer graphene in the
presence of a strong displacement field is a much more homogeneous system. In
the case of aligned bilayer graphene/BN there are extra elements to consider
such as the in-plane atomic displacement calling for further theoretical
developments.

To conclude, our experimental results show the existence of non-identical
moirés in bilayer graphene aligned with BN. For $0^{\circ}$ of alignment we
observe a fully developed valley Hall effect that is completely absent at
$60^{\circ}$. We attribute this difference in the atomic configuration in the
commensurated state. These results present major progress in the understanding
of the effects of moiré superlattice and layer strain in the electronic
properties of van der Waals heterostructure. Although the origin of the valley
Hall effect is still not fully explained, we hope that our experimental and
numerical results inspire further theoretical and experimental developments to
address this puzzling problem.

## References

  * (1) D. A. Abanin, A. V. Shytov, L. S. Levitov, and B. I. Halperin.  Nonlocal charge transport mediated by spin diffusion in the spin Hall effect regime.  Phys. Rev. B, 79(3):035304, Jan 2009.
  * (2) A. Aharon-Steinberg, A. Marguerite, D. J. Perello, K. Bagani, T. Holder, Y. Myasoedov, L. S. Levitov, A. K. Geim, and E. Zeldov.  Long-range nontopological edge currents in charge-neutral graphene - Nature.  Nature, 593:528–534, May 2021.
  * (3) M. Beconcini, F. Taddei, and M. Polini.  Nonlocal topological valley transport at large valley Hall angles.  Phys. Rev. B, 94(12):121408, Sep 2016.
  * (4) C. Brüne, A. Roth, E. G. Novik, M. König, H. Buhmann, E. M. Hankiewicz, W. Hanke, J. Sinova, and L. W. Molenkamp.  Evidence for the ballistic intrinsic spin Hall effect in HgTe nanostructures.  Nat. Phys., 6:448–454, June 2010.
  * (5) K. Endo, K. Komatsu, T. Iwasaki, E. Watanabe, D. Tsuya, K. Watanabe, T. Taniguchi, Y. Noguchi, Y. Wakayama, Y. Morita, and S. Moriyama.  Topological valley currents in bilayer graphene/hexagonal boron nitride superlattices.  Appl. Phys. Lett., 114(24):243105, Jun 2019.
  * (6) N. R. Finney, M. Yankowitz, L. Muraleetharan, K. Watanabe, T. Taniguchi, C. R. Dean, and J. Hone.  Tunable crystal symmetry in graphene–boron nitride heterostructures with coexisting moiré superlattices.  Nat. Nanotechnol., 14:1029–1034, Nov. 2019.
  * (7) R. V. Gorbachev, J. C. W. Song, G. L. Yu, A. V. Kretinin, F. Withers, Y. Cao, A. Mishchenko, I. V. Grigorieva, K. S. Novoselov, L. S. Levitov, and A. K. Geim.  Detecting topological currents in graphene superlattices.  Science, Oct 2014.
  * (8) B. Hunt, J. D. Sanchez-Yamagishi, A. F. Young, M. Yankowitz, B. J. LeRoy, K. Watanabe, T. Taniguchi, P. Moon, M. Koshino, P. Jarillo-Herrero, and R. C. Ashoori.  Massive Dirac Fermions and Hofstadter Butterfly in a van der Waals Heterostructure.  Science, 340(6139):1427–1430, Jun 2013.
  * (9) J. Jung, A. M. DaSilva, A. H. MacDonald, and S. Adam.  Origin of band gaps in graphene on hexagonal boron nitride.  Nat. Commun., 6(6308):1–11, Feb 2015.
  * (10) A. N. Kolmogorov and V. H. Crespi.  Registry-dependent interlayer potential for graphitic systems.  Phys. Rev. B, 71(23):235415, Jun 2005.
  * (11) K. Komatsu, Y. Morita, E. Watanabe, D. Tsuya, K. Watanabe, T. Taniguchi, and S. Moriyama.  Observation of the quantum valley Hall state in ballistic graphene superlattices.  Sci. Adv., May 2018.
  * (12) I. Leven, T. Maaravi, I. Azuri, L. Kronik, and O. Hod.  Interlayer Potential for Graphene/h-BN Heterostructures.  J. Chem. Theory Comput., 12(6):2896–2905, Jun 2016.
  * (13) Y. Li, M. Amado, T. Hyart, Grzegorz. P. Mazur, and J. W. A. Robinson.  Topological valley currents via ballistic edge modes in graphene superlattices near the primary Dirac point - Communications Physics.  Commun. Phys., 3(224):1–7, Dec 2020.
  * (14) L. Lindsay and D. A. Broido.  Optimized Tersoff and Brenner empirical potential parameters for lattice dynamics and phonon thermal transport in carbon nanotubes and graphene.  Phys. Rev. B, 81(20):205441, May 2010.
  * (15) P. Moon and M. Koshino.  Electronic properties of graphene/hexagonal-boron-nitride moir\’e superlattice.  Phys. Rev. B, 90(15):155406, Oct 2014.
  * (16) L. A. Ponomarenko, R. V. Gorbachev, G. L. Yu, D. C. Elias, R. Jalil, A. A. Patel, A. Mishchenko, A. S. Mayorov, C. R. Woods, J. R. Wallbank, M. Mucha-Kruczynski, B. A. Piot, M. Potemski, I. V. Grigorieva, K. S. Novoselov, F. Guinea, V. I. Fal’ko, and A. K. Geim.  Cloning of Dirac fermions in graphene superlattices.  Nature, 497:594–597, May 2013.
  * (17) J. Renard, M. Studer, and J. A. Folk.  Origins of Nonlocality Near the Neutrality Point in Graphene.  Phys. Rev. Lett., 112(11):116601, Mar 2014.
  * (18) R. Ribeiro-Palau, S. Chen, Y. Zeng, K. Watanabe, T. Taniguchi, J. Hone, and C. R. Dean.  High-Quality Electrostatically Defined Hall Bars in Monolayer Graphene.  Nano Lett., 19(4):2583–2587, Apr 2019.
  * (19) R. Ribeiro-Palau, C. Zhang, K. Watanabe, T. Taniguchi, J. Hone, and C. R. Dean.  Twistable electronics with dynamically rotatable heterostructures.  Science, Aug 2018.
  * (20) Y. Shimazaki, M. Yamamoto, I. V. Borzenets, K. Watanabe, T. Taniguchi, and S. Tarucha.  Generation and detection of pure valley current by electrically induced Berry curvature in bilayer graphene - Nature Physics.  Nat. Phys., 11:1032–1036, Dec 2015.
  * (21) W. Song Justin C., S. Polnop, and S. Levitov Leonid.  Topological Bloch bands in graphene superlattices.  Proc. Natl. Acad. Sci. U.S.A., 112(35):10879–10883, Sep 2015.
  * (22) M. Sui, G. Chen, L. Ma, W.-Y. Shan, D. Tian, K. Watanabe, T. Taniguchi, X. Jin, W. Yao, D. Xiao, and Y. Zhang.  Gate-tunable topological valley transport in bilayer graphene - Nature Physics.  Nat. Phys., 11:1027–1031, Dec 2015.
  * (23) X. Sun, S. Zhang, Z. Liu, H. Zhu, J. Huang, K. Yuan, Z. Wang, K. Watanabe, T. Taniguchi, X. Li, M. Zhu, J. Mao, T. Yang, J. Kang, J. Liu, Y. Ye, Z. V. Han, and Z. Zhang.  Correlated states in doubly-aligned hBN/graphene/hBN heterostructures.  Nat. Commun., 12(7196):1–8, Dec. 2021.
  * (24) G. Trambly de Laissardière, D. Mayou, and L. Magaud.  Localization of Dirac Electrons in Rotated Graphene Bilayers.  Nano Lett., 10(3):804–808, Mar 2010.
  * (25) L. Wang, Y. Gao, B. Wen, Z. Han, T. Taniguchi, K. Watanabe, M. Koshino, J. Hone, and C. R. Dean.  Evidence for a fractional fractal quantum Hall effect in graphene superlattices.  Science, 350(6265):1231–1234, Dec 2015.
  * (26) L. Wang, S. Zihlmann, M.-H. Liu, P. Makk, K. Watanabe, T. Taniguchi, A. Baumgartner, and C. Schönenberger.  New Generation of Moiré Superlattices in Doubly Aligned hBN/Graphene/hBN Heterostructures.  Nano Lett., 19(4):2371–2376, Apr. 2019.
  * (27) Z. Wang, C. Tang, R. Sachs, Y. Barlas, and J. Shi.  Proximity-Induced Ferromagnetism in Graphene Revealed by the Anomalous Hall Effect.  Phys. Rev. Lett., 114(1):016603, Jan. 2015.
  * (28) C. R. Woods, L. Britnell, A. Eckmann, R. S. Ma, J. C. Lu, H. M. Guo, X. Lin, G. L. Yu, Y. Cao, R. V. Gorbachev, A. V. Kretinin, J. Park, L. A. Ponomarenko, M. I. Katsnelson, Yu. N. Gornostyrev, K. Watanabe, T. Taniguchi, C. Casiraghi, H.-J. Gao, A. K. Geim, and K. S. Novoselov.  Commensurate–incommensurate transition in graphene on hexagonal boron nitride.  Nat. Phys., 10:451–456, Jun 2014.
  * (29) Z. Wu, B. T. Zhou, X. Cai, P. Cheung, G.-B. Liu, M. Huang, J. Lin, T. Han, L. An, Y. Wang, S. Xu, G. Long, C. Cheng, K. T. Law, F. Zhang, and N. Wang.  Intrinsic valley Hall transport in atomically thin MoS2.  Nat. Commun., 10(611):1–8, Feb 2019.
  * (30) M. Yamamoto, Y. Shimazaki, I. V. Borzenets, and S. Tarucha.  Valley Hall Effect in Two-Dimensional Hexagonal Lattices.  J. Phys. Soc. Jpn., 84(12):121006, Nov 2015.
  * (31) J. Yin, C. Tan, D. Barcons-Ruiz, I. Torre, K. Watanabe, T. Taniguchi, J. C. W. Song, J. Hone, and F. H. L. Koppens.  Tunable and giant valley-selective Hall effect in gapped bilayer graphene.  Science, 375(6587):1398–1402, Mar 2022.

## Acknowledgements

The authors acknowledge discussions with Ulf Gennser, Marco Polini, Herve
Aubin and Justin Song. This research was supported the Flag-Era JTC project
“TATTOOS” (N∘ R.8010.19). R.R.-P. also acknowledge the ERC starting grant
TWISTRONICS. This work was partly supported by the French RENATECH network and
the DIM-SIRTEC. V.-H.N. and J.-C.C. acknowledge financial support from the
Fédération Wallonie-Bruxelles through the ARC Grant (N∘ 21/26-116), from the
European Union’s Horizon 2020 Research Project and Innovation Program -
Graphene Flagship Core3 (N∘ 881603), from the Flag-Era JTC project “TATTOOS”
(N∘ R.8010.19), from the EOS project “CONNECT” (N∘ 40007563) and from the
Belgium F.R.S.-FNRS through the research project (N∘ T.029.22F). Computational
resources have been provided by the CISM supercomputing facilities of
UCLouvain and the CÉCI consortium funded by F.R.S.-FNRS of Belgium (N∘
2.5020.11). V.-H.N. thanks Dr. Xuan-Hoang TRINH for his helps in
implementation of numerical codes to compute the lattice relaxation.

## Author Contributions

R.R.-P. and D.M. designed the experiment. E.A. fabricated the devices for
electron transport measurements. G.M. and M.DL. fabricated the samples for
structural characterization and performed the AFM measurements. E.A. and R.R-P
performed the electron transport experiments and analyzed the data. T.T. and
K.W. grew the crystals of hexagonal boron nitride. V.-H.N. and J.-C.C.
performed the numerical simulations and participated to the data analysis. All
authors participated to writing the paper.

## Supplementary Information

### Mean free path and mobility at different angular alignments

We calculate the mean free path versus temperature from the measurements of
resistance as a function of gate voltage at different temperatures and using:

$l_{\mathrm{mfp}}=\frac{\sigma h}{2e^{2}\sqrt{\pi n}}$ (S1)

From Fig. S1 we can see that the sample is in a ballistic regime for $T<10$ K,
shadow area and we also remark that the temperature dependence of the mean
free path for the misaligned case turns our to be different for the misaligned
case ($30^{\circ}$). The understanding of this regime is out of the scope of
our manuscript and it will be addressed in a different scientific publication.

Figure S1: Mean free path calculated from the resistance as a function of
$V_{\mathrm{g}}$ at different temperatures for a density of $+0.64\times
10^{12}$ cm-2 filled symbols and $-0.64\times 10^{12}$ cm-2 (empty symbols)
for different angular alignments. Dashed horizontal line represents the width
of our sample. Shadow area signals the values of temperature where the mean
free path is not varying with temperature.

### Non-local measurements configurations

Measuring non-local voltages has proven to be a challenging task [refs]. To
avoid any common grounds and spurious signals we have followed the measurement
scheme developed in [ref]. This one consist of using an operational amplifier
to keep the voltage of our sample balanced and high impedance amplifiers
(CELIANS ADD REF) to avoid any current leak, see Figure S2.

Figure S2: Non-local measurements. Non-local configuration measurement using
an operational amplifier and high impedance voltage amplifiers.

Our non-local measurements are performed at low current (10 nA)Joule heating
at the injector would cause heat to flow into and past the detector region.
The resulting temperature gradient, along the detector region would give rise
to a non-local voltage across the detector contacts via the Nernst effect,
quantified by the transverse thermopower coefficient. This temperature
gradient is proportional to heating power, quadratic in current, and therefore
contributes to the non-local voltage only at the second harmonic of the
excitation frequency. As a result, Joule heating would not affect the first
harmonic data presented in this report Renard2014

### Hall resistance measurement

Satellite peaks and CNP are accompanied by inversion of the sign in the Hall
resistance measurement. For the $30^{\circ}$ of alignment, full misalignment,
there are not satellite peaks and therefore the Hall resistance stays close to
zero.

Figure S3: Hall resistance measurement. Measurements at low temperature and
0.2 T of the Hall resistance for the three crystallographic alignments
$0^{\circ}$, $30^{\circ}$ and $60^{\circ}$.

### Determining $\lambda$ by magneto transport measurements

By using magneto transport measurements we can determine the magnetic field at
which one flux quantum threads the supperlattice unit cell. To obtain this we
plot the longitudinal resistance as a function of the gate voltage (or carrier
density) and 1$/B$, Fig. S4. In this plot we can see the appearance of equally
spaces horizontal lines that appear when the Landau fans coming from the
charge neutrality point and the satellite peaks intercept each other. It is
important to remark that this is a purely geometrical effect and gives
therefore a very accurate estimation of the moiré superlattice size Hunt2013 .

Figure S4: Magneto transport measurements. Longitudinal resistance, in log
scale, as a function of the applied gate voltage and the inverse magnetic
field for a $0^{\circ}$ and b $60^{\circ}$ of alignment. Horizontal dashed
lines are guides for the eyes.

We obtain a space between these horizontal lines corresponding to
$B_{0}=24.16$ T and $B_{0}=22.69$ T, this corresponds to
$\lambda$=14.1$\pm$0.4 nm and $\lambda$=14.5$\pm$0.3 nm for $0^{\circ}$ and
$60^{\circ}$, respectively.

### Temperature dependence and thermally activated regimes for the local and
non-local resistance

In our local measurements we observe different transport regimes at the CNP:
first a thermally activated regimen, fit with an Arrhenious law and a hopping
regime for temperatures lower than $<10$ K. In Fig. S5, we plot an example for
the activation regime of the CNP for the measurements at $0^{\circ}$ in local
and non-local configuration. The values for the local energy gaps at the CNP
for different alignments are summarized in Table 1.

Figure S5: Temperature dependence local and non-local measurements at the CNP
for $0^{\circ}$. and b, Aarhenious plot for the energy gap at the CNP for the
local and non-local measurements, respectively. c Energy gaps for the sample
described in the main text, the non-optimized sample (sampel II) and numerical
simulations.

Differently at the satellite peaks we do not observe any activation regime,
Fig. S6. This might be due to the presence of a very small energy gap that is
washed out by disorder, a short discussion of this can be found in the main
text.

Figure S6: Temperature dependence local measurements at the satellite peaks
for $0^{\circ}$. a-b measurements at the satellite peak for negative (a) and
positive (b) values of the gate voltage applied to the graphite gate, while
the Si gate is kept at high doping. Energy gaps
---
Angle |  CNP-Local |  CNP-Non-local |  Ratio
|  (K) |  (K) |
0∘ |  51 $\pm$ 11 |  144 $\pm$ 18 |  2.82
30∘ |  87 $\pm$ 17 |  161 $\pm$ 20 |  1.85
60∘ |  59 $\pm$ 16 |  61 $\pm$ 10 |  1.04
120 ∘ |  65 $\pm$ 10 |  144 $\pm$ 20 |  2.21
Table 1: Energy gaps. Summary of the energy gaps for the local and non-local
configurations at different crystallographic alignments.

If we compare the energy gap obtained from the local and no-local
measurements, for the $0^{\circ}$ and $120^{\circ}$ measurements, we can
easily see that, as expected, there is about factor of three of difference
between these energy gaps Yamamoto2015 ; Endo2019 . It is also important to
notice that the extracted values for the energy gaps, local and non-local, are
at least a factor of two larger than previously reported Endo2019 , this
reaffirms the high quality of our samples.

### Non-local measurements at 120 degrees

We have performed the same local versus non-local measurements for the
configuration at $120^{\circ}$. However, at the moment when we performed these
measurements we have lost one contact (#6 in Fig.S7a). This clearly impacted
the saturation regime of our measurements, since it changes the geometrical
configuration of our sample. However, the cubic relation is still visible,
Fig. S7b.

Figure S7: Non-local versus local at $120^{\circ}$. a, measurements
configuration after one contact got broken. b, non-local resistance versus
local resistivity for $0^{\circ}$, $30^{\circ}$, $60^{\circ}$ and
$120^{\circ}$.

### Angular dependence in a non-optimized sample

Our first generation of samples had a global graphite gate instead of a local
one. This generated a very large contact resistance when passing through the
CNP, since there is a part of our graphene flake that is exposed and therefore
not as clean as the part covered by the BN handles. This sample architecture
was improved by the ones presented in the main text.

In this non-optimized samples we can still observe a close to cubic relation
of the non-local and local resistance, $R_{\mathrm{NL}}\propto\rho^{2.7}$.
This relation is strongly modified for $30^{\circ}$ and $60^{\circ}$ of
alignment, Fig. S8.

Figure S8: Non-local versus local resistance in a non-optimized sample.

It is important to point out that these measurements were performed in a dry
$4$ K cryostat with a resistive coil that goes up to 300 mT. As we should
expect the resistive coil have no remnant field and therefore the non-local
signal cannot be attributed to remnant field in the magnet.

Notice that the exponent for this sample in the case of zero degrees is not
exactly the same as the one described in the main text. This difference
highlights the effects of disorder in the samples. It is important to remark
that this small difference do not change the fact that there is a strong
dependence on the layer alignment between $0^{\circ}$ and $60^{\circ}$ .

### Sample for structural characterization

We have prepared a sample to measure with the AFM the out of plane atomic
displacement, solitons formed as a consequence of the atomic in plane
displacement inside the moire cell. This sample consist of BN and bilayer
graphene on surface. The two crystals have been prealigned at the moment of
the preparation by aligning their crystallographic edges. A picture of this
sample and its cross section are shown in Fig. S9.

Figure S9: Sample for structural characterization. Optical micrograph of the
sample (left). Dashed line represents the graphene area. Scale bar 20 $\mu$m.
Cross section of the sample (right).

The bilayer character of the sample has been tested by Raman spectroscopy. The
AFM measurement in this sample, shown in the main text

### Electronic models

To compute the electronic structure of the considered graphene/hBN systems, we
employed the $p_{z}$ tight-binding Hamiltonian, similar to those presented in
TramblydeLaissardiere2010 ; Moon2014 . In particular, the Hamiltonian is
written as

$H_{tb}=\sum_{n}V_{n}a^{\dagger}_{n}a_{n}+\sum_{n,m}t_{nm}a^{\dagger}_{n}a_{m}$

where the on-site energies $V_{n}$ = 0, 3.34 eV, and -1.4 eV for carbon,
boron, and nitride atoms, respectively. The hopping energies $t_{nm}$ are
determined using the standard Slater-Koster formula

$\displaystyle t_{nm}(r_{nm})$ $\displaystyle=$ $\displaystyle
V_{pp\pi}\sin^{2}\phi_{nm}+V_{pp\sigma}\cos^{2}\phi_{nm},$ $\displaystyle
V_{pp\pi}$ $\displaystyle=$ $\displaystyle
V^{0}_{pp\pi}\exp\left((a_{0}-r_{nm})/r_{0}\right),$ $\displaystyle
V_{pp\sigma}$ $\displaystyle=$ $\displaystyle
V^{0}_{pp\sigma}\exp\left((d_{0}-r_{nm})/r_{0}\right)$

where the direction cosine of ${\vec{r}}_{nm}$ along Oz axis is
$\cos\phi_{nm}=z_{nm}/r_{nm}$, $r_{0}=0.184a$, $a_{0}=a/\sqrt{3}$, and
$d_{0}=3.35$Å while $a\simeq 2.49$Å.

### Atomic structure relaxation and electronic band structures

More details related to atomic structure relaxation are presented in Fig. S10,
including additionally the variation of interlayer distance between
graphene/hBN layers and the modification of stacking structure, that is due to
atomic reconstruction effects. In Fig. S11, the electronic band structures
computed for both unrelaxed and relaxed lattices are displayed. While it is
almost negligible in the unrelaxed case, the presented results clearly
demonstrate that the effects of atomic structure relaxation essentially govern
the significant difference in the electronic structures of $0^{\circ}$\- and
$60^{\circ}$\- alignments.

Figure S10: Atomic structure relaxation. Bilayer graphene/hBN are aligned at
$0^{\circ}$ (left) and $60^{\circ}$ (right). The variation of interlayer
distance $d_{z}$ between the second graphene layer (see Fig. 4) and hBN one is
presented in the two top images. The four images below illustrate the
variations of the in-plane displacements $D_{xy}$ in both graphene layers. At
last, the two bottom images display the modification of stacking configuration
between two graphene layers, by illustrating the variation of stacking vector
t. t is determined as the translation vector applied locally to a C-ring of
one graphene layer to recover the AA stacking configuration at the considered
position, i.e., $|$t$|$ = $a_{0}$ and 0 for the local AB and AA stacking
configurations, respectively. Figure S11: Electronic band structures. Bilayer
graphene/hBN are aligned at $0^{\circ}$ (top) and $60^{\circ}$ (bottom). Ideal
pristine atomic structure of the graphene layers are preserved on the left,
while atomic structure relaxations are considered on the right.

### Berry curvature:

if we consider the existence of a Berry curvature, $\Omega(\mathbf{k})$
Gorbachev2014 ; Komatsu2018 ; Shimazaki2015 ; Sui2015 ; Li2020 , which acts as
a pseudo-magnetic field, this will give rise to an “anomalous” velocity,
perpendicular to the external longitudinal electric field $\mathbf{E}$. The
electron velocity, $\mathbf{v}$, becomes then:

$\mathbf{v}(\mathbf{k})=\frac{\partial\epsilon(\mathbf{k})}{\hbar\partial\mathbf{k}}-\frac{q}{\hbar}\mathbf{E}\times\Omega(\mathbf{k}),$
(S2)

where $k$, is the wavevector, $\epsilon(\mathbf{k})$ is the band energy, $q$
is the carrier charge and $\hbar$ is the reduced Planck constant. From
equation S2 here we can see that the larger the Berry curvature, stronger will
be the anomalous velocity. The relation between the Berry curvature and the
energy gap is given by Yin2022 :

$\Omega(\mathbf{p})=\frac{2\hbar^{2}\nu^{4}\gamma\Delta|\mathbf{p}|^{2}}{(\nu^{4}\mathbf{p}^{4}+\gamma^{2}\Delta^{2})^{3/2}},$
(S3)

where $\Delta$ is the value of the energy gap, $\gamma$ is the interlayer
coupling, $\nu$ is the velocity in monolayer and $\mathbf{p}$ is momentum.
This means that the Berry curvature will spike to a maximum at very small
values of the energy gap and then reduce rapidly as the energy gap increases