# Supplementary Information for
Imaging tunable quantum Hall broken-symmetry orders
in graphene

Alexis Coissard1,∗, David Wander1,∗, Hadrien Vignaud1, Adolfo G. Grushin1,
Cécile Repellin2, Kenji Watanabe3, Takashi Taniguchi4, Frédéric Gay1,
Clemens B. Winkelmann1, Hervé Courtois1, Hermann Sellier1 & Benjamin
Sacépé${}^{1,\,\textrm{{\char 0\relax}}}$

  * 1

Univ. Grenoble Alpes, CNRS, Grenoble INP, Institut Néel, 38000 Grenoble,
France

  * 2

Univ. Grenoble-Alpes, CNRS, LPMMC, 38000 Grenoble, France

  * 3

Research Center for Functional Materials, National Institute for Materials
Science, 1-1 Namiki, Tsukuba 305-0044, Japan

  * 4

International Center for Materials Nanoarchitectonics, National Institute for
Materials Science, 1-1 Namiki, Tsukuba 305-0044, Japan

  * *

These authors contributed equally to this work.

  * ${}^{\textrm{{\char 0\relax}}}$

email<EMAIL_ADDRESS>
##  I Sample fabrication

Graphene/hBN heterostructures were assembled from exfoliated flakes with the
van der Waals pick-up technique using a polypropylene carbonate (PPC) polymer
[1]. Stacks were deposited using the methods described in Ref. [2] (for sample
STO07) or in Ref. [3] (for samples AC04, AC23 and AC24), on either highly
doped Si wafers with a $285$ nm thick SiO2 layer, or on 500$\>\mu$m thick
SrTiO3 (100) substrates cleaned with hydrofluoric acid buffer solution before
deposition of the graphene/hBN heterostructures (a Ti/Au bilayer was deposited
later on the other side of the SrTiO3 substrate to enable back-gate effect).
The geometrical parameters of the samples are listed in Table S1. Electron-
beam lithography using a PMMA resist was used to pattern a guiding markerfield
on the whole $5\times 5\>$mm2 substrate to drive the STM tip toward the
device. Cr/Au electrodes contacting the graphene flake were also patterned by
electron-beam lithography and metalized by e-gun evaporation. Samples were
thermally annealed at $350\,^{\circ}$C in vacuum under an halogen lamp to
remove resist residues and clean graphene, before being mounted into the STM
where they were heated in situ during the cooling to $4.2\;$K.

##  II Measurements

Experiments were performed with a home-made hybrid scanning tunneling
microscope (STM) and atomic force microscope (AFM) operating at a temperature
of $4.2\;$K in magnetic fields up to $14\;$T. The AFM mode is used to guide
the tip toward the graphene device. The sensor consists of a hand-cut PtIr tip
glued on the free prong of a tuning fork, the other prong being glued on a
Macor plate. Once mounted inside the STM, the tip is roughly aligned over the
sample at room temperature and then guided toward the graphene in AFM mode at
low temperature using the guiding markerfield. Scanning tunneling spectroscopy
was performed using a lock-in amplifier technique with a modulation frequency
of $263\;$Hz and rms modulation voltage between $1-5\;$mV depending on the
spectral range of interest. Imaging of the zLL lattice-scale orders was
carried out in STM constant-height mode. Starting from a tunneling contact at
($V_{\text{b}}=300\>$mV, $I_{\text{t}}=1\>$nA) with the Z-regulation on, we
switch off the regulation and lower the bias voltage to either energies
corresponding to the LL${}_{0^{\pm}}$ peaks, which drastically decreases the
tunneling current. We then manually approach the tip toward graphene until the
recovery of a tunneling current of a few nA. STM images of the tunneling
current measured while scanning at constant tip-sample distance subsequently
yield atomically resolved images of the honeycomb lattice or lattice-scale
orders.

##  III Samples summary

We studied two different types of samples (see Fig. 1d of the main text) :
first graphene on hBN/SiO2 substrate, the usual substrate used in transport
measurements, and secondly graphene on hBN/SrTiO3 substrate with a thin hBN
layer of thickness $d_{\text{BN}}$. Optical images of the samples are shown in
Fig. S1 and their geometrical parameters are listed in Table S1. Both types of
samples are equipped with a back-gate electrode.

Table S1: Geometrical parameters of the four measured samples. Fig. S1:
Graphene samples. Optical pictures of the studied samples listed in Table S1.
The dashed blue lines outline the hBN flakes, while the white dashed lines
outline the graphene flakes. For every image, the scale bar is $10\>\mu$m.

##  IV Dielectric constant of SrTiO3 in magnetic field.

Fig. S2: Estimation of the dielectric constant of SrTiO3 and rescaling of the
gate map. a, Line cut at $V_{\text{b}}=0\>$V, averaged on a range of $\pm
20\>$mV, of the d$I_{\text{t}}$/d$V_{\text{b}}$ gate map in Fig. 2a. b,
Estimation from the filling factors obtained in a of the charge carrier
density $n$ (blue dots), its polynomial fit (blue curve), and computed values
of $\epsilon_{\text{r}}\simeq\epsilon_{\text{STO}}$ (red curve), as a function
of gate voltage. The fit yields $V^{\text{CNP}}\simeq 13.5\>$V. c, Rescaling
of the gate map of Fig. 2a as a function of $\nu$.

We show here how to estimate the SrTiO3 dielectric constant,
$\epsilon_{\text{STO}}$, from tunneling conductance gate maps.
$\epsilon_{\text{STO}}$ is related to the global dielectric constant of the
back gate, $\epsilon_{\text{r}}$, which can be obtained by modeling the back-
gate capacitance $C_{\text{g}}$ as the sum of the series capacitances of
SrTiO3 and hBN assuming plane capacitors :
$\frac{1}{C_{\text{g}}}=\frac{1}{C_{\text{STO}}}+\frac{1}{C_{\text{BN}}}\Rightarrow\frac{d_{\text{STO}}+d_{\text{BN}}}{\epsilon_{\text{r}}}=\frac{d_{\text{STO}}}{\epsilon_{\text{STO}}}+\frac{d_{\text{BN}}}{\epsilon_{\text{BN}}}.$
Since $d_{\text{BN}}\sim 10\>\text{nm}\ll d_{\text{STO}}=500\>\mu$m, we write
the gate dielectric constant as :

$\epsilon_{\text{r}}=\epsilon_{\text{STO}}\left(1+\frac{d_{\text{BN}}}{d_{\text{STO}}}\frac{\epsilon_{\text{STO}}}{\epsilon_{\text{BN}}}\right)^{-1}$
(1)

Numerically, $\epsilon_{\text{BN}}\simeq 3.6$ and $\epsilon_{\text{STO}}\sim
10^{4}$ at low temperature [4], so that
$d_{\text{BN}}\epsilon_{\text{STO}}/d_{\text{STO}}\epsilon_{\text{BN}}\sim
0.1$. We can thus assume that
$\epsilon_{\text{r}}\simeq\epsilon_{\text{STO}}$.

In order to estimate $\epsilon_{\text{r}}$ as a function of the back-gate
voltage, $V_{\text{g}}$, we consider the tunneling conductance gate map of
Fig. 2a (sample STO07) from which we can extract some values of the back-gate
voltage at specific filling factors $\nu$. Note that the electron-hole
asymmetry visible in this gate map stems from the non-linear behavior of
$\epsilon_{\text{STO}}$ with gate voltage [5, 6, 7]. We plot in Fig. S2a the
line cut of the gate map at zero bias, averaged on a range of $\pm 20\>$meV
around this value. We clearly observe the different non-zero conductance
plateaus forming when $E_{\text{F}}$ is pinned inside one LL, with gaps in-
between. As the gate voltages at the middle of those gaps correspond to
completely filled and empty LLs, we identify the positions in gate voltage of
the integer filling factors $\nu=-10,-6,-2,2$. Those values are converted into
charge carrier density values $n$ in Fig. S2b using $\nu=n\phi_{0}/B$ with
$\phi_{0}=h/e$ the flux quantum and $n$ the charge carrier density. A
polynomial of degree 5 fits and interpolates the evolution of $n$ with
$V_{\text{g}}$. From this fit, charge neutrality at $n=0$ is achieved at
$V^{\text{CNP}}=13.5\>$V. We then straightforwardly obtain the
$V_{\text{g}}$-dependence of $\epsilon_{\text{r}}$ via :

$\epsilon_{\text{r}}=\frac{d_{\text{STO}}}{\epsilon_{0}}\frac{en}{V_{\text{g}}-V^{\text{CNP}}}$
(2)

The red curve in Fig. S2b shows the resulting $\epsilon_{\text{r}}$, which
decreases with increasing gate voltage and ranges between $12\>500$ and
$3\>500$. A similar $\epsilon_{\text{r}}(V_{\text{g}})$ profile but with
slightly weaker values is obtained for sample AC23
($3\>000<\epsilon_{\text{r}}<11\>500$). Finally, using the fit of the filling
factor $\nu$, we can rescale the gate map as a function of $\nu$ as shown in
Fig. S2c. In particular, note that in the rescaled map the interaction-induced
gap is maximal at charge neutrality $\nu=0$, as expected considering that the
exchange interaction is maximal at half-filling of the zeroth Landau level.

##  V Landau levels of Dirac fermions

We show in Fig. S3a a tunneling conductance d$I_{\text{t}}$/d$V_{\text{b}}$
spectrum taken on sample STO07 (SrTiO3) at $B=14\>$T (same than Fig. 2b of the
main text), where graphene was brought to charge neutrality with a gate
voltage $V_{\text{g}}=13\>$V using the d$I_{\text{t}}$/d$V_{\text{b}}$ gate
maps from Fig. 2a. Additionally, Fig. S3b displays a spatially-averaged
d$I_{\text{t}}$/d$V_{\text{b}}$ spectrum on a $100\times 100\>$nm2 area around
the same position as Fig. S3a, where we clearly see well-resolved Landau
levels in the local density of states up to $N=\pm 6$ as well as both peaks of
the LL0 broken-symmetry state, with a gap of $\Delta E^{0}\simeq 36\>$meV. We
fit in Fig. S3c the positions $E_{N}$ of the Landau levels as a function of
$\left(|N|B\right)^{1/2}$ ($E_{0}\approx 0.7\>$meV is taken at the middle of
the LL0 peaks). We obtain an excellent agreement with the theoretical
dispersion relation for graphene :

$E_{N}=E_{\text{D}}+\text{sign}(N)v_{\text{F}}\sqrt{2\hbar e|N|B}$ (3)

confirming the massless behavior of charge carriers in graphene. The fit
yields a Fermi velocity of $v_{\text{F}}=(1.403\pm 0.005)\times
10^{6}\>$m.s-1. This value is much greater than the expected theoretical one
of $1.0\times 10^{6}\>\text{m.s}^{-1}$, which we attribute to the enhancement
of electron-electron interactions at charge neutrality [8, 9, 10].

Fig. S3: Relativistic Landau levels in graphene in sample STO07 on SrTiO3. a,
Individual tunneling differential conductance d$I_{\text{t}}$/d$V_{\text{b}}$
spectrum at $B=14\>$T and $V_{\text{g}}=13\>$V showing well-defined Landau
level peaks. The zeroth Landau level is split into two peaks
LL${}_{0}{{}^{\pm}}$. b, Spatially-averaged d$I_{\text{t}}$/d$V_{\text{b}}$
spectrum under the same conditions on a $100\times 100\>$nm2 area. c, Fit
using Equation (3) of the positions of Landau levels $E_{N}$ as a function of
$(|N|B)^{1/2}$ with $N$ the Landau level index. An excellent agreement is
obtained and the fit yields a Fermi velocity $v_{\text{F}}=(1.403\pm
0.005)\times 10^{6}\>$m.s-1.

Interestingly, we note that this energy spacing of Landau levels is symmetric
with respect to electron and hole levels. This indicates that a possible tip-
induced local doping due to the work-function difference between the tip and
graphene does not yield a significant spatial confinement on the Landau
levels, which could have modified the Landau level spectrum [11].

##  VI Tip-induced gating

Probing graphene density of states in tunneling spectroscopy induces an
unintentional local doping of graphene due to the field effect of the tip. The
total doping of graphene from field effects of both the back gate and the tip
is therefore :

$n=\frac{1}{e}\left(C_{\text{g}}V_{\text{g}}-C_{\text{tip}}V_{\text{b}}\right)$
(4)

with the capacitances $C_{\text{g}}$ between graphene and the back gate and
$C_{\text{tip}}$ between the tip and graphene, and the voltages $V_{\text{g}}$
applied to the back gate and $V_{\text{b}}$ applied to the sample (note that
$V_{\text{b}}\ll V_{\text{g}}$) while the tip is grounded, as shown in Fig.
S4a.

Fig. S4: Tip-induced gating. a, Gating of graphene from both the back gate and
the tip. b, Tunneling conductance gate map, same as Fig. 2a from the main
text. The isodensity lines, represented in red, are tilted and their slope
gives an estimation of the capacitances ratio $C_{\text{g}}/C_{\text{tip}}$.

Constant-carrier-density lines in tunneling conductance gate maps thus do not
appear as vertical lines at constant back-gate voltages, but rather as tilted
lines, see the red lines in Fig. S4b which show that the Landau levels
staircase transitions do not occur at the same back-gate voltage. The slopes
of the isodensity lines are given by :

$\frac{\text{d}n}{\text{d}V_{\text{g}}}=0\Rightarrow\frac{\text{d}V_{\text{b}}}{\text{d}V_{\text{g}}}=\frac{C_{\text{g}}}{C_{\text{tip}}}$
(5)

Since the back-gate capacitance (per unit area) $C_{\text{g}}$ can be
estimated as $C_{\text{g}}=\epsilon_{0}\epsilon_{\text{r}}/d$, with $d$ the
thickness of the back-gate insulator and $\epsilon_{\text{r}}$ its relative
dielectric constant, and the slope $C_{\text{g}}/C_{\text{tip}}$ of the
isodensity lines can be obtained from the gate maps, following Ref. [10], it
is possible to deduce the capacitance $C_{\text{tip}}$. The variation of the
carrier density due to tip-induced gating at constant back-gate voltage is
then :

$\delta n_{\text{tip}}=\frac{1}{e}C_{\text{tip}}V_{\text{b}}$ (6)

and the resulting tip-induced variation of the filling factor is :

$\delta\nu_{\text{tip}}=\frac{\delta n_{\text{tip}}\phi_{0}}{B}$ (7)

with $\phi_{0}$ the flux quantum and $B$ the magnetic field.

###  VI.1 SiO2 sample

The back-gate insulator is here comprised of the SiO2 layer and the hBN flake.
Its thickness is $d=d_{\text{SiO}_{2}}+d_{\text{BN}}\simeq 300\>$nm (see Table
S1). Its dielectric constant is written as (see Section IV for a similar
derivation in the case of SrTiO3 samples) :

$\epsilon_{\text{r}}=\epsilon_{\text{SiO}_{2}}\left(1+\frac{d_{\text{BN}}}{d_{\text{SiO}_{2}}}\frac{\epsilon_{\text{SiO}_{2}}}{\epsilon_{\text{BN}}}\right)^{-1}$
(8)

For $\epsilon_{\text{SiO}_{2}}=3.9$ and $\epsilon_{\text{BN}}\simeq 3.6$, we
obtain $\epsilon_{\text{r}}\simeq 3.7$. The back-gate capacitance is thus :

$C_{\text{g}}=\frac{\epsilon_{0}\epsilon_{\text{r}}}{d}\simeq
10.9\>\text{nF.cm}^{-2}$ (9)

From several tunneling conductance gate maps we compute the slopes of the
isodensity lines around charge-neutrality to be of the order of 0.1 (the tip
capacitance changes with the tip apex shape). We then obtain an estimation of
the tip capacitance as :

$\frac{C_{\text{g}}}{C_{\text{tip}}}\approx 0.1\Rightarrow
C_{\text{tip}}\approx 110\>\text{nF.cm}^{-2}$ (10)

We now focus on the tip-induced gating at $\nu=0$ at the bias positions of the
two peaks of the split LL0, positions that we used to estimate the $\nu=0$ gap
$\Delta E^{0}(B)$. The tip-induced gating at those biases is given by :

$\left|\delta n_{\text{tip}}\right|=\frac{C_{\text{tip}}}{e}\frac{\Delta
E^{0}(B)}{2}$ (11)

which yields for the tip-induced variation of the filling factor ($\nu=0$ at
$V_{\text{b}}=0$, set by the gate voltage) :

$\left|\delta\nu_{\text{tip}}\right|=\frac{\phi_{0}C_{\text{tip}}\Delta
E^{0}(B)}{2eB}$ (12)

Numerically, at $B=14\>$T, we measured $\Delta E^{0}\simeq 55\>$meV (see Fig.
2e from the main text), which gives $\left|\delta\nu_{\text{tip}}\right|\simeq
0.06$. At $B=1.5\>$T, we measured $\Delta E^{0}\simeq 15\>$meV, which gives
$\left|\delta\nu_{\text{tip}}\right|\simeq 0.14$. In both cases, the tip-
induced variation of the filling factor is negligible and we can assume it
does not alter the physics of the $\nu=0$ state we are probing. Note also that
since $\Delta E^{0}$ scales as the Coulomb energy
$\mathcal{E}_{\text{C}}\propto\sqrt{B}$, this yields
$\left|\delta\nu_{\text{tip}}\right|\propto B^{-1/2}$ and the effect of the
tip on the filling factor decreases at higher magnetic fields.

###  VI.2 SrTiO3 samples

The back-gate insulator is here comprised of the SrTiO3 substrate and the hBN
flake. Its thickness is $d\simeq d_{\text{STO}}=500\>\mu$m and its dielectric
constant $\epsilon_{\text{r}}$ can be determined as explained in Section IV.
From Fig. S2b we can estimate $\epsilon_{\text{r}}\approx 7\>500$ at charge-
neutrality for sample STO07. We then compute the back-gate capacitance as :

$C_{\text{g}}=\frac{\epsilon_{0}\epsilon_{\text{r}}}{d}\simeq
13.3\>\text{nF.cm}^{-2}$ (13)

We obtain from several tunneling conductance gate maps a similar ratio of the
capacitances of the order of :

$\frac{C_{\text{g}}}{C_{\text{tip}}}\approx 0.1\Rightarrow
C_{\text{tip}}\approx 133\>\text{nF.cm}^{-2}$ (14)

Using Equation (12), we estimate at $B=14\>$T where we measured $\Delta
E^{0}\simeq 40\>$mV a tip-induced variation of the filling factor around
$\nu=0$ of $\left|\delta\nu_{\text{tip}}\right|\simeq 0.05$. At $B=6\>$T where
$\Delta E^{0}\simeq 15\>$meV, we have
$\left|\delta\nu_{\text{tip}}\right|\simeq 0.04$. The effect of the tip on the
filling factor is therefore also negligible.

##  VII Contact potential difference measurements

We show here contact potential difference $V^{\text{CPD}}$ measurements
performed in AFM mode on the graphene as a function of the tip-sample
distance, which enables us to extract the exact band-bending beneath the tip
due to the work-function difference between the tip and graphene. We start
from the contact in AFM mode and then retract the tip at a position
$z_{\text{scan}}$. The frequency shift $\Delta f$ of the tuning fork is
measured as the sample voltage $V_{\text{b}}$ is swept from $-5\>$V to $+5\>$V
while the tip is kept grounded. We obtain a parabola described by the
following equation :

$\Delta f=C\left(V_{\text{b}}-V^{\text{CPD}}\right)+\Delta f_{\text{min}}$
(15)

such as the one shown in Fig. S5a. $V^{\text{CPD}}$ is estimated as the
position of the maximum of the parabola, using a quadratic function to fit the
curve. The variation of $V^{\text{CPD}}$ with the tip-sample distance is shown
in Fig. S5b. A linear dependence is observed, and the intrinsic
$V^{\text{CPD}}_{0}$ between the tip and graphene can be determined by the
intercept of the linear fit with $z=0$, where $z$ is the absolute distance
between the tip and graphene, given by $z=z_{\text{scan}}-z_{0}$, with $z_{0}$
the tip-sample contact position.

Fig. S5: Contact potential difference measurements on sample AC04. a,
Frequency shift $\Delta f$ as a function of sample bias $V_{\text{b}}$ (blue
points) at a tip-sample distance of $z=115\>$nm from graphene.
$V^{\text{CPD}}(z)$ is determined as the maximum of the fitting parabola in
dashed yellow. b, $V^{\text{CPD}}$ as a function of tip-sample distance $z$
from graphene (blue points). The CPD between the tip and graphene is
determined by the intercept of the linear fit (black line) with $z=0$. c,
Evolution of the curvature $C$ of the parabolas as a function of the tip-
sample distance $z$ from graphene (blue points). The fit using Equation (16)
(black line) yields the tip-sample contact position $z_{0}$, from which we
compute the tip-sample distance $z=z_{\text{scan}}-z_{0}$.

In order to compute $z_{0}$, we fit in Fig. S5c the evolution of the curvature
coefficient $C$ of the parabola with the tip-sample distance using the
following equation, see Ref. [12] :

$C=\frac{\beta}{\left(z_{\text{scan}}-z_{0}\right)^{2}}$ (16)

with $z_{0}$ and $\beta$, the force probe calibration parameter, as the two
fitting parameters. The obtained value of $z_{0}$ enables us to plot Figs. S5b
and c as a function of the absolute tip-sample distance $z$. From Fig. S5b we
obtain a negligible contact potential difference of
$V^{\text{CPD}}_{0}=0.02\>$V for a tip-sample distance in tunnel mode, that
is, of the order of $1\>$nm. This very small potential can be accounted for by
the fact that the tip apex is covered with gold due to the tip reshaping on
the gold contact, and is not made of PtIr unlike the bulk of the wire. The
work-function of gold is about $4.8\>$eV [13] while graphene is given at
$4.6\>$eV [14]. Our measurements show that these two work-functions for our
tip apex and graphene on hBN are actually very close and thus do not induce a
significant band-bending. This must be compared to the work-function of PtIr
which is $5.7\>$eV [15], which would lead to a much more unfavorable situation
with a significant band-bending.

##  VIII Tunneling conductance gate maps and interaction-induced gap in sample
AC04

We show in Fig. S6 some examples of $\Delta E^{0}$ estimations at different
magnetic fields $B$ for the sample AC04 (Gr/hBN/SiO2) which were reported in
Fig. 2e of the main text. $\Delta E^{0}$ is computed as the maximum value of
the peak-to-peak energy between LL${}_{0^{\pm}}$ extracted from the individual
spectra of the gate map.

Fig. S6: Evolution of $\bm{\Delta E^{0}}$ with the magnetic field on hBN/SiO2,
in sample AC04. a-l, Tunneling conductance, d$I_{\text{t}}$/d$V_{\text{b}}$,
gate maps at different magnetic fields. b,c were performed at the same
position, same for h,l. m, Some d$I_{\text{t}}$/d$V_{\text{b}}$ spectra,
centered around zero bias, used to estimate $\Delta E^{0}$ from
d$I_{\text{t}}$/d$V_{\text{b}}$ gate maps at different magnetic fields $B$.

##  IX Tunneling conductance gate maps and interaction-induced gap in sample
STO07

Fig. S7: Tunneling conductance gate maps as a function of the magnetic field
$\bm{B}$ on hBN/SrTiO3, in sample STO07 and at the same position.

Figure S7 shows the evolution of the tunneling conductance
d$I_{\text{t}}$/d$V_{\text{b}}$ gate maps for decreasing magnetic fields from
$B=14\>$T to $B=1\>$T in sample STO07 (SrTiO3 substrate). For each magnetic
field, we perform the same gate sweep from $V_{\text{g}}=80\>$V to
$V_{\text{g}}=-20\>$V in order to keep the hysteresis cycle of the SrTiO3
substrate constant. The full d$I_{\text{t}}$/d$V_{\text{b}}$ gate map at
$B=14\>$T is shown in panel (a), and the zoom on the yellow rectangle centered
on the LL0 gap is displayed in panel (b). The next panels are zooms on the
same area as (b). All those d$I_{\text{t}}$/d$V_{\text{b}}$ gate maps were
acquired at the same position on graphene, up to the magnetic field drift.

Let us start at $B=1\>$T in the bottom right panel. Many peaks dispersing
negatively with the gate voltage are visible, among them we can already
distinguish LL0 and LL1. With increasing magnetic fields, other resonant peaks
eventually merge with each other to form LLs, as described in Ref. [16]. For
instance, LL2 is formed at $B=3\>$T while LL-1 becomes distinguishable at
$B=5\>$T. We also notice at $B=1\>$T that there is no pinning effect of
$E_{\text{F}}$ inside LLs, which thus disperse continuously with the gate
voltage. LLs start to pin the Fermi level at $B=3\>$T with the formation of a
small plateau for LL0 at $V_{\text{g}}=27.5\>$V. However note that the
splitting of LL0 in panel (m) is mostly due to the lifting of the orbital
degeneracy [17], such that the apparent gap at zero bias (indicated by the
green arrow) may be different from the interaction-induced gap we are aiming
for. This orbital splitting is maximum at $B=4\>$T and then decreases at
higher magnetic fields.

The gap $\Delta E^{0}$ finally opens at $B=5\>$T, see the yellow arrow in
panel (k). Since the density of states of LL0 grows with $B$, the Fermi level
stays pinned inside LL0 for a wider range of gate voltage with increasing $B$,
and as a result the gap develops on a larger LL0 plateau. At $B=8\>$T the gap
features a maximum due to its enhancement by exchange interactions, see the
red arrow in panel (h). At $B=10\>$T, a second maximum appears on the right
extremity of the LL0 plateau, see the blue arrow in panel (f), while at
$B=12\>$T we clearly distinguish two lobes marked by the white arrows in panel
(d). The left lobe corresponds to the opening of the $\nu=0$ gap whereas the
right one is due to the opening of the $\nu=1$ gap (however the $\nu=-1$ gap
is not visible).

##  X Asymmetry of the KB pattern.

We show in Fig. S8a a STM image of an asymmetric KB pattern. Fig. S8b displays
the corresponding 2D Fourier transform (2D-FT). The 2D-FT is mainly comprised
of three hexagons, defined by the yellow, red and blue encircled peaks. To
decipher the 2D-FT, we filter the STM image by considering certain peaks only.
The yellow peaks alone yield the usual honeycomb lattice in Fig. S8c. The red
peaks give the image shown in Fig. S8d, which features a triangular lattice.
When we superimpose the KB lattice drawing, we notice that each bright point
of the triangular lattice in Fig. S8d falls either on the strong white bonds
of the Kekulé lattice or at the center of the hexagons devoid of strong bond:
the addition of both images yields the bond-density wave shown in Fig. S8e
where we have filtered the STM image by considering this time both yellow and
red peaks and mostly recovered the original KB pattern. This also justifies
why the hexagon devoid of strong bond in the KB pattern appears brighter than
the neighboring hexagons comprised of three strong bonds, similarly visible in
Fig. 4a of the main text. Note that the presence of two red peaks with halved
amplitude in one direction is responsible of the slight asymmetry that is
already visible in Fig. S8e.

Fig. S8: 2D-FT decomposition of the asymmetric Kekulé distortion. a, $3\times
3\>$nm2 image showing an asymmetric KB pattern, measured at $B=14\>$T and
$V_{\text{b}}=2\>$mV. b, 2D Fourier transform (2D-FT) of the STM image in a.
Yellow circles indicate peaks of the honeycomb lattice, red and blue circles
indicate peaks of the bond-density wave. c-h, Filtered images obtained by
considering certain peaks of the FFT as indicated in the top right corner of
each panel. The Kekulé lattice is drawn in white for reference. The KB order
is mostly retrieved by considering only the yellow and red peaks. The
asymmetry of the KB pattern is encoded in the blue peaks whose two of them are
twice as high as the others due to the K-CDW order. The scale bar is the same
for all figures : $500\>$pm.

We show in Fig. S8f the image obtained after filtering using only the blue
peaks. We observe a strongly asymmetric triangular lattice encoding the Kekulé
spatial modulation at $\sqrt{3}$ times the graphene lattice parameter. The
asymmetry arises from a large asymmetry between the blue peaks in the 2D-FT,
where two peaks in one direction are twice as high as the others. This yields
dissimilar weights to the bond-density wave, as shown in Fig. S8g where we
have filtered considering red and blue peaks, and explains the strong
asymmetry we observe in the KB pattern, which is fully recovered in Fig. S8h
where we have filtered with the yellow, red and blue peaks. We conjecture that
this strong asymmetry of the 2D-FT originates from the existence of the K-CDW
order whose contribution is visible in Fig. S8f, since in symmetric KB pattern
(where this K-CDW order is not visible) there is no such asymmetry between the
blue and red peaks.

##  XI Bias dependence of the KB order.

We show in Fig. S9 constant-height mode STM images where we have changed the
sample bias $V_{\text{b}}$ during scanning. The red and blue arrows on the
right of each image show the direction of the slow scan axis, and their color
corresponds to the actual sample bias which is indicated in the bottom
figures. In Fig. S9a we clearly observe a contrast inversion when switching
the sample bias from LL${}_{0^{+}}$ to LL${}_{0^{-}}$, with the continuity of
the KB pattern at the interface. In Fig. S9b we see the transition from the
usual honeycomb lattice to the KB pattern when switching the sample bias from
LL1 to LL${}_{0^{+}}$. This indicates that both the KB order and the
underlying K-CDW disappear when the bias voltage is away from the zLL peaks,
which rules out a tip-induced artifact as the origin of the observed KB order.

Fig. S9: Contrast inversion and emergence of the Kekulé bond order. $3\times
3\>$nm2 STM images during which we changed the bias voltage as shown in the
bottom insets (the current color bars are tuned separately for each half of
the images). a, We start (bottom) at $V_{\text{b}}=32\>$mV (LL${}_{0^{+}}$)
and switch (top) to $V_{\text{b}}=-12\>$mV (LL${}_{0^{-}}$) to observe a
contrast inversion of the KB lattice. b, We start (top) at
$V_{\text{b}}=200\>$mV (LL1) and switch (bottom) to $V_{\text{b}}=20\>$mV
(LL${}_{0^{+}}$) and observe the emergence of the KB order from the honeycomb
lattice. Scale bars for both images are $500\>$pm.

##  XII K-CDW configuration change.

Fig. S10: Asymmetry reversal of the Kekulé pattern. $3\times 3\>$nm2 STM
images measured at $B=14\>$T, $V_{\text{b}}=2\>$mV and at the same position.
The three images were measured successively (scanning time : 1 min). A jump
occurs in b at the scan line indicated by the red arrows, leading to an
inversion of the asymmetry of the Kekulé pattern. The slow scan axis direction
is indicated by the blue arrows on the left of each image. Scale bars for the
three images are $500\>$pm.

We illustrate here the time-varying nature of the K-CDW. We show in Fig. S10
three successive images acquired in a row at the same position and
$V_{\text{b}}=2\>$mV. The vertical arrow on the left of each image indicates
the direction of the slow scan axis. Figure 4d belongs to the same set of
image acquisition. The lattice in overlay describes the asymmetric KB pattern,
with the white links being the strong bonds of the KB order, whereas the
asymmetry that comes from the K-CDW order makes the hexagons with blue weak
bonds brighter than the hexagons with red weak bonds. The next image in panel
(b) (duration of each image : 53 seconds) starts from the bottom, where we
observe the same KB pattern. However, a jump occurs at the line indicated by
the red arrows, and, after that, in the top part of the image, the asymmetry
of the KB pattern is reversed. Using the lattice in overlay as a guide for the
eye, we see that the red hexagons are brighter (due to the three strong white
bonds almost merging together), such that the new pattern is the mirror of the
previous one. Eventually, the next image in panel (c) displays this new
pattern with brighter red hexagons on the whole area, and the next images we
realized during several minutes happened to be identical. This indicates that
the K-CDW order transited and reversed the asymmetry of the KB phase. Note
that the pattern of the strong white bonds, which defines the KB order, stays
unchanged in the three images (in opposition to the KB order transition shown
in Figs. 4f and g).

Fig. S11: Change of the Kekulé asymmetry. $10\times 10\>$nm2 STM images
measured at $B=14\>$T and $V_{\text{b}}=25\>$mV. In b, the asymmetry pattern
changes at the scan line indicated by the red arrows. Scale bars for both
images are $1\>$nm.

The asymmetry reversal of the KB order due to the K-CDW transition is well
seen in the 2D-FT of both images in Figs. S10a and c, see top insets. For Fig.
S10a, the K-CDW appears in the inner hexagon, where the two peaks encircled in
green are twice as bright as the other four peaks. On the contrary, for Fig.
S10c, this is now the yellow peaks that are brighter than the other four, with
the amplitude of the green peaks lowered. The change of the direction of the
two brighter peaks induces the change of the asymmetry pattern of the KB
order. Interestingly, the outermost hexagon, which corresponds to the bond-
density wave, also features a change in the intensity of its peaks : in Fig.
S10a the blue peaks are halved in amplitude while in Fig. S10c the red peaks
are halved. This does not yield any significant change of the KB pattern but
this may mean that the bond-density wave and the K-CDW are entangled.

Therefore, the asymmetry of the KB patterns we observed depends on the K-CDW
order and its fluctuations. Moreover, Fig. S10b shows an image in which the
K-CDW switched from that of Fig. S10a to that of Fig. S10c during the
acquisition. This change that occurred during the scan indicates that the
K-CDW switches instantaneously on the time-scale of the scan speed. This can
also indicates either a change of the K-CDW on the entire sample in case of a
homogeneous K-CDW, or the displacement of domains with different K-CDW
configurations separated by domain walls. Note that this only concerns this
anomalous K-CDW, which coexists with the KB order, the latter being unchanged
in the three images (the bright bonds pattern remains the same).

Such transitions of the K-CDW happened a few times during our measurements. In
Fig. S11a we show a $10\times 10\>$nm2 image of an asymmetric KB pattern with
the circle-like pattern formed by the merging of the strong bonds inside one
hexagon of the KB order unit cell. Imaging the same area a few minutes later
in Fig. S11b unveils a spontaneous change of K-CDW configuration, similar to
that in Fig. S10b, which occurred during scanning on the line indicated by the
red arrows: in the top part, the three strong bonds merge together inside
another hexagon of the KB order unit cell with respect to the bottom part (see
the white dashed line which intercepts the circles in the top part of the
image and, conversely, passes between the circles in the bottom part). As
previously, the KB order lattice itself does not change.

We point out that we cannot exclude a K-CDW configuration change induced by
the action of the scanning tip. Still, such a tip-induced change also implies
that the K-CDW is not pinned and can be subject to fluctuations.

##  XIII Induced $\mathbf{t_{2}}$ asymmetry in the charge-density-wave state.

We discuss here how a second nearest-neighbour hopping asymmetry gaps the
zeroth Landau level (zLL) of graphene. We consider the spinor
$\psi=(\psi_{AK},\psi_{BK},\psi_{AK^{\prime}},\psi_{BK^{\prime}})$, where
$\psi_{\sigma\tau}$ is a zLL single-particle wavefunction in sublattice
$\sigma$ and valley $\tau$. In this basis, both the sublattice imbalance
$\Delta n=n_{A}-n_{B}$ and the second nearest-neighbor hopping asymmetry
$\Delta t_{2}=t_{2A}-t_{2B}$ (see Fig. S12) enter the low-energy Hamiltonian
close to the Dirac point with the matrix $\tau_{0}\otimes\sigma_{z}$ in valley
($\tau$) and sublattice space ($\sigma$). This matrix structure implies that
both perturbations gap out the $K$ and $K^{\prime}$ points of graphene, with a
gap given by [30] :

$E_{\text{g}}=\Delta n+\dfrac{3}{2}\Delta t_{2},$ (17)

which is of equal sign for both valleys. We can visualize the effect of
$\Delta t_{2}$ and $\Delta n$ on the zeroth Landau level by diagonalizing the
graphene Hamiltonian in the presence of a magnetic field. The spectrum with
$\Delta t_{2}\neq 0$ and $\Delta n\neq 0$ are shown in Figs. S12a and b
respectively, obtained with the kwant package [18]. We can confirm numerically
that the gap is given by Equation (17) and that when $\Delta
n=-\dfrac{3}{2}\Delta t_{2}$ the gap closes, confirming that both
perturbations enter the Hamiltonian with the same matrix structure.

Fig. S12: Effect of sublattice charge imbalance and a $\bm{t_{2}}$ asymmetry
on the zeroth Landau level. a shows that the effect of a finite charge
imbalance $\Delta n=n_{A}-n_{B}$ is to gap the zeroth Landau Level of
graphene. b shows that a hopping asymmetry $\Delta t_{2}=t_{2,A}-t_{2,B}$ also
opens up a gap, that depends on momentum $k$ as we move away from the $K$ and
$K^{\prime}$ points. The parameters are chosen so that $E_{\text{g}}$ is the
same on both plots at the $K$ and $K^{\prime}$ points, according to Equation
(17). Simulations were performed using the kwant software [18] for a $41\times
41$ hexagonal lattice with $\phi=0.003$ flux per plaquette, in units of the
flux quantum. Energies are measured in units of the nearest-neighbor hopping
$t$. For a, $\Delta t_{2}=0$ and $\Delta n=0.045$ , while for b, $\Delta
t_{2}=0.015$ and $\Delta n=0$.

The above argument suggests that interactions that induce a finite $\Delta n$
will generically induce a finite $\Delta t_{2}\neq 0$, as they both enter with
the same matrix structure. To exemplify this generic behavior we use the
Hamiltonian of graphene in the presence of nearest-neighbor interactions
$V_{1}$ :

$H=-t\sum_{\left\langle
ij\right\rangle}\left(c^{\dagger}_{i}c_{j}+\mathrm{h.c.}\right)+V_{1}\sum_{\left\langle
ij\right\rangle}n_{i}n_{j},$ (18)

where the sums are taken over nearest neighbors of the honeycomb lattice. Note
in particular that the second-nearest neighbor hopping is explicitly zero in
the Hamiltonian. In the limit of infinitely large interaction, the ground
state of $H$ at half-filling is a charge-density wave with one fully occupied
and one fully empty sublattice, a state characterized by $\Delta n=1$. The
bond asymmetry $\Delta t_{2}$ is expected to be exactly zero in this limit,
since all sites on one sublattice are completely full and thus no states are
available to hop to. Similarly, all sites on the other sublattice are
completely empty such that no states are available to hop from. At
sufficiently large (but finite) $V_{1}/t$ the ground state is a charge-density
wave with partial sublattice imbalance, as we numerically show in Fig. S13. In
Fig. S13a, we show the expectation value $\Delta n=\left\langle
c^{\dagger}_{iA}c_{iA}-c^{\dagger}_{iB}c_{iB}\right\rangle$, where $c_{i\tau}$
and $c^{\dagger}_{i\tau}$ act on unit cell $i$ of sublattice $\tau$, in the
ground state of $H$ for different $V_{1}$ obtained by using the infinite
density matrix renormalization group (iDMRG), implemented using the tenpy
package [19], as explained in Ref. [28]. For small interactions $\Delta n$ is
close to zero, and grows continuously to one as $V_{1}$ is increased [28],
signaling a second-order phase transition (see e.g. Ref. [28] for a
discussion). As shown in Fig. S13b, we observe a concomitant second nearest-
neighbor bond asymmetry, defined as $\Delta
t_{2}=\mathrm{Re}\left[\left\langle
c^{\dagger}_{iA}c_{i+1A}-c^{\dagger}_{iB}c_{i+1B}\right\rangle\right]$, that
develops at intermediate values of $V_{1}$, as expected based on our previous
symmetry discussion. As $V_{1}$ increases, $\Delta t_{2}$ increases until
reaching a maximum, and then decreases as $V_{1}$ becomes larger, for all
cylinder circumferences $L_{y}=6,8,10$. It is possible to check numerically
that setting $t=0$ in Equation (18), i.e. in the limit $V_{1}/t\to\infty$,
leads to $\Delta t_{2}=0$ and $\Delta n=1$, as discussed above.

The above results support that a charge-density-wave order with a partial
sublattice imbalance, i.e. $0<\Delta n<1$, is generically accompanied by a
second-nearest-neighbor bond asymmetry, $\Delta t_{2}\neq 0$, as argued in the
main text.

Fig. S13: Induced $\bm{\Delta t_{2}}$ asymmetry by interactions. a shows that
a sublattice charge imbalance $\Delta n=\left\langle
n_{A}-n_{B}\right\rangle\neq 0$ develops as $V_{1}$ increases. b shows the
concomitant emergence of a second nearest-neighbor bond asymmetry $\Delta
t_{2}=\left\langle t_{2A}-t_{2B}\right\rangle\neq 0$, peaking at intermediate
values of $V_{1}$. The simulations are carried out for cylinder circumferences
of $L_{y}=6,8,10$ sites, all with bond-dimension $\chi=1000$, using the tenpy
package [19].

In the zLL of graphene the wavefunctions at each valley live in different
sublattices and thus a full valley polarization implies a full sublattice
polarization $\Delta n=1$, in which case $\Delta t_{2}=0$. When the sublattice
polarization is not maximal then a finite $\Delta n$ and $\Delta t_{2}$ are
expected, consistent with what is observed in experiment (see Fig. 5). As
mentioned in the main text, this effect can originate from Landau level mixing
since the sublattice index is not locked to valley index beyond the zLL
[31-33].

##  XIV Influence of the Moiré superlattice on the CDW phase.

Fig. S14: Disappearance of the charge-density wave at low magnetic field in
sample AC23. a, CDW at $B=9\>$T. b, CDW at $B=7\>$T. The Moiré superlattice is
visible but does not perturb the CDW pattern. c,d, Honeycomb lattice with no
CDW at $B=4\>$T. e, Honeycomb lattice at $B=4\>$T with residual traces of CDW,
see the zoom in f of the white rectangle. Scale bars for all figures are
$500\>$pm.

The sample AC23 displays a weak Moiré superlattice (weak in the sense that is
not always visible in our images). This rises the question of whether the CDW
phase we observed was induced by the Moiré pattern, which could also break the
sublattice symmetry, or not. In such case, we may expect the CDW pattern to
rely on that of the Moiré, with the sublattice polarization depending on the
position inside the Moiré superlattice (due to the periodic potential it
induces in graphene). Fig. S14b displays a CDW phase observed at $B=7\>$T. The
Moiré pattern is barely visible but appears as bright spots, such as the ones
indicated by dashed white circles. The CDW pattern itself is homogeneous over
the whole Moiré lattice.

STM images of the similar sample AC24 on SrTiO3 did not exhibit any Moiré
superlattice. Still, in the same conditions, at charge neutrality, we observed
signatures of a CDW phase, shown in Fig. S15b, which indicates that the CDW we
observe in our hBN/SrTiO3 samples is indeed an intrinsic consequence of many-
body interactions at charge neutrality and not due to extrinsic substrate-
induced sublattice symmetry breaking.

##  XV Disappearance of the CDW order at low magnetic field.

Figures S14a and b show that the CDW persists at $B=9\>$T and $B=7\>$T on
sample AC23. However, Figs. S14c to e at $B=4\>$T does not display the CDW
anymore but only the usual sublattice-unpolarized honeycomb lattice.

Fig. S15: Charge-density-wave order in sample AC24. a, Honeycomb lattice at
$B=14\>$T and $V_{\text{b}}=-350\>$mV observed at $\nu=0$. b, CDW under the
same conditions but at $V_{\text{b}}=-18\>$mV. Scale bars for both figures are
$500\>$pm.

Finally, it is theoretically expected [10] that graphene undergoes a first-
order phase transition from the CDW to the F phase. Such a transition should
induce the formation of domains in graphene with the coexistence of both
phases around the magnetic field at which the transition occurs. When taking a
closer look to Fig. S14e taken at $B=4\>$T, see its zoom in Fig. S14f, one can
discern some bright dots in some parts of the image (see the blue circles).
This residual asymmetry of the honeycomb lattice is reminiscent of the charge-
density wave. It is possible that it may constitute a signature of such
domains around the transition between the CDW and F phases.

## References

  * [1] Wang, L. _et al._ One-Dimensional Electrical Contact to a Two-Dimensional Material.  _Science_ 342, 614–617 (2013).
  * [2] Li, X.-X. _et al._ Gate-controlled reversible rectifying behaviour in tunnel contacted atomically-thin MoS2 transistor.  _Nature Communications_ 8, 970 (2017).
  * [3] Choi, Y. _et al._ Electronic correlations in twisted bilayer graphene near the magic angle.  _Nature Physics_ 15, 1174–1180 (2019).
  * [4] Sakudo, T. & Unoki, H.  Dielectric Properties of SrTiO3 at Low Temperatures.  _Phys. Rev. Lett._ 26, 851–853 (1971).
  * [5] Hemberger, J., Lunkenheimer, P., Viana, R., Böhmer, R. & Loidl, A.  Electric-field-dependent dielectric constant and nonlinear susceptibility in SrTiO3.  _Phys. Rev. B_ 52, 13159 (1995).
  * [6] Sachs, R., Lin, Z. & Shi, J.  Ferroelectric-like SrTiO3 surface dipoles probed by graphene.  _Scientific Reports_ 4, 3657 (2014).
  * [7] Chen, S., Chen, X., Duijnstee, E. A., Sanyal, B. & Banerjee, T.  Unveiling Temperature-Induced Structural Domains and Movement of Oxygen Vacancies in SrTiO3 with Graphene.  _ACS Appl. Mater. Interfaces_ 12, 52915–52921 (2020).
  * [8] Das Sarma, S., Hwang, E. H. & Tse, W.-K.  Many-body interaction effects in doped and undoped graphene: Fermi liquid versus non-fermi liquid.  _Phys. Rev. B_ 75, 121406(R) (2007).
  * [9] Luican, A., Li, G. & Andrei, E. Y.  Quantized Landau level spectrum and its density dependence in graphene.  _Phys. Rev. B_ 83, 041405(R) (2011).
  * [10] Chae, J. _et al._ Renormalization of the Graphene Dispersion Velocity Determined from Scanning Tunneling Sprectroscopy.  _Phys. Rev. Lett._ 109, 116802 (2012).
  * [11] Ren, Y.-N. _et al._ Spatial and magnetic confinement of massless Dirac fermions.  _Phys. Rev. B_ 104, L161408 (2021).
  * [12] Laurent, J., Sellier, H., Mosset, A., Huant, S. & Chevrier, J.  Casimir force measurements in Au-Au and Au-Si cavities at low temperature.  _Phys. Rev. B_ 85, 035426 (2012).
  * [13] Anderson, P. A.  Work Function of Gold.  _Phys. Rev._ vol. 115, 553 (1959).
  * [14] Yan, R. _et al._ Determination of graphene work function and graphene-insulator-semiconductor band alignment by internal photoemission spectroscopy.  _Appl. Phys. Lett._ 101, 022105 (2012).
  * [15] Morán-Meza, J. A. _et al._ Reverse electrochemical etching method for fabricating ultra-sharp platinum/iridium tips for combined scanning tunneling microscope/atomic force microscope based on a quartz tuning fork.  _Current Applied Physics_ vol. 15, 1015–1021 (2015).
  * [16] Gutiérrez, C. _et al._ Interaction-driven quantum Hall wedding cake-like structures in graphene quantum dots.  _Science_ 361, 789–794 (2018).
  * [17] Luican-Mayer, A. _et al._ Screening Charged Impurities and Lifting the Orbital Degeneracy in Graphene by Populating Landau levels.  _Phys. Rev. Lett._ 112, 036804 (2014).
  * [18] Groth, C. W., Wimmer, M., Akhmerov, A. R. & Waintal, X.  Kwant: a software package for quantum transport.  _New Journal of Physics_ 16, 063065 (2014).
  * [19] Hauschild, J. & Pollmann, F.  Efficient numerical simulations with Tensor Networks: Tensor Network Python (TeNPy).  _SciPost Phys. Lect. Notes_ 5 (2018).