# Landau-Level Mixing and SU(4) Symmetry Breaking in Graphene

Nemin Wei  Department of Physics, Yale University, New Haven, CT 06520, USA
Guopeng Xu  Department of Physics and Astronomy, University of Kentucky,
Lexington, Kentucky 40506-0055, USA  Inti Sodemann Villadiego  Institut fur
Theoretische Physik, Universitat Leipzig, Bruderstrabe 16, 04103, Leipzig,
Germany  Chunli Huang  Department of Physics and Astronomy, University of
Kentucky, Lexington, Kentucky 40506-0055, USA

(August 28, 2024)

###### Abstract

Recent scanning tunneling microscopy experiments on graphene at charge
neutrality under strong magnetic fields have uncovered a ground state
characterized by Kekulé distortion (KD). In contrast, non-local spin and
charge transport experiments in double-encapsulated graphene, which has a
higher dielectric constant, have identified an antiferromagnetic (AF) ground
state. We propose a mechanism to reconcile these conflicting observations, by
showing that Landau-level mixing can drive a transition from AF to KD with the
reduction of the dielectric screening. Our conclusion is drawn from studying
the effect of Landau-level mixing on the lattice-scale, valley-dependent
interactions to leading order in graphene’s fine structure constant
$\kappa=e^{2}/(\hbar v_{F}\epsilon)$. This analysis provides three key
insights: 1) Valley-dependent interactions remain predominantly short-range
with the $m=0$ Haldane pseudopotential being at least an order of magnitude
greater than the others, affirming the validity of delta-function
approximation for these interactions. 2) The phase transition between the AF
and KD states is driven by the microscopic process in the double-exchange
Feynman diagram. 3) The magnitudes of the coupling constants are significantly
boosted by remote Landau levels. Our model also provides a theoretical basis
for numerical studies of fractional quantum Hall states in graphene.

Introduction: The rich interplay between interactions, topology, and the
approximate spin-valley SU(4) symmetry in graphene quantum Hall physics
continues to spark puzzles even after two decades of research Halperin and
Jain (2020); Alicea and Fisher (2006); Nomura and MacDonald (2006); Yang  _et
al._ (2006); Checkelsky  _et al._ (2009); Nomura  _et al._ (2009); Dean  _et
al._ (2011); Goerbig (2011); Sodemann and MacDonald (2014); Spanton  _et al._
(2018); Zibrov  _et al._ (2018); Zhou  _et al._ (2020); Veyrat  _et al._
(2020); Farahi  _et al._ (2023) and could provide insights into more complex
correlated states such as those emerging in moiré superlattice flatbands Lu
_et al._ (2023); Sharpe  _et al._ (2019); Serlin  _et al._ (2020). A recent
notable discovery is the Kekulé distorted (KD) state identified by scanning
tunneling microscopy (STM) experiments in charge neutral monolayer graphene Li
_et al._ (2019); Liu  _et al._ (2022); Coissard  _et al._ (2022), which is in
tension with the antiferromagnetic (AF) state supported by earlier non-local
spin and charge transport experiments Zhou  _et al._ (2022); Wei  _et al._
(2018); Stepanov  _et al._ (2018); Young  _et al._ (2014); Fu  _et al._
(2021); Wei  _et al._ (2021). We present a model to reconcile these diverging
observations.

Kharitonov Kharitonov (2012) proposed that the phase diagram of neutral
graphene under a strong magnetic field can be conveniently represented in a
two-dimensional parameter space defined by $u_{\perp}$ and $u_{z}$. These
phenomenological parameters, describing the strength of valley-exchange
($u_{\perp}$) and valley-antisymmetric ($u_{z}$) interactions that breaks the
SU(4) symmetry in the zeroth Landau level (LL), are influenced by electron-
electron and electron-phonon potentials. Given their complex origins, they are
often regarded as experimental fitting parameters Zibrov  _et al._ (2018);
Zhou  _et al._ (2022); Liu  _et al._ (2022).

Figure 1: The flow of $u_{\perp}(\kappa)$ and $u_{z}(\kappa)$ with graphene’s
fine structure constant, $\kappa$, indicates a ground state transition from
the AF phase to the KD phase as $\kappa$ increases from $0$ to $2$. The dashed
lines serve as visual guides and angles in the inset are
$\arctan(u_{z}/u_{\perp})$. In the absence of Zeeman splitting, the AF-KD
phase boundary is marked by the black $135^{\circ}$ line. For comparison, we
note $\Delta_{Z}2\pi l_{B}^{2}=492\text{meV}\text{nm}^{2}$ where
$\Delta_{Z}=2\mu_{B}B$ is the Zeeman energy.

In this study, we present a comprehensive analysis of the effect of Landau
level mixing on the $(u_{\perp},u_{z})$ to leading order of graphene’s fine-
structure constant, $\kappa=e^{2}/(\hbar v_{F}\epsilon)$. The constant
$\kappa$, which depends on the electric charge $e$, Planck’s constant $\hbar$,
the Fermi velocity $v_{F}$, and the adjustable environmental dielectric
constant $\epsilon$, controls the strength of the SU(4) symmetric long-range
Coulomb interactions and the extent of Landau-Level mixing, as noted in
previous studies Peterson and Nayak (2014); Sodemann and MacDonald (2014,
2013). Our study reveals that as $\kappa$ increases (i.e., Landau level mixing
becomes stronger), $u_{\perp}(\kappa)$ decreases more rapidly than
$u_{z}(\kappa)$ increases. This triggers a transition from AF to a KD, as
illustrated in Fig. 1. Thus, the variation in Landau-level mixing across
different samples leads to the emergence of either the KD or AF state. This
trend is consistent with experimental findings: double-encapsulated samples
used in transport experiments, with their larger $\epsilon$ (and consequently
smaller $\kappa$), exhibit AF order, whereas un-encapsulated samples in STM
experiments display KD order.

The Standard Model of Monolayer Graphene: At the atomic scale, interactions
between electrons in graphene exhibit weak sublattice dependence consistent
with its $C_{6v}$ symmetry. In the dilute doping regime, a specific subset of
Bloch states –those with energies below an energy cutoff $\Lambda$ and
concentrated at the Brillouin zone corners – form the basis for constructing
the long-wavelength field theory. A practical method to incorporate atomic-
scale interactions into this long-wavelength theory is to model them as Dirac
delta functions with corresponding bare coupling constants computed from the
three dimensional atomic wavefunctions. This methodology leads to the
formulation of what we refer to as ”graphene’s standard model”, encompassing
the band Hamiltonian $H_{0}$, the SU(4) symmetric long-range Coulomb
interaction $H_{s}$, and the lattice-scale interactions $H_{a}$: Lemonik  _et
al._ (2012); Aleiner  _et al._ (2007)

$H_{0}=\sum_{i=x,y}v_{F}\int\psi^{\dagger}(\bm{r})\hat{\Pi}_{i}\sigma^{i}\psi(\bm{r})\,d^{2}{r}$
(1a)
$H_{s}=\frac{e^{2}}{2\epsilon}\iint\frac{\psi^{\dagger}(\bm{r})\psi^{\dagger}(\bm{r}^{\prime})\psi(\bm{r}^{\prime})\psi(\bm{r})}{|\bm{r}-\bm{r^{\prime}}|}\,d^{2}rd^{2}r^{\prime}$
(1b) $H_{a}=\frac{1}{2}\sum_{\begin{subarray}{c}\alpha,\mu=\\\
0,x,y,z\end{subarray}}g_{\alpha\mu}^{0}\int:\left[\psi^{\dagger}(\bm{r})\tau^{\alpha}\sigma^{\mu}\psi(\bm{r})\right]^{2}:\;d^{2}r$
(1c)

Here $v_{F}$ is the Fermi velocity,
$\hat{\Pi}_{i}=-i\hbar\partial_{i}+eA_{i}(\bm{r})$ is the kinetic momentum
operator with the vector potential $A_{i}$, $\epsilon$ is dielectric constant,
and the basis $\psi=(\psi_{\uparrow},\psi_{\downarrow})^{T}$ with
$\psi_{s}=(\psi_{KAs},\psi_{KBs},\psi_{K^{\prime}Bs},-\psi_{K^{\prime}As})^{T}$
where $s$ can be $\uparrow$ or $\downarrow$. $\tau^{\alpha},\sigma^{\mu}$ are
Pauli matrices operating in the valley and sublattice subspace, respectively.

Eq. (1c) contains 16 valley-and sublattice-dependent interactions. However,
under $C_{6v}$ and lattice translation symmetry, $g_{\alpha x}^{0}=g_{\alpha
y}^{0}\equiv g_{\alpha\perp}^{0}$, $g_{x\mu}^{0}=g_{y\mu}^{0}\equiv
g_{\perp\mu}^{0}$ and we neglect $g_{00}^{0}$ since it is dwarfed by long-
range Coulomb interaction. Thus, Eq. (1c) is characterized by 8 independent
parameters or coupling-constantsAleiner  _et al._ (2007):

$g_{0z}^{0}\,,\,g_{0\perp}^{0}\,,\,g_{z0}^{0}\,,\,g_{z\perp}^{0}\,,\,g_{zz}^{0}\,,\,g_{\perp
0}^{0}\,,\,g_{\perp\perp}^{0}\,,\,g_{\perp z}^{0}.$ (2)

The coupling constants $g_{\alpha\mu}^{0}$ are influenced by electron-electron
repulsion and electron-phonon interactions, with the superscript $0$ denoting
their role as ‘bare’ inputs in the theory with cutoff $\Lambda$. A key finding
of our work is that Landau-level mixing significantly renormalizes these
atomic-scale interactions when reducing the energy from $\Lambda$ to the
cyclotron energy.

In the Supplementary Material, we provide detailed calculations of estimates
for the bare coupling constants, based on the assumption that the Bloch
functions are linear combinations of (3D) atomic $p_{z}$ orbitals. It is
important to note that convergence in real space, when summing atomic
orbitals, can be notably slow and may lead to inaccuracies Knothe and Fal’ko
(2020); sup . To circumvent these issues, performing these calculations in
momentum space has proven more effective. Our calculation suggests that the
intra-sublattice inter-valley-exchange interaction $g_{\perp\perp}^{0}$ is the
dominant lattice-scale interaction, in particular
$g_{\perp\perp}^{0}>g_{zz}^{0}$ sup ; Raines  _et al._ (2021). Upon
incorporating electron-phonon contributions, we obtain the following estimates
(constants not listed below are found to be negligible):

$(g_{\perp\perp}^{0},g_{zz}^{0},g_{\perp
z}^{0},\,g_{z\perp}^{0})=(269,184,-26,-27)\text{meV}\cdot\text{nm}^{2}.$ (3)

Figure 2: (a)-(c) Feynman diagrams representing the $\mathcal{O}(\kappa)$
corrections to two-body valley-dependent Haldane pseudopotentials. Solid lines
are fermionic propagators and the external lines are in the $n=0$ Landau
level. The curly and dashed lines stand for the long-range SU(4) Coulomb and
the contact sublattice-and valley-dependent interactions, respectively.
Swapping the curly and dashed lines in (a) and (b) does not change the
corrections to Haldane pseudopotentials and for (c) either vertex of the
dashed line can be dressed. (d)-(f) Examples of diagrams that do not
contribute to the valley anisotropic interactions. (d) and (e) vanish. (f)
yields SU(4) symmetric corrections and thus can be neglected.

Landau-level mixing and the modifications of valley-dependent interactions: In
a constant magnetic field (i.e. $\nabla\times\bm{A}=-B\hat{z}$ in Eq.(1a)),
the eigenstates of $H_{0}$ for spin-valley flavor are denoted as
$\ket{\overline{\xi n},m}$, where $\overline{\xi n}$ represents the Landau
levels $n=0,1,2...$, $\xi=\pm 1$, and $m$ enumerates the guiding center
degeneracy of each Landau level. The wavefunction is a sublattice spinor,

$\ket{\overline{\xi n}}=(\sqrt{2})^{\delta_{n,0}-1}\begin{pmatrix}\ket{n}\\\
(1-\delta_{n,0})\xi\ket{n-1}\end{pmatrix},$ (4)

where each state has a kinetic energy of $\epsilon_{\overline{\xi n}}=\xi\hbar
v_{F}\sqrt{2n}/l_{B}$. The quantum state $\ket{n}\equiv a^{\dagger
n}\ket{0}/\sqrt{n!}$ and $\ket{m}\equiv b^{\dagger m}\ket{0}/\sqrt{m!}$ are
constructed from two sets of ladder operators constructed from kinematic
momentum and guiding centers:
$a^{\dagger}=(\hat{\Pi}_{x}-i\hat{\Pi}_{y})l_{B}/\sqrt{2}$ and
$b^{\dagger}=(\hat{\eta}_{x}+i\hat{\eta}_{y})/\sqrt{2}l_{B}$ where
$\hat{\bm{\eta}}=\bm{r}-l_{B}^{2}\bm{z}\times\hat{\bm{\Pi}}$. In the symmetric
gauge, we can loosely call $m$ the angular momentum although the physical
angular momentum is $m-n$.

Since any translationally and rotationally invariant two-body interaction
Hamiltonian $\hat{V}$ can be decomposed into the Haldane pseudopotentials
using the center of mass and relative coordinates
$|NM\rangle_{c}|nm\rangle_{r}$ Haldane (1983),

$\displaystyle{}_{c}\bra{N^{\prime}M^{\prime}}{}_{r}\bra{n^{\prime}m^{\prime}}\hat{V}\ket{NM}_{c}\ket{nm}_{r}$
$\displaystyle={}_{r}\bra{n^{\prime}m^{\prime}}\hat{V}\ket{nm}_{r}\delta_{m^{\prime}-n^{\prime},m-n}\delta_{N^{\prime},N}\delta_{M^{\prime},M}.$
(5)

For two interacting electrons in the zeroth LL of our interest,
$n=n^{\prime}=0$. In the spirit of Larmor’s theorem, we can eliminate the
Zeeman field from the equation of motion in a rotating frame and recover the
spin SU(2) symmetry. Thus, the effective interactions can be represented in a
spin-independent form. Moreover, the $C_{6v}$, lattice momentum conservation,
and the valley-sublattice locking of the zeroth LL, enforce effectively an
U(1) valley conservation, leading to the Haldane pseudopotentials
$\prescript{}{r}{\bra{0m}}V\ket{0m}_{r}=V_{0}(m)+V_{a}(m)$, where $V_{0}(m)$
is spin-valley SU(4) Haldane pseudopotential that has been studied previously
Peterson and Nayak (2013) and

$\hat{V}_{a}(m)=\frac{u_{z}(m)}{2}\tau_{i}^{z}\tau_{j}^{z}+\frac{u_{\perp}(m)}{2}(\tau_{i}^{x}\tau_{j}^{x}+\tau_{i}^{y}\tau_{j}^{y}),$
(6)

are the valley dependent pseudopotentials arising from $H_{sv}$ interactions
in Eq. (1c). At the zeroth order in interactions, the effective Hamiltonian is
$\hat{V}^{(1)}=\mathcal{P}(\hat{H}_{s}+\hat{H}_{a})\mathcal{P}$ with
$\mathcal{P}$ being the projector into the degenerate non-interacting ground-
state manifold where $n<0$ LLs are fully occupied and $n>0$ LLs empty. Because
a pair of electrons of a nonzero relative angular momentum do not experience
contact interactions sup , $u_{\alpha}^{(1)}(m)=u_{\alpha}\delta_{m,0}$, where
$u_{\alpha}=(g_{\alpha 0}^{0}+g_{\alpha z}^{0})/2\pi l_{B}^{2}$ is the
notation used in Ref. Kharitonov (2012). However, Landau level mixing
generates perturbative corrections to these Haldane pseudotentials and one of
our main goals is to compute such corrections.

The leading order correction to $\hat{V}^{(1)}$ can be derived from second-
order perturbation theory Sodemann and MacDonald (2013). Our focus here is on
the effect of dielectric screening on valley-dependent interactions, which is
captured by the following Hamiltonian, comprising one $\hat{H}_{s}$ and one
$\hat{H}_{a}$:

$\displaystyle\hat{V}_{sa}^{(2)}$
$\displaystyle=-\mathcal{P}\hat{H}_{s}\mathcal{P}_{\perp}\frac{1}{\hat{H}_{0}-E_{0}}\mathcal{P}_{\perp}\hat{H}_{a}\mathcal{P}+h.c.$
(7)

Here $E_{0}$ is the energy of degenerate noninteracting ground-states and
$\mathcal{P}_{\perp}=\mathcal{I}-\mathcal{P}$. Because three-body interactions
are forbidden in the zeroth LL of monolayer graphene by the particle-hole
symmetry Peterson and Nayak (2013) and the one-body self-energy generated by
$\hat{V}_{sa}^{(2)}$ is SU(4) symmetric, the only nontrivial terms in Eq. (7)
are the two-body interactions which can be delineated into Feynman diagrams in
Fig. 2 sup . We call the three non-vanishing Feynman diagrams (a) the
particle-particle channel, (b) the double-exchange channel, and (c) the vertex
correction. Their contributions are parametrized as follows:

$u_{\alpha}^{(2)}(m)=\frac{\kappa}{2\pi
l_{B}^{2}}\sum_{\mu=0,x,y,z}(c_{m,a}^{\mu}+c_{m,b}^{\mu}+c_{m,c}^{\mu})g_{\alpha\mu}^{0}$
(8)

Because of our choice of pseudospin components (see text below Eqs. (1)), both
valleys sharing identical orbital wave functions as per Eq. (4), the
coefficients $c_{m,a/b/c}^{\mu}$ are independent of the valley index $\alpha$.

Evaluating the Feynman diagram in Fig. 2a, we find that the valley-sublattice
locking unique to the zeroth LL of graphene,
$\sigma^{z}\ket{\overline{0}m}=\ket{\overline{0}m}$, leads to the vanishing
interacting matrix elements
$\bra{\bar{n}_{i},\bar{n}_{j}}\sigma_{i}^{+}\sigma_{j}^{-}\ket{\bar{0}_{i},\bar{0}_{j}}=0$,
and $c_{m,a}^{\perp}\equiv c_{m,a}^{x}+c_{m,a}^{y}=0$. Additionally, due to
the contact valley-dependent interaction, the particle-particle ladder diagram
only induces corrections to $u_{\alpha}^{(2)}(m=0)$, with sup

$\displaystyle\frac{c_{m=0,a}^{0,z}}{4\pi}=$
$\displaystyle-2\sum_{n=1}^{2N_{c}}\prescript{}{r}{\bra{00}}\delta(\bm{r}_{ij})\ket{nn}_{r}\prescript{}{r}{\bra{nn}}r_{ij}^{-1}\ket{00}_{r}$
$\displaystyle\quad\times\sum_{n_{i,j}=-\infty}^{\infty}\frac{\Theta(n_{i}n_{j})}{|\epsilon_{\overline{n_{i}}}+\epsilon_{\overline{n_{j}}}|}|\prescript{}{c}{\langle}0|\prescript{}{r}{\langle}n||n_{i}|,|n_{j}|\rangle|^{2},$
$\displaystyle=$
$\displaystyle\sum_{n=1}^{2N_{c}}\frac{\Gamma(n+\frac{1}{2})}{2^{n+1}n!}\sum_{n_{1}=0}^{n}\frac{\binom{n}{n_{1}}}{\sqrt{2n_{1}}+\sqrt{2(n-n_{1})}}.$
(9)

Here, $\Theta(x)$ is the Heaviside function with ${\Theta(0)=1}$ and we
measure the length in the unit of $l_{B}$ an the kinetic energy in the unit of
$\hbar v_{F}/l_{B}$ to shorten the expressions. Asymptotically,
$\Gamma(n+\frac{1}{2})/n!\sim n^{-\frac{1}{2}}$ and therefore $c_{0,a}^{0,z}$
is logarithmically divergent with respect to the LL cutoff $N_{c}$, which is
expected from the renormalization group analysis Kharitonov (2012).

The Feynman diagrams in Fig. 2b and c involve particle-hole excitations in the
intermediate states and need to be calculated by a different method. We follow
Ref. Sodemann and MacDonald (2013); Haldane (1987) and express the bare two-
body interactions as Fourier integrals, $r_{ij}^{-1}=\int
d\bm{q}e^{i\bm{q}\cdot\bm{r}_{ij}}/2\pi q$ and $\delta(\bm{r}_{ij})=\int
d\bm{q}^{\prime}e^{i\bm{q}^{\prime}\cdot\bm{r}_{ij}}/(2\pi)^{2}$. For each
$(\bm{q},\bm{q}^{\prime})$, the scattering amplitude is diagonalized in the
relative angular momentum basis $\\{\ket{m}_{r}\\}$ sup . The Haldane
pseudopotentials can therefore be expressed as an integral, with

$\displaystyle\frac{c_{m,b}^{\mu}}{4\pi}=\iint_{\bm{q},\bm{q}^{\prime}}$
$\displaystyle\sum_{n,n^{\prime}=-\infty}^{\infty}\frac{\Theta(n)-\Theta(n^{\prime})}{\epsilon_{\overline{n}}-\epsilon_{\overline{n^{\prime}}}}\mathcal{V}_{n0,0n^{\prime}}^{\mu}(\bm{q}^{\prime})\mathcal{V}_{0n^{\prime},n0}(\bm{q})$
$\displaystyle\times\prescript{}{r}{\langle
m|}e^{i(\bm{q}-\bm{q}^{\prime})\cdot\bm{\eta}_{ij}}|m\rangle_{r}+c.c.,$ (10)
$\displaystyle\frac{c_{m,c}^{\mu}}{4\pi}=\iint_{\bm{q},\bm{q}^{\prime}}$
$\displaystyle\sum_{n,n^{\prime}=-\infty}^{\infty}\frac{\Theta(n^{\prime})-\Theta(n)}{\epsilon_{\overline{n}}-\epsilon_{\overline{n^{\prime}}}}\mathcal{V}_{0n,0n^{\prime}}^{\mu}(\bm{q}^{\prime})\mathcal{V}_{0n^{\prime},0n}(\bm{q})$
$\displaystyle\times e^{i\bm{q}^{\prime}\wedge\bm{q}}\prescript{}{r}{\langle
m|}e^{-i\bm{q}^{\prime}\cdot\bm{\eta}_{ij}}|m\rangle_{r}+c.c..$ (11)

Here $\int_{\bm{q}}=\int d^{2}q/(2\pi)^{2}$ and the matrix elements for the
long-range Coulomb interaction
$\mathcal{V}_{nn^{\prime},n_{1}n_{1}^{\prime}}(\bm{q})=\frac{2\pi}{q}F_{n,n_{1}}^{0}(\bm{q})F_{n^{\prime},n_{1}^{\prime}}^{0}(-\bm{q})$
and the bare sublattice-and valley-dependent interactions
$\mathcal{V}_{nn^{\prime},n_{1}n_{1}^{\prime}}^{\mu}(\bm{q})=F_{n,n_{1}}^{\mu}(\bm{q})F_{n^{\prime},n_{1}^{\prime}}^{\mu}(-\bm{q})$
with the form factors
$F_{n^{\prime},n}^{\mu}(\bm{q})=\bra{\overline{n^{\prime}}}\sigma^{\mu}e^{i\bm{q}\cdot(\hat{z}\times\bm{\Pi})}\ket{\overline{n}}$
for $\mu=0,x,y,z$. Due to the valley-sublattice locking in the zeroth LL,
$c_{m,c}^{\perp}=0$ as in the particle-particle channel.

Table 1: Coefficients of the first-order-in-$\kappa$ corrections to the
valley-dependent interactions receive contributions from three non-vanishing
Feynman diagrams in Fig. 2a (p-p), b (DE), and c(vertex), see Eq. (8). (a)
$c_{m=0,a/b/c}^{\mu}$ are listed in this table. (b)
$c_{m}^{\mu}=c_{m,a}^{\mu}+c_{m,b}^{\mu}+c_{m,c}^{\mu}$ for $m\leq 5$. The
numbers in red depend logarithmically on the cutoff $N_{c}$ of the Landau
levels. Here, $N_{c}=160$.

(a) $c_{m=0,a}^{\mu},\ \ \ c_{m=0,b}^{\mu},\ \ \ c_{m=0,c}^{\mu}$
---
$\mu$ | p-p | DE | vertex
$0$ | $-1.6773$ | $1.5609$ | $0.3915$
$z$ | $-1.6773$ | $1.5609$ | $2.7386$
$\perp$ | 0 | $-2.5900$ | 0

(b)
---
$\mu$ m | 0 | 1 | 2 | 3 | 4 | 5
$0$ | $0.2799$ | $-0.1094$ | $-0.1283$ | $-0.0144$ | $-0.0173$ | $-0.0051$
$z$ | $2.6223$ | $0.1128$ | $-0.0627$ | $0.0149$ | $-0.0011$ | $0.0051$
$\perp$ | $-2.5900$ | $0.1064$ | $-0.0141$ | $0.0032$ | $-0.0001$ | $0.0006$

Discussion of main results:– The main results of this work, encapsulated in
Eq. (S37), (10), and (11), required numerical computations and they lead to
the following three important conclusions:

1) Table 1(a) shows the contributions of three Feynman diagrams to the $m=0$
component of $c_{m}^{\mu}(c_{m}^{\perp}\equiv c_{m}^{x}+c_{m}^{y})$. Notably,
the inter-valley channel $c_{0}^{\perp}$, which is entirely derived from the
double-exchange diagram, is found to be attractive ($c_{0}^{\perp}<0$). This
observation is significant because it implies that a repulsive valley-and
sublattice-exchange scattering $g_{\perp\perp}^{0}>0$ can induce a negative
valley-exchange Haldane pseudopotential $V_{\perp}$ in the zeroth LL. This
contrasts sharply with conventional fermions with parabolic dispersion, where
repeated exchange scattering enhances opposite-spin repulsion and suppresses
spin-singlet superconductivity near a ferromagnetic critical point Berk and
Schrieffer (1966). In graphene, we found that opposite-valley scattering
becomes attractive due to the opposite chirality of electron and hole
wavefunctions 111While our conclusion is based on double-exchange diagram, the
Berk-Schrieffer diagram consists of an infinite series of repeated exchange
diagrams. In the context of Landau levels, the second sublattice components of
wavefunctions (Eq.(4)) for intermediate electron ($\xi>0$) and hole ($\xi<0$)
have opposite signs, leading to an attractive interaction matrix elements
$\mathcal{V}_{n0,0n^{\prime}}^{\perp}<0$.

2) The red coefficients in Table I(a) highlight that these values
logarithmically increase with the increasing number of Landau levels, a
consequence of the linearly increasing density-of-states in our low-energy
theory. This is an important effect as they amplify the magnitude of lattice-
scale interactions, making them 3 to 5 times stronger than the Zeeman energy,
consistent with experimental observations sup ; Zibrov  _et al._ (2018).

3) Table 1(b) details the values of $u_{\alpha}^{(2)}(m)$ for $m=0,...,5$,
incorporating contributions from the aformentioned Feynman diagrams. Notably,
the $m=0$ component is significantly larger, by at least an order of
magnitude, than the $m\neq 0$ component. This indicates that even when contact
interactions are dressed by long-range Coulomb interactions in second-order
perturbation theory, the resulting interaction remains primarily short-ranged.
While this is evident in the particle-particle channel ($c_{m\neq
0,a}^{\mu}=0$) the particle-hole channels require further explanation. In the
particle-hole channel, the bubble diagram contribution is exactly zero leaving
the exchange scattering diagrams as the primary contributors. The bubble
diagram vanishes because the susceptibility of valley or sublattice
polarization induced by charge perturbation vanishes under $C_{6v}$ symmetry.
The exchange diagrams fail to generate long-range interactions, as they
require substantial spatial overlap of the wavefunctions. This result suggests
that the coexistence phase of the KD and AF states, as proposed in Ref. Das
_et al._ , 2022; De  _et al._ , 2023, does exist, but it occupies a very small
sliver of the phase diagram near the boundary of KD and AF Stefanidis and
Villadiego (2023a, b).

Focusing on the $m=0$ component, we obtain the final expression for the
corrections to the lattice-scale interactions to leading order in $\kappa$ are
given by:

$\displaystyle u_{\perp}(\kappa)$ $\displaystyle\equiv
u_{\perp}+u_{\perp}^{(2)}(0)=\frac{g_{\perp z}^{0}\left(1+\kappa
c^{z}_{0}\right)-\kappa|c^{\perp}_{0}|g_{\perp\perp}^{0}}{2\pi l_{B}^{2}},$
(12) $\displaystyle u_{z}(\kappa)$ $\displaystyle\equiv
u_{z}+u_{z}^{(2)}(0)=\frac{g_{zz}^{0}\left(1+\kappa
c^{z}_{0}\right)-\kappa|c^{\perp}_{0}|g_{z\perp}^{0}}{2\pi l_{B}^{2}}.$ (13)

By incorporating the estimates Eq. (3), we derive the main results shown in
Fig. 1. At $\kappa=0$ (i.e., without Landau-level mixing), $u_{\perp}$ and
$u_{z}$ are determined by attractive $K$-optical-phonon potential ($g_{\perp
z}^{0}$) and repulsive sublattice asymmetry electron-electron repulsion
($g_{zz}^{0}$) respectively, with the latter being more dominant, leading to
$u_{z}>-u_{\perp}$ and thus an AF ground state. At finite $\kappa$, the
respective contributions of $g_{\perp z}^{0}$ and $g_{zz}^{0}$ to $u_{\perp}$
an $u_{z}$ are enhanced by the same factor $1+\kappa c_{0}^{z}$. However, the
double-exchange Feynman diagram introduces into $u_{\perp}$ the on-site
valley-exchange scattering ($g_{\perp\perp}^{0}$), the strongest type of bare
lattice-scale interactions in monolayer graphene sup ; Raines  _et al._
(2021), making $u_{\perp}$ decrease faster than the increases in $u_{z}$ and
ultimately tipping the balance to $u_{z}<-u_{\perp}$ and the ground-state
transition to the KD state, as shown in Fig. 1.

Conclusions and outlook: In summary, our study reveals a transition from
antiferromagnetic to Kekulé distorted states in graphene as a function of the
fine-structure constant $\kappa$ occurring around $\kappa\sim 0.5$. This
transition is driven by the enhancement of Landau-level mixing in the double-
exchange Feynman diagram. This diagram uniquely produces an attractive
$u_{\perp}$ from repulsive valley-exchange scattering $g_{\perp\perp}$ due to
matrix-element effects. The logarithmic divergence of $u_{z}$ and $u_{\perp}$
agrees qualitatively with prior renormalization group analyses conducted
without a magnetic field Kharitonov (2012); Aleiner  _et al._ (2007); sup .
This correspondence suggests that the notable enhancement of lattice-scale
interactions is largely attributed to the influence of remote Landau levels
sup . Our research advances the understanding of the phase diagram of charge
neutral graphene under the strong magnetic field. The calculated valley-
dependent Haldane pseudopotentials provide a foundation for accurate numerical
studies of SU(4) symmetry breaking in graphene’s fractional quantum Hall
states Peterson and Nayak (2013). Furthermore, our study sheds light on the
complex interplay between dielectric screening, valley- and sublattice-
dependent interactions, and unveils a novel mechanism for generating
attractive valley-exchange interactions from bare lattice-scale electron-
electron repulsion. This mechanism could be relevant in other contexts, such
as superconductivity in strongly-correlated multilayer graphene systems
Chatterjee  _et al._ (2022); Qin  _et al._ (2023); Huang  _et al._ (2022).

Acknowledgement We thank Jeanie Lau, J.I.A Li, X. Liu, A. H. MacDonald, G.
Murthy, A. Yazdani, A. Young, and H. Zhou for discussions. N. W. acknowledges
the hospitality of University of Kentucky where part of this work has been
performed. We acknowledge support from the Deutsche Forschungsgemeinschaft
(DFG) through research grant project number 518372354.

## References

  * Halperin and Jain (2020) B. I. Halperin and J. K. Jain,  _Fractional Quantum Hall Effects: New Developments_ (World Scientific, 2020).
  * Alicea and Fisher (2006) J. Alicea and M. P. A. Fisher, Phys. Rev. B 74, 075422 (2006).
  * Nomura and MacDonald (2006) K. Nomura and A. H. MacDonald, Phys. Rev. Lett. 96, 256602 (2006).
  * Yang  _et al._ (2006) K. Yang, S. Das Sarma, and A. H. MacDonald, Phys. Rev. B 74, 075423 (2006).
  * Checkelsky  _et al._ (2009) J. G. Checkelsky, L. Li, and N. P. Ong, Physical Review B 79, 115434 (2009).
  * Nomura  _et al._ (2009) K. Nomura, S. Ryu, and D.-H. Lee, Phys. Rev. Lett. 103, 216801 (2009).
  * Dean  _et al._ (2011) C. R. Dean, A. F. Young, P. Cadden-Zimansky, L. Wang, H. Ren, K. Watanabe, T. Taniguchi, P. Kim, J. Hone, and K. Shepard, Nature Physics 7, 693 (2011).
  * Goerbig (2011) M. O. Goerbig, Rev. Mod. Phys. 83, 1193 (2011).
  * Sodemann and MacDonald (2014) I. Sodemann and A. H. MacDonald, Phys. Rev. Lett. 112, 126804 (2014).
  * Spanton  _et al._ (2018) E. M. Spanton, A. A. Zibrov, H. Zhou, T. Taniguchi, K. Watanabe, M. P. Zaletel, and A. F. Young, Science 360, 62 (2018).
  * Zibrov  _et al._ (2018) A. Zibrov, E. Spanton, H. Zhou, C. Kometter, T. Taniguchi, K. Watanabe, and A. Young, Nature Physics 14, 930 (2018).
  * Zhou  _et al._ (2020) H. Zhou, H. Polshyn, T. Taniguchi, K. Watanabe, and A. Young, Nature Physics 16, 154 (2020).
  * Veyrat  _et al._ (2020) L. Veyrat, C. Déprez, A. Coissard, X. Li, F. Gay, K. Watanabe, T. Taniguchi, Z. Han, B. A. Piot, H. Sellier, _et al._ , Science 367, 781 (2020).
  * Farahi  _et al._ (2023) G. Farahi, C.-L. Chiu, X. Liu, Z. Papic, K. Watanabe, T. Taniguchi, M. P. Zaletel, and A. Yazdani, “Broken symmetries and excitation spectra of interacting electrons in partially filled landau levels,” (2023), arXiv:2303.16993 [cond-mat.mes-hall] .
  * Lu  _et al._ (2023) Z. Lu, T. Han, Y. Yao, A. P. Reddy, J. Yang, J. Seo, K. Watanabe, T. Taniguchi, L. Fu, and L. Ju, arXiv preprint arXiv:2309.17436 (2023).
  * Sharpe  _et al._ (2019) A. L. Sharpe, E. J. Fox, A. W. Barnard, J. Finney, K. Watanabe, T. Taniguchi, M. Kastner, and D. Goldhaber-Gordon, Science 365, 605 (2019).
  * Serlin  _et al._ (2020) M. Serlin, C. Tschirhart, H. Polshyn, Y. Zhang, J. Zhu, K. Watanabe, T. Taniguchi, L. Balents, and A. Young, Science 367, 900 (2020).
  * Li  _et al._ (2019) S.-Y. Li, Y. Zhang, L.-J. Yin, and L. He, Phys. Rev. B 100, 085437 (2019).
  * Liu  _et al._ (2022) X. Liu, G. Farahi, C.-L. Chiu, Z. Papic, K. Watanabe, T. Taniguchi, M. P. Zaletel, and A. Yazdani, Science 375, 321 (2022).
  * Coissard  _et al._ (2022) A. Coissard, D. Wander, H. Vignaud, A. G. Grushin, C. Repellin, K. Watanabe, T. Taniguchi, F. Gay, C. B. Winkelmann, H. Courtois, _et al._ , Nature 605, 51 (2022).
  * Zhou  _et al._ (2022) H. Zhou, C. Huang, N. Wei, T. Taniguchi, K. Watanabe, M. P. Zaletel, Z. Papić, A. H. MacDonald, and A. F. Young, Phys. Rev. X 12, 021060 (2022).
  * Wei  _et al._ (2018) D. S. Wei, T. van der Sar, S. H. Lee, K. Watanabe, T. Taniguchi, B. I. Halperin, and A. Yacoby, Science 362, 229 (2018).
  * Stepanov  _et al._ (2018) P. Stepanov, S. Che, D. Shcherbakov, J. Yang, R. Chen, K. Thilahar, G. Voigt, M. W. Bockrath, D. Smirnov, K. Watanabe, _et al._ , Nature Physics 14, 907 (2018).
  * Young  _et al._ (2014) A. Young, J. Sanchez-Yamagishi, B. Hunt, S. Choi, K. Watanabe, T. Taniguchi, R. Ashoori, and P. Jarillo-Herrero, Nature 505, 528 (2014).
  * Fu  _et al._ (2021) H. Fu, K. Huang, K. Watanabe, T. Taniguchi, and J. Zhu, Phys. Rev. X 11, 021012 (2021).
  * Wei  _et al._ (2021) N. Wei, C. Huang, and A. H. MacDonald, Phys. Rev. Lett. 126, 117203 (2021).
  * Kharitonov (2012) M. Kharitonov, Physical Review B 85, 155439 (2012).
  * Peterson and Nayak (2014) M. R. Peterson and C. Nayak, Phys. Rev. Lett. 113, 086401 (2014).
  * Sodemann and MacDonald (2013) I. Sodemann and A. H. MacDonald, Phys. Rev. B 87, 245425 (2013).
  * Lemonik  _et al._ (2012) Y. Lemonik, I. Aleiner, and V. Fal’Ko, Physical review b 85, 245451 (2012).
  * Aleiner  _et al._ (2007) I. Aleiner, D. Kharzeev, and A. Tsvelik, Physical Review B 76, 195415 (2007).
  * Knothe and Fal’ko (2020) A. Knothe and V. Fal’ko, Physical Review B 101, 235423 (2020).
  * (33) See Supplemental Material for the estimates of bare coupling and discussions on the perturbative calculations.
  * Raines  _et al._ (2021) Z. M. Raines, V. I. Fal’ko, and L. I. Glazman, Phys. Rev. B 103, 075422 (2021).
  * Haldane (1983) F. D. M. Haldane, Phys. Rev. Lett. 51, 605 (1983).
  * Peterson and Nayak (2013) M. R. Peterson and C. Nayak, Phys. Rev. B 87, 245129 (2013).
  * Haldane (1987) F. D. M. Haldane, in  _The Quantum Hall Effect_ , Vol. 1, edited by R. Prange and S. Girvin (Springer-Verlag New York, 1987) Chap. 8.
  * Berk and Schrieffer (1966) N. F. Berk and J. R. Schrieffer, Phys. Rev. Lett. 17, 433 (1966).
  * Note (1) While our conclusion is based on double-exchange diagram, the Berk-Schrieffer diagram consists of an infinite series of repeated exchange diagrams.
  * Das  _et al._ (2022) A. Das, R. K. Kaul, and G. Murthy, Phys. Rev. Lett. 128, 106803 (2022).
  * De  _et al._ (2023) S. J. De, A. Das, S. Rao, R. K. Kaul, and G. Murthy, Physical Review B 107, 125422 (2023).
  * Stefanidis and Villadiego (2023a) N. Stefanidis and I. S. Villadiego, Phys. Rev. B 108, 235137 (2023a).
  * Stefanidis and Villadiego (2023b) N. Stefanidis and I. S. Villadiego, Phys. Rev. B 107, 045132 (2023b).
  * Chatterjee  _et al._ (2022) S. Chatterjee, T. Wang, E. Berg, and M. P. Zaletel, Nature communications 13, 6013 (2022).
  * Qin  _et al._ (2023) W. Qin, C. Huang, T. Wolf, N. Wei, I. Blinov, and A. H. MacDonald, Physical Review Letters 130, 146001 (2023).
  * Huang  _et al._ (2022) C. Huang, N. Wei, W. Qin, and A. H. MacDonald, Physical Review Letters 129, 187001 (2022).
  * Herring (1966) C. Herring,  _Magnetism: Exchange interactions among itinerant electrons_ (Academic Press, 1966).
  * Wu  _et al._ (2018) F. Wu, A. H. MacDonald, and I. Martin, Physical review letters 121, 257001 (2018).
  * Clementi and Raimondi (1963) E. Clementi and D.-L. Raimondi, The Journal of Chemical Physics 38, 2686 (1963).
  * Note (2) We take the convention $\sigma^{\pm}=(\sigma^{x}\pm i\sigma^{y})/2$.
  * Gradshteyn and Ryzhik (2014) I. S. Gradshteyn and I. M. Ryzhik,  _Table of integrals, series, and products_ (Academic press, 2014).

Supplementary Material: Landau-Level Mixing and SU(4) Symmetry Breaking in
Graphene

##  I Estimation of the bare g-values from tight-binding calculation

To compute the contribution of electron-electron repulsion to
$g_{\alpha\mu}^{0}$, we adopt a basic tight-binding model, assuming that the
Bloch function of graphene is a linear combination of the $p_{z}$ atomic
orbitals. For this calculation, we only require two Bloch functions situated
at the non-equivalent corners, $\bm{\tau}=\pm\bm{K}$, of the Brillouin zone
and their two sublattice projections. They can be grouped into a vector,

$v=\left[u_{KA},u_{KB},u_{K^{\prime}B},-u_{K^{\prime}A}\right]^{T},$ (S1)
$\langle\bm{r}\ket{u_{\tau\sigma}}=\sqrt{\frac{A_{0}}{N}}\sum_{\bm{R}}e^{i\bm{\tau}\cdot(\bm{R}_{\sigma}+\bm{R})}\phi(\bm{r}-\bm{R}_{\sigma}-\bm{R}).$
(S2)

Here $\bm{R}_{\sigma}$ represents the sublattice vector so that $\bm{R}_{A}=0$
and $\bm{R}_{B}=\bm{R}_{AB}$. In our coordinate system, the origin coincides
with sublattice $A$ and and $\bm{R}_{AB}$ is used to denote the relative
displacement between the two carbon sublattices inside the unit cell.
Consequently, the summation over $\bm{R}_{A}$ covers the set:
$\\{\bm{R}_{A}|m\bm{a}_{1}+n\bm{a}_{2}\,\forall m,n\in\mathbb{Z}$} and the
summation over $\bm{R}_{B}$ covers the set:
$\\{\bm{R}_{B}|\bm{R}_{AB}+m\bm{a}_{1}+n\bm{a}_{2}\,\forall m,n\in\mathbb{Z}$}
where $\bm{a}_{1}$ and $\bm{a}_{2}$ are the primitive lattice vectors.

The normalized wavefunction for $p_{z}$ orbitals is given by

$\phi(r,\theta,\phi)=\frac{2\tilde{Z}^{2}\sqrt{\tilde{Z}}}{\sqrt{96a_{B}^{5}}}re^{-\frac{\tilde{Z}r}{2a_{B}}}Y_{1}^{0}(\theta,\phi),$
(S3)

where $\int|\phi(\bm{r})|^{2}d^{3}r=1$. Here, $a_{B}$ is the Bohr radius and
$\tilde{Z}$ is the effective nuclear charge of carbon atom. With this basis in
place, $g_{\alpha\mu}^{0}$ can be directly computed as,

$\displaystyle g_{\alpha\mu}^{0}$
$\displaystyle=\frac{1}{4A_{0}}\sum_{ijkl}V_{ik,jl}(\tau^{\alpha}\sigma^{\mu})_{ij}(\tau^{\alpha}\sigma^{\mu})_{kl},$
(S4)

where the Coulomb matrix elements are described by the six-dimensional
integrals,

$V_{ik,jl}=e^{2}\int\int_{u.c.}\frac{v^{*}_{i}(\bm{r}_{1})v_{j}(\bm{r}_{1})v^{*}_{k}(\bm{r}_{2})v_{l}(\bm{r}_{2})}{|\bm{r}_{1}-\bm{r}_{2}|}\;d^{3}r_{1}d^{3}r_{2}.$
(S5)

Here $A_{0}$ is the area of the graphene unit cell. Next, we express Eq. (S4)
as the Coulomb potential energy between overlap charges:

$\displaystyle g_{zz}$
$\displaystyle=\frac{e^{2}}{4A_{0}}\sum_{\tau,\sigma,\tau^{\prime},\sigma^{\prime}}sgn(\sigma)sgn(\sigma^{\prime})\int\int_{u.c.}d^{3}\bm{r}_{1}\int
d^{3}\bm{r}_{2}\frac{\rho_{\tau\sigma,\tau\sigma}(\bm{r}_{1})\rho_{\tau^{\prime}\sigma^{\prime},\tau^{\prime}\sigma^{\prime}}(\bm{r}_{2})}{|\bm{r}_{1}-\bm{r}_{2}|}$
(S6) $\displaystyle g_{\perp\perp}$
$\displaystyle=\frac{e^{2}}{4A_{0}}\sum_{\tau\sigma}\int\int_{u.c.}d^{3}\bm{r}_{1}\int
d^{3}\bm{r}_{2}\frac{\rho_{\tau\sigma,\bar{\tau}\sigma}(\bm{r}_{1})\rho_{\bar{\tau}\sigma,\tau\sigma}(\bm{r}_{2})}{|\bm{r}_{1}-\bm{r}_{2}|}$
(S7) $\displaystyle g_{z\perp}$
$\displaystyle=\frac{e^{2}}{4A_{0}}\sum_{\tau\sigma}2\int\int_{u.c.}d^{3}\bm{r}_{1}\int
d^{3}\bm{r}_{2}\frac{\rho_{\tau\bar{\sigma},\tau\sigma}(\bm{r}_{1})\rho_{\tau\sigma,\tau\bar{\sigma}}(\bm{r}_{2})}{|\bm{r}_{1}-\bm{r}_{2}|}$
(S8) $\displaystyle g_{\perp z}$
$\displaystyle=\frac{e^{2}}{4A_{0}}\sum_{\tau\sigma}2\int\int_{u.c.}d^{3}\bm{r}_{1}\int
d^{3}\bm{r}_{2}\frac{\rho_{\tau\sigma,\bar{\tau}\bar{\sigma}}(\bm{r}_{1})\rho_{\bar{\tau}\bar{\sigma},\tau\sigma}(\bm{r}_{2})}{|\bm{r}_{1}-\bm{r}_{2}|}.$
(S9)

We define the overlap charge density $\rho_{ij}(\bm{r})$ as
$v^{*}_{i}(\bm{r})v_{j}(\bm{r})$, where $i,j$ are indices representing
$KA,KB,K^{\prime}A,K^{\prime}B$. Note that the overlap charge within the same
sublattice, termed as the intra-sublattice overlap charge, is significantly
larger than the overlap charge between different sublattices, termed as the
inter-sublattice overlap charge. Consequently, we anticipate that the
magnitudes of $g^{0}_{zz}$ and $g^{0}_{\perp\perp}$ will be greater compared
to those of $g^{0}_{z\perp}$ and $g^{0}_{\perp z}$. When considering the
intra-sublattice overlap charge, this term can be approximated to leading
order by the separation distance between two $p_{z}$ atomic orbitals:

$\displaystyle\rho_{\tau_{i}\sigma,\tau_{j}\sigma}(\bm{r})\equiv
v^{*}_{\tau_{i}\sigma}(\bm{r})v_{\tau_{j}\sigma}(\bm{r})$
$\displaystyle=\frac{A_{0}}{N}\sum_{\bm{R}_{1},\bm{R}_{2}}e^{-i\bm{\tau}_{i}\cdot(\bm{R}_{1}+\bm{R}_{\sigma})+i\bm{\tau}_{j}\cdot(\bm{R}_{2}+\bm{R}_{\sigma})}\phi(\bm{r}-\bm{R}_{1}-\bm{R}_{\sigma})\phi(\bm{r}-\bm{R}_{2}-\bm{R}_{\sigma})$
(S10) $\displaystyle\approx
A_{0}\sum_{\bm{R}}e^{-i(\bm{\tau}_{i}-\bm{\tau}_{j})\cdot(\bm{R}+\bm{R}_{\sigma})}\phi^{2}(\bm{r}-\bm{R}-\bm{R}_{\sigma})$

The second line above is equivalent to selecting the $N$ combinations within
the summations where $\bm{R}_{1}$ is equal to $\bm{R}_{2}$. We have
numerically confirmed that contributions from subsequent $R$ shells are
notably smaller, showing a reduction of at least an order of magnitude. This
significant decrease is mainly due to the compact nature of the $p_{z}$
orbitals in comparison to the distance between carbon atoms.

For inter-sublattice overlap charge, the two $p_{z}$ orbitals are separated by
the carbon-carbon distance $\bm{R}_{AB}$. Consequently, the leading
approximation is given by the following:

$\displaystyle\rho_{\tau_{i}\bar{\sigma},\tau_{j}\sigma}(\bm{r})\equiv
v^{*}_{\tau_{i}\bar{\sigma}}(\bm{r})v_{\tau_{j}\sigma}(\bm{r})$
$\displaystyle=\frac{A_{0}}{N}\sum_{\bm{R}_{1},\bm{R}_{2}}e^{-i\bm{\tau}_{i}\cdot(\bm{R}_{1}+\bm{R}_{\bar{\sigma}})+i\bm{\tau}_{j}\cdot(\bm{R}_{2}+\bm{R}_{\sigma})}\phi(\bm{r}-\bm{R}_{1}-\bm{R}_{\bar{\sigma}})\phi(\bm{r}-\bm{R}_{2}-\bm{R}_{\sigma})$
(S11) $\displaystyle\approx
A_{0}\sum_{\bm{R}}e^{-i(\bm{\tau}_{i}-\bm{\tau}_{j})\cdot\bm{R}}e^{-i\tau_{i}\cdot\bm{R}_{\bar{\sigma}}+i\tau_{j}\cdot\bm{R}_{\sigma}}\left(\phi(\bm{r}-\bm{R}-\bm{R}_{\bar{\sigma}})\phi(\bm{r}-\bm{R}-\bm{R}_{\sigma})\right.$
$\displaystyle+\phi(\bm{r}-\bm{R}-\bm{a}_{1}-\bm{R}_{\bar{\sigma}})\phi(\bm{r}-\bm{R}-\bm{R}_{\sigma})e^{-i\bm{\tau}_{i}\cdot\bm{a}_{1}}+\phi(\bm{r}-\bm{R}-\bm{a}_{2}-\bm{R}_{\bar{\sigma}})\phi(\bm{r}-\bm{R}-\bm{R}_{\sigma})e^{-i\bm{\tau}_{i}\cdot\bm{a}_{2}}\bigg{)}$

We have also numerically verified that the contribution from the next order in
the $R$-shell is smaller by an order of magnitude. After substituting the
expressions for intra- and intersublattice overlap charges into Eq.(S6), we
have two independent summations over lattice vectors, resulting from the
squaring of the overlap charge. We then apply a change of variables to the
center-of-mass coordinates, denoted as $\bm{R}_{cm}$, and relative
coordinates, $\bm{R}$. The summation over $\bm{R}_{cm}$ efficiently replaced
the integration over the unit cell with an integration across the entire
spatial domain. Consequently, we are primarily left with the task of dealing
with the summation over $\bm{R}$:

$\displaystyle g_{zz}=$ $\displaystyle
8\sum_{\bm{R}}\left[\frac{e^{2}A_{0}}{4}\int\frac{d^{3}\bm{q}}{(2\pi)^{3}}\frac{4\pi|\rho^{a}(\bm{q})|^{2}}{\bm{q}^{2}}(1-e^{i\bm{q}\cdot\bm{R}_{AB}})e^{i\bm{q}\cdot\bm{R}}\right],$
(S12) $\displaystyle g_{\perp\perp}=$ $\displaystyle
4\sum_{\bm{R}}\left[\frac{e^{2}A_{0}}{4}\int\frac{d^{3}\bm{q}}{(2\pi)^{3}}\frac{4\pi|\rho^{a}(\bm{q})|^{2}}{\bm{q}^{2}}e^{i\bm{q}\cdot\bm{R}}e^{i\bm{R}\cdot\Delta\bm{K}}\right],$
(S13) $\displaystyle g_{z\perp}=$ $\displaystyle
8\sum_{\bm{R}}\left[\frac{e^{2}A_{0}}{4}\int\frac{d^{3}\bm{q}}{(2\pi)^{3}}\frac{4\pi|\rho^{e}(\bm{q})|^{2}}{\bm{q}^{2}}e^{i\bm{q}\cdot\bm{R}}\right],$
(S14) $\displaystyle g_{\perp z}=$ $\displaystyle
8\sum_{\bm{R}}\left[\frac{e^{2}A_{0}}{4}\int\frac{d^{3}\bm{q}}{(2\pi)^{3}}\frac{4\pi|\rho^{e}(\bm{q})|^{2}}{\bm{q}^{2}}e^{i\bm{q}\cdot\bm{R}}\right],$
(S15) $\displaystyle\rho^{a}(\bm{q})$ $\displaystyle=\int
d^{3}\bm{r}\phi(\bm{r})\phi(\bm{r})e^{-i\bm{q}\cdot\bm{r}}$ (S16)
$\displaystyle\rho^{e}(\bm{q})$ $\displaystyle=\int
d^{3}\bm{r}\bigg{(}\phi(\bm{r})\phi(\bm{r}-\bm{R}_{ab})+e^{i\bm{K}\cdot\bm{a}_{2}}\phi(\bm{r})\phi(\bm{r}-\bm{a}_{2}-\bm{R}_{ab})+e^{i\bm{K}\cdot\bm{a}_{1}}\phi(\bm{r})\phi(\bm{r}-\bm{a}_{1}-\bm{R}_{ab})\bigg{)}e^{-i\bm{q}\cdot\bm{r}}.$
(S17)

Figure S1: The intrasublattice overlap charge $\rho^{a}(\bm{G},q_{z})$ versus
$q_{z}$. As the effective nuclear charge $\tilde{Z}$ increases, leading to a
more localized wave-function, the Fourier components decay slower. Figure S2:
The intersublattice overlap charge $\rho^{e}(\bm{G}=\bm{G_{1s}},q_{z})$ versus
$q_{z}$ for $\bm{G}$. Note the imaginary part is an order of magnitude
stronger than the real part.

Note the terms inside the square brackets $[...]$ decay very slowly as
$\bm{R}$ increases. For example,

$\int\frac{d^{3}\bm{q}}{(2\pi)^{3}}\frac{|\rho^{a}(\bm{q})|^{2}}{\bm{q}^{2}}e^{i\bm{q}\cdot\bm{R}}=\int
d^{3}\bm{r}_{1}d^{3}\bm{r}_{2}\frac{\phi^{2}(\bm{r}_{1})\phi^{2}(\bm{r}_{2})}{|\bm{r}_{1}-\bm{r}_{2}+\bm{R}|}\xrightarrow[]{\bm{R}\rightarrow\infty}{}\frac{1}{|\bm{R}|}$
(S18)

Therefore, given this poor convergence, it becomes more pragmatic to evaluate
them in the Fourier space. Using
$\sum_{\bm{R}}e^{i\bm{q}\cdot\bm{R}}=\frac{(2\pi)^{2}}{A_{0}}\sum_{\bm{G}}\delta^{(2)}(\bm{q}-\bm{G})$
where $\bm{G}$ is the 2D reciprocal lattice vectors, we arrive at the
following expressions:

$\displaystyle
g^{0}_{\perp\perp}=4\sum_{\bm{G}}\frac{e^{2}}{4}\int\frac{dq_{z}}{2\pi}\frac{4\pi|\rho^{a}(\bm{G}-\Delta\bm{K},q_{z})|^{2}}{(\bm{G}-\Delta\bm{K})^{2}+q_{z}^{2}}$
(S19) $\displaystyle
g^{0}_{zz}=8\sum_{\bm{G}}\frac{e^{2}}{4}\int\frac{dq_{z}}{2\pi}\frac{4\pi|\rho^{a}(\bm{G},q_{z})|^{2}}{\bm{G}^{2}+q_{z}^{2}}(1-\cos(\bm{G}\cdot\bm{R}_{ab}))$
(S20) $\displaystyle
g_{z\perp}^{0}=8\sum_{\bm{G}}\frac{e^{2}}{4}\int\frac{dq_{z}}{2\pi}\frac{4\pi|\rho^{e}(\bm{G},q_{z})|^{2}}{\bm{G}^{2}+q_{z}^{2}}$
(S21) $\displaystyle g_{\perp
z}^{0}=8\sum_{\bm{G}}\frac{e^{2}}{4}\int\frac{dq_{z}}{2\pi}\frac{4\pi|\rho^{e}(\bm{G}-\Delta\bm{K},q_{z})|^{2}}{(\bm{G}-\Delta\bm{K})^{2}+q_{z}^{2}}$
(S22)

Here $\Delta\bm{K}=\bm{K}-\bm{K}^{\prime}$ is the momentum transfer. From the
above equations, it is evident that they are all positive definite.
Specifically, the positivity of $g_{\perp\perp}^{0}$ and $g_{z\perp}^{0}$
originates from the properties of the exchange integral, as discussed Ref.
Herring (1966). We note that the estimates presented in Ref. Knothe and Fal’ko
(2020) seem to be inaccurate, especially regarding their result indicating
$g_{\perp\perp}^{0}<0$. This discrepancy might stem from their approach of
performing the summation in real space instead of in reciprocal space. Fig. S2
and Fig. S3 show the Fourier components of the overlap charge vs $q_{z}$ for
different effective nuclear charge $\tilde{Z}$. Note that $\rho^{a}$ is real
while $\rho^{e}$ is complex.

Next, we introduce truncation in $\bm{G}$. For $g^{0}_{\perp\perp}$ and
$g^{0}_{\perp z}$, which involve inter-valley scattering with a momentum
transfer of $\Delta\bm{K}$, the primary contributions to the reciprocal
lattice sum arise mainly from the $\bm{G}=0$ and the two corners in the first
shell closest to $\Delta\bm{K}$, as indicated by the dashed lines in Fig.
(S3). In contrast, for $g^{0}_{zz}$ and $g^{0}_{z\perp}$, corresponding to
zero-momentum forward scattering, the major contributions arise from the six
corners in the first shell, defined as $\bm{G}_{1s}$ in Fig. (S3). For
$g^{0}_{zz}$, the absence of contribution from $\bm{G}=0$ can be attribute to
the $1-\cos(\bm{G}\cdot\bm{R}_{ab})$ factor, a result stemming from AA-
scattering minus AB-scattering. For $g^{0}_{z\perp}$ the absence of
contribution from $\bm{G}=0$ is consequence of $C_{3}$ rotation symmetry,
which implies $\rho^{e}(0)=0$. Numerical verification has confirmed that
contributions from the next shell are significantly smaller, by an order of
magnitude, compared to those from the first shell. Consequently, evaluating
the overlap charge using the leading contributions in $\bm{R}^{\prime}$ and
computing $g^{0}_{\alpha\mu}$ to leading contributions in $\bm{G}$ yield the
results summarized in the table below:

Figure S3: The Brillouin zone of monolayer graphene. Here
$\bm{b}_{1}=\frac{2\pi}{3a}(1,\sqrt{3})$ and
$\bm{b}_{2}=\frac{2\pi}{3a}(1,-\sqrt{3})$, $a\approx 1.42\text{\AA}$ is
carbon-carbon distance. The $\Gamma$ point corresponds to $\bm{G}=0$ and the
first shell corresponds to $\bm{G}_{1s}=\\{G|\text{for
$\bm{G}$}\in(\pm\bm{b}_{1},\pm\bm{b}_{2},\pm(\bm{b}_{1}+\bm{b}_{2}))$}. The
two corners in the first shell, which are connected by a dashed line, are
equidistant from the zero shell at $\Delta\bm{K}$.

$meVnm^{2}$ | $\tilde{Z}=1$ | $\tilde{Z}=1.5$ | $\tilde{Z}=2$ | $\tilde{Z}=2.5$ | $\tilde{Z}=3$ | $\tilde{Z}=4$
---|---|---|---|---|---|---
$g_{\perp\perp}^{0}$ | 8.88 | 61.5 | 156.4 | 268.9 | 384.3 | 600.6
$g_{zz}^{0}$ | 0.55 | 12.55 | 67.46 | 184.03 | 353.7 | 781.15
$g_{z\perp}^{0}$ | 0.234 | 4.41 | 15.84 | 25.11 | 25.38 | 12.42
$g_{\perp z}^{0}$ | 5.94 | 29.7 | 46.08 | 42.93 | 30.42 | 9.72

Note that the intersublattice ( $g_{\perp z}^{0}$ and $g_{z\perp}^{0}$)
contribution is almost an order of magnitude smaller than the intra-sublattice
($g_{\perp\perp}^{0}$ and $g_{zz}^{0}$) contribution. However, the values for
$g_{\perp z}^{0}$ and $g_{z\perp}^{0}$ actually has a strong contribution
coming from electron-optical-phonon scattering. The estimates for these
interactions Wu  _et al._ (2018) are:

$(g_{\perp
z}^{0},\,g_{z\perp}^{0})_{phonon}=(-69,-52)\text{meV}\cdot\text{nm}^{2}$ (S23)

By setting the effective nucleus charge $\tilde{Z}=2.5$, a value slightly
reduced from the $\tilde{Z}=3$ typically associated with an isolated carbon
atom Clementi and Raimondi (1963), and factoring in the electron-phonon
contributions, we derived the following estimates:

$(g_{\perp\perp}^{0},g_{zz}^{0},g_{\perp
z}^{0},\,g_{z\perp}^{0})=(269,184,-26,-27)\text{meV}\cdot\text{nm}^{2}.$ (S24)

It is important to point out that our central result, namely that an increase
in $\kappa$ leads to a more pronounced decrease in $u_{\perp}$ relative to the
increase of $u_{z}$, remains consistent regardless of the choice of
$\tilde{Z}$. To leading order in electron-electron and electron-phonon
contribution, $g_{\perp 0}^{0}=g_{0\perp}^{0}=g_{z0}^{0}=g_{0z}^{0}=0$
Kharitonov (2012).

##  II Calculation Details of Landau level mixing

In the weakly interacting systems, the perturbative expansion of low-energy
effective Hamiltonian reads that
$\hat{H}=\hat{H}_{0}+\hat{V}^{(1)}+\hat{V}^{(2)}+...$, where $\hat{H}_{0}$ is
the non-interacting Hamiltonian, $\hat{V}^{(1)}$ is the projection of the bare
interactions onto the degenerate non-interactng ground state manifold of
energy $E_{0}$. In this study, the non-interacting ground states have two
electrons in the empty $n=0$ Landau level (LL) along with fully
empty(occupied) $n>0(n<0)$ LLs Sodemann and MacDonald (2013).The next order
correction to the effective Hamiltonian is obtained from the second order
perturbation theory,

$\displaystyle\hat{V}^{(2)}$
$\displaystyle=-\mathcal{P}(\hat{H}_{s}+\hat{H}_{a})\mathcal{P}_{\perp}\frac{1}{\hat{H}_{0}-E_{0}}\mathcal{P}_{\perp}(\hat{H}_{s}+\hat{H}_{a})\mathcal{P}$
$\displaystyle=-\mathcal{P}\hat{H}_{s}\mathcal{P}_{\perp}\frac{1}{\hat{H}_{0}-E_{0}}\mathcal{P}_{\perp}\hat{H}_{s}\mathcal{P}-\left(\mathcal{P}\hat{H}_{s}\mathcal{P}_{\perp}\frac{1}{\hat{H}_{0}-E_{0}}\mathcal{P}_{\perp}\hat{H}_{a}\mathcal{P}+h.c.\right)-\mathcal{P}\hat{H}_{a}\mathcal{P}_{\perp}\frac{1}{\hat{H}_{0}-E_{0}}\mathcal{P}_{\perp}\hat{H}_{a}\mathcal{P}$
$\displaystyle\equiv\hat{V}_{ss}^{(2)}+\hat{V}_{sa}^{(2)}+\hat{V}_{aa}^{(2)}.$
(S25)

The first term $\hat{V}_{ss}^{(2)}$ has been studied in previous works and
could potentially play a crucial role in certain fractional quantum Hall
states Peterson and Nayak (2013, 2014). However, due to its spin-valley SU(4)
symmetry, it cannot differentiate competing spin-valley orders and will not be
further discuss here. $\hat{V}_{sa}$ gives rise to first-order-in-$\kappa$
corrections to the valley-dependent effective interactions, while the
quadratic-in-$g_{\alpha\mu}^{0}$ corrections $\hat{V}_{aa}$ is independent of
$\kappa$ (i.e., dielectric screening environment).

As we briefly explained in the main text, $\hat{V}^{(2)}$ admits the
decomposition $\hat{V}^{(2)}=\hat{V}_{3b}+\hat{V}_{2b}+\hat{V}_{1b}+constant$
based on the number $n_{v}$ of the virtually excited single-particle states.
Note that $n_{v}\neq 0$ since the intermediate many-body states must have
energy different from $E_{0}$. The $n_{v}=1$ and $2$ terms give rise to three-
body interactions $\hat{V}_{3b}$ and two-body interactions $\hat{V}_{2b}$
Sodemann and MacDonald (2013), respectively. $n_{v}=3$ corresponds to a SU(4)
symmetric single-particle self-energy $\hat{V}_{1b}$ and can be omitted
together with the constant energy shift generated by the $n_{v}=4$ term.
$\hat{V}_{3b}=0$ in the zeroth LL of monolayer graphene thanks to the
particle-hole symmetric bare Hamiltonian and the particle-hole symmetry
between the $N>0$ and $N<0$ remote LLs Peterson and Nayak (2013). The
corrections to valley-dependent effective interactions are entirely
contributed by $\hat{V}_{2b}$.

To study the two-body interactions preserving the translational and rotational
symmetries, it is convenient to use center of mass and relative coordinates,

$A=\frac{a_{i}+a_{j}}{\sqrt{2}},\ B=\frac{b_{i}+b_{j}}{\sqrt{2}},\
a=\frac{a_{i}-a_{j}}{\sqrt{2}},\ b=\frac{b_{i}-b_{j}}{\sqrt{2}},$ (S26)

where $a_{i},b_{i}$ are the lowering operators in the LL and guiding center
space, respectively. We construct states with well defined center of mass and
relative coordinate quantum numbers,

$|NM\rangle_{c}|nm\rangle_{r}=\frac{A^{\dagger N}B^{\dagger M}a^{\dagger
n}b^{\dagger m}}{\sqrt{N!M!n!m!}}|0\rangle.$ (S27)

For a two-body interaction $\hat{V}$ preserving the translational and
rotational symmetries,

${}_{c}\bra{N^{\prime}M^{\prime}}{}_{r}\bra{n^{\prime}m^{\prime}}\hat{V}\ket{NM}_{c}\ket{nm}_{r}={}_{r}\bra{n^{\prime}m^{\prime}}\hat{V}\ket{nm}_{r}\delta_{m^{\prime}-n^{\prime},m-n}\delta_{N^{\prime},N}\delta_{M^{\prime},M},$
(S28)

where the subscript ‘c’ and ‘r’ stands for the center of mass and relative
coordinates, respectively. ${}_{r}\langle 0m|\hat{V}|0m\rangle_{r}$ are the
Haldane pseudopotentials in the zeroth LL. We follow Ref. Sodemann and
MacDonald (2013) to evaluate in the next two sections the corrections to the
valley-dependent Haldane pseudopotentials, $u_{\alpha}^{(2)}(m)/2$ induced by
$V_{sa}^{(2)}$ and $u_{\alpha}^{(2)^{\prime}}(m)/2$ induced by $V_{aa}^{(2)}$.

###  II.1 $V_{sa}-$induced Haldane pseudopotential corrections
$u_{\alpha}^{(2)}(m)$

$u_{\alpha}^{(2)}(m)$ receive contributions from three types of non-vanishing
Feynman diagrams depicted in Fig.2: (a) the particle-particle ladder diagram,
(b) the double-exchange diagram, and (c) the vertex correction to the valley-
dependent interactions. As discussed in the main text,

$u_{\alpha}^{(2)}(m)=\frac{\kappa}{2\pi
l_{B}^{2}}\sum_{\mu=0,x,y,z}(c_{m,a}^{\mu}+c_{m,b}^{\mu}+c_{m,c}^{\mu})g_{\alpha\mu}^{0}.$
(S29)

Because $g_{\alpha x}^{0}=g_{\alpha y}^{0}=g_{\alpha\perp}^{0}$, it will be
convenient to define $c_{m,a/b/c}^{\perp}=c_{m,a/b/c}^{x}+c_{m,a/b/c}^{y}$.

It is straightforward to derive $c_{m,a}^{0,z}$ under two-particle bases
($c_{m,a}^{\perp}=0$ due to the valley-sublattice locking in the zeroth LL,
see the main text.):

$\frac{c_{m,a}^{0,z}}{4\pi}=-2\sum_{n=1}^{2N_{c}}\prescript{}{r}{\bra{0m}}\delta(\bm{r}_{ij})\ket{n\
n+m}_{r}\prescript{}{r}{\bra{n\
n+m}}\frac{1}{r_{ij}}\ket{0m}_{r}\sum_{n_{i,j}=-\infty}^{\infty}\frac{\Theta(n_{i}n_{j})}{|\epsilon_{\overline{n_{i}}}+\epsilon_{\overline{n_{j}}}|}|\prescript{}{c}{\langle}0|\prescript{}{r}{\langle}n||n_{i}|,|n_{j}|\rangle|^{2}.$
(S30)

We measure the length in the unit of $l_{B}$ and the kinetic energy in the
unit of $\hbar v_{F}/l_{B}$ in this section. The interaction matrix elements
are derived as follows:

$\displaystyle{}_{r}\bra{n(m+n)}\frac{1}{r_{ij}}\ket{0m}_{r}$
$\displaystyle=\int\frac{d^{2}q}{(2\pi)^{2}}\frac{2\pi}{|\bm{q}|}\mathcal{F}_{n,0}(\sqrt{2}\bm{q})\mathcal{F}_{m+n,m}^{*}(-\sqrt{2}\bm{q})$
(S31)
$\displaystyle=\frac{(-1)^{n}}{2}\sqrt{\frac{(n+m)!}{n!m!}}\frac{\Gamma(n+\frac{1}{2})}{n!}{}_{2}F_{1}(-m,n+\frac{1}{2};n+1;1),$
(S32) $\displaystyle{}_{r}\bra{n(m+n)}\delta(\bm{r}_{ij})\ket{0m}_{r}$
$\displaystyle=\int\frac{d^{2}q}{(2\pi)^{2}}\mathcal{F}_{n,0}(\sqrt{2}\bm{q})\mathcal{F}_{m+n,m}^{*}(-\sqrt{2}\bm{q})$
$\displaystyle=\frac{(-1)^{n}}{4\pi}\sqrt{\frac{(n+m)!}{n!m!}}{}_{2}F_{1}(-m,n+1;n+1;1)=\frac{(-1)^{n}}{4\pi}\delta_{m,0}.$
(S33)

In Eq. (S31), we used the density form factor
$\mathcal{F}_{n^{\prime},n}(\bm{q})=\bra{n^{\prime}}e^{i\bm{q}\cdot(\hat{z}\crossproduct\bm{\Pi})}\ket{n}$
of the LLs of the two dimensional electron gas and an analog relation for the
guiding center coordinates
$\bra{m^{\prime}}e^{i\bm{q}\cdot\bm{R}}\ket{m}=\mathcal{F}_{m^{\prime},m}^{*}(-\bm{q})$.
$\mathcal{F}_{n^{\prime},n}(\bm{q})=\mathcal{F}_{n,n^{\prime}}^{*}(-\bm{q})$
and for $n^{\prime}\geq n$,

$\mathcal{F}_{n^{\prime},n}(\bm{q})=\sqrt{\frac{n!}{n^{\prime}!}}\left(\frac{q_{x}+iq_{y}}{\sqrt{2}}\right)^{n^{\prime}-n}L_{n}^{n^{\prime}-n}\left(\frac{|q|^{2}}{2}\right)e^{-|q|^{2}/4},$
(S34)

with $L_{n}^{m}(x)$ the generalized Laguerre polynomial. The $\sqrt{2}$ factor
in Eq. (S31) is because we are using the relative coordinates Eq. (S26).
${}_{2}F_{1}$ is the hypergeometric function. Eq. (II.1) indicates that a pair
of electrons with nonzero relative angular momentum do not interact via
contact interaction because their wave functions vanish at $\bm{r}_{ij}=0$.

The last term in Eq. (S30) is obtained from the relation

$\left|n_{1},n_{2}\right\rangle=\sum_{n=0}^{n_{1}+n_{2}}R_{n_{2},n}^{n_{1}+n_{2}}\times\left|n_{1}+n_{2}-n\right\rangle_{c}|n\rangle_{r},\quad(n_{1,2}\geq
0)$ (S35)

where $R^{L}$ is a $L\times L$ real, symmetric and orthogonal matrix whose
explicit form reads that

$R_{n^{\prime},n}^{L}=\sqrt{\frac{\binom{L}{n^{\prime}}}{2^{L}\binom{L}{n}}}\sum_{\nu=\text{max}(0,n+n^{\prime}-L)}^{\text{min}(n,n^{\prime})}\binom{L-n^{\prime}}{n-\nu}\binom{n^{\prime}}{\nu}(-1)^{\nu}.$
(S36)

Only $R_{n^{\prime},n}^{n}=(-1)^{n}\sqrt{\binom{n}{n^{\prime}}/2^{n}}$ is
needed in Eq. (S30) because the center of mass kinetic quantum number
$n_{1}+n_{2}-n=0$. Combining Eq. (S30)-(S36), we obtain $c_{m,a}^{0,z}=0$ for
$m\neq 0$ and as shown in Eq.(9) in the main text,

$\frac{c_{m=0,a}^{0,z}}{4\pi}=\sum_{n=1}^{2N_{c}}\frac{\Gamma(n+\frac{1}{2})}{2^{n+1}n!}\sum_{n_{1}=0}^{n}\frac{\binom{n}{n_{1}}}{\sqrt{2n_{1}}+\sqrt{2(n-n_{1})}}.$
(S37)

Here the cut off of the LL numbers is $2N_{c}$ because $n=|n_{i}|+|n_{j}|$ and
$|n_{i,j}|\leq N_{c}$.

The expressions for $c_{m,b/c}^{\mu}$ are given by Eqs.(10)-(11) in the main
text. In the interaction matrix elements $\mathcal{V}$ and
$\mathcal{V}^{\mu}$, we introduced the graphene form factors
$F_{\xi^{\prime}n^{\prime},\xi
n}^{\mu}(\bm{q})=\bra{\overline{\xi^{\prime}n^{\prime}}}\sigma^{\mu}e^{i\bm{q}\cdot(\hat{z}\times\bm{\Pi})}\ket{\overline{\xi
n}}$ for $\mu=0,\pm,z$ 222We take the convention $\sigma^{\pm}=(\sigma^{x}\pm
i\sigma^{y})/2$,

$\displaystyle F_{\xi^{\prime}n^{\prime},\xi
n}^{0/z}(\bm{q})=\frac{1}{2^{1-(\delta_{\xi,0}+\delta_{\xi^{\prime},0})/2}}\left(\mathcal{F}_{n^{\prime},n}(\bm{q})\pm\xi\xi^{\prime}\mathcal{F}_{n^{\prime}-1,n-1}(\bm{q})\right),$
(S38) $\displaystyle F_{\xi^{\prime}n^{\prime},\xi
n}^{+}(\bm{q})=\frac{1}{2^{1-(\delta_{\xi,0}+\delta_{\xi^{\prime},0})/2}}\xi\mathcal{F}_{n^{\prime},n-1}(\bm{q}),$
(S39) $\displaystyle F_{\xi^{\prime}n^{\prime},\xi
n}^{-}(\bm{q})=\frac{1}{2^{1-(\delta_{\xi,0}+\delta_{\xi^{\prime},0})/2}}\xi^{\prime}\mathcal{F}_{n^{\prime}-1,n}(\bm{q}),$
(S40)

where $\xi,\xi^{\prime}=\pm 1$ and $n,n^{\prime}\geq 0$. Here we require that
the density form factor $\mathcal{F}_{n^{\prime},n}(\bm{q})=0$ if
$n^{\prime}<0$ or $n<0$. Using the relation
$\prescript{}{r}{\bra{m}}e^{i\bm{q}\cdot\bm{\eta}_{ij}}\ket{m}_{r}=e^{-\frac{\bm{q}^{2}}{2}}L_{m}(\bm{q}^{2})$,
we find that

$\displaystyle\frac{c_{m,b}^{z}}{4\pi}$
$\displaystyle=\iint\frac{d^{2}\bm{q}_{1}d^{2}\bm{q}_{2}}{(2\pi)^{4}}\frac{2\pi}{|\bm{q}_{2}|}e^{-\bm{q}_{1}^{2}-\bm{q}_{2}^{2}+\bm{q}_{1}\cdot\bm{q}_{2}}L_{m}(|\bm{q}_{1}-\bm{q}_{2}|^{2})\sideset{}{{}^{\prime}}{\sum}_{n,n^{\prime}=0}^{N_{c}}\frac{1}{\sqrt{2n}+\sqrt{2n^{\prime}}}\frac{1}{n!n^{\prime}!}\left(\frac{\bar{q}_{1}\bar{q}_{2}^{*}}{2}\right)^{n}\left(\frac{\bar{q}_{1}^{*}\bar{q}_{2}}{2}\right)^{n^{\prime}}$
$\displaystyle=\int_{0}^{\infty}\int_{0}^{\infty}\frac{dq_{1}dq_{2}}{(2\pi)^{2}}\int_{0}^{2\pi}d\theta\
q_{1}e^{-q_{1}^{2}-q_{2}^{2}-q_{1}q_{2}\cos\theta}L_{m}(q_{1}^{2}+q_{2}^{2}-2q_{1}q_{2}\cos\theta)\sideset{}{{}^{\prime}}{\sum}_{n,n^{\prime}=0}^{N_{c}}\frac{1}{\sqrt{2n}+\sqrt{2n^{\prime}}}\frac{1}{n!n^{\prime}!}\left(\frac{q_{1}q_{2}}{2}\right)^{n+n^{\prime}}e^{i(n^{\prime}-n)\theta},$
(S41) $\displaystyle\frac{c_{m,b}^{\perp}}{4\pi}$
$\displaystyle=-2\iint\frac{d^{2}\bm{q}_{1}d^{2}\bm{q}_{2}}{(2\pi)^{4}}\frac{2\pi}{|\bm{q}_{2}|}e^{-\bm{q}_{1}^{2}-\bm{q}_{2}^{2}+\bm{q}_{1}\cdot\bm{q}_{2}}L_{m}(|\bm{q}_{1}-\bm{q}_{2}|^{2})\sum_{n,n^{\prime}=1}^{N_{c}}\frac{1}{\sqrt{2n}+\sqrt{2n^{\prime}}}\frac{\sqrt{nn^{\prime}}}{n!n^{\prime}!}\left(\frac{\bar{q}_{1}\bar{q}_{2}^{*}}{2}\right)^{n}\left(\frac{\bar{q}_{1}^{*}\bar{q}_{2}}{2}\right)^{n^{\prime}}\left(\frac{\bm{q}_{1}^{2}}{2}\right)^{-1}$
$\displaystyle=-4\int_{0}^{\infty}\int_{0}^{\infty}\frac{dq_{1}dq_{2}}{(2\pi)^{2}}\int_{0}^{2\pi}d\theta\
\frac{1}{q_{1}}e^{-q_{1}^{2}-q_{2}^{2}-q_{1}q_{2}\cos\theta}L_{m}(q_{1}^{2}+q_{2}^{2}-2q_{1}q_{2}\cos\theta)\sum_{n,n^{\prime}=1}^{N_{c}}\frac{1}{\sqrt{2n}+\sqrt{2n^{\prime}}}\frac{\sqrt{nn^{\prime}}}{n!n^{\prime}!}\left(\frac{q_{1}q_{2}}{2}\right)^{n+n^{\prime}}e^{i(n^{\prime}-n)\theta},$
(S42)

with $\bar{q}_{i}\equiv q_{i,x}+iq_{i,y}$.

$\displaystyle c_{m,c}^{0/z}=-\text{Re}$
$\displaystyle\int_{0}^{\infty}\int_{0}^{\infty}\frac{dq_{1}dq_{2}}{(2\pi)^{2}}\int_{0}^{2\pi}d\theta\
q_{1}e^{-q_{1}^{2}-\frac{q_{2}^{2}}{2}-iq_{1}q_{2}\sin\theta}L_{m}(q_{1}^{2})\left[\sideset{}{{}^{\prime}}{\sum}_{n^{\prime}\geq
n\geq
0}\frac{1}{\sqrt{2n}+\sqrt{2n^{\prime}}}\frac{1}{2^{n^{\prime}}n^{\prime}!}q_{1}^{n^{\prime}-n}q_{2}^{n+n^{\prime}}e^{i(n^{\prime}-n)\theta}\right.$
$\displaystyle\times\left.\left(L_{n}^{n^{\prime}-n}(\frac{q_{1}^{2}}{2})\mp
L_{n-1}^{n^{\prime}-n}(\frac{q_{1}^{2}}{2})\right)+\sum_{n>n^{\prime}\geq
0}\frac{1}{\sqrt{2n}+\sqrt{2n^{\prime}}}\frac{1}{2^{n}n!}(-q_{1})^{n-n^{\prime}}q_{2}^{n+n^{\prime}}e^{i(n^{\prime}-n)\theta}\left(L_{n^{\prime}}^{n-n^{\prime}}(\frac{q_{1}^{2}}{2})\mp
L_{n^{\prime}-1}^{n-n^{\prime}}(\frac{q_{1}^{2}}{2})\right)\right]$
$\displaystyle=-\text{Re}$
$\displaystyle\int_{0}^{\infty}\frac{dq_{1}}{2\pi}q_{1}e^{-q_{1}^{2}}L_{m}(q_{1}^{2})\sideset{}{{}^{\prime}}{\sum}_{n,n^{\prime}=0}^{N_{c}}\frac{1}{\sqrt{2n}+\sqrt{2n^{\prime}}}\frac{1}{2^{l}l!}q_{1}^{|n-n^{\prime}|}I_{n,n^{\prime}}(q_{1})\left(L_{j}^{|n-n^{\prime}|}(\frac{q_{1}^{2}}{2})\mp
L_{j-1}^{|n-n^{\prime}|}(\frac{q_{1}^{2}}{2})\right).$ (S43)

In the last line, we defined $j=\min(n,n^{\prime})$, $l=\max(n,n^{\prime})$,
and the integral Gradshteyn and Ryzhik (2014)

$I_{n,n^{\prime}}(q_{1})=\int_{0}^{\infty}dq_{2}q_{2}^{n+n^{\prime}}J_{|n-n^{\prime}|}(q_{1}q_{2})e^{-\frac{q_{2}^{2}}{2}}=\frac{2^{l+\frac{1}{2}}\Gamma\left(l+\frac{1}{2}\right)}{2^{|n-n^{\prime}|+1}\Gamma(|n-n^{\prime}|+1)}\,_{1}F_{1}\left(l+\frac{1}{2};|n-n^{\prime}|+1;-\frac{q_{1}^{2}}{2}\right)q_{1}^{|n-n^{\prime}|}.$
(S44)

Here, $\ {}_{1}F_{1}$ is a confluent hypergeometric function. The summations
in Eq. (S41) and (S43) exclude the $n=n^{\prime}=0$ case. Our numerical
results of Eqs. (S41)-(S43), summarized in Tables I(a) and (b) in the main
text, indicate that $c_{0,b}^{\mu}$and $c_{0,c}^{z}$ logarithmically diverge
with respect to the LL cutoff $N_{c}$ as $c_{0,a}^{0/z}$ do, whereas
$c_{0,c}^{0}$ and $c_{m\neq 0,b/c}^{\mu}$ converge.

Figure S4: The logarithmic dependence of $|c_{0}^{z,\perp}|$ on the LL cutoff
$N_{c}$. The blue dots and red triangles represent numerical results for
$c_{0}^{z}$ and $-c_{0}^{\perp}$, respectively. The dashed lines are linear
fit in the semi-log plot, which indicate $|c_{0}^{z,\perp}|\approx
0.11+0.5\ln{N_{c}}$.

The logarithmic divergence of the contact valley-dependent interactions in
graphene have been studied previously in the zero-field with the
renormalization group (RG) analysis Kharitonov (2012). It was shown that the
sublattice-and valley-dependent short-range interactions obey the RG flow
equations $d\bar{g}_{\alpha 0}/d\xi=0$ and $d\bar{g}_{\alpha
z}/d\xi=4F(w)(\bar{g}_{\alpha z}-\bar{g}_{\alpha\perp})$, where
$F(w)\approx\kappa/4$ for small $\kappa$, $\xi=\ln{l/a_{0}}$ and $l$ is the
running length scale. Ref. Kharitonov (2012) then estimates the effective
valley-dependent interactions in the zeroth LL in the magnetic field as $2\pi
l_{B}^{2}u_{\alpha}=\bar{g}_{\alpha 0}(l=l_{B})+\bar{g}_{\alpha z}(l=l_{B})$.
In the limit $\kappa\ln{l_{B}/a_{0}}\ll 1$, the RG equation yields

$\displaystyle 2\pi l_{B}^{2}u_{\alpha}$ $\displaystyle\approx g_{\alpha
0}+g_{\alpha z}(1+\kappa\ln\frac{l_{B}}{a_{0}})-\kappa
g_{\alpha\perp}\ln\frac{l_{B}}{a_{0}}$ (S45) $\displaystyle=g_{\alpha
0}+g_{\alpha z}(1+0.5\kappa\ln{N_{c}})-0.5\kappa
g_{\alpha\perp}\ln{N_{c}}+const.,$ (S46)

To arrive at the second line, we note that $N_{c}$ depends on the out-of-plane
magnetic field such that $\sqrt{N_{c}}/l_{B}\propto$ a fixed energy cutoff
$\Lambda$. The coefficient of the first-order-in-$\kappa\ln{N_{c}}$ term
agrees with our results, see Fig. S4. We note that the ratio $u_{z}/u_{\perp}$
is affected by both $\kappa$ and $l_{B}$. Thus, there can be a quantum phase
transition between the KD and AF states induced by either by out-of-plane
magnetic field (via $l_{B}$) or the dielectric screening environment (via
$\kappa$). In the main text, we primarily explore the transition induced by
changes in the dielectric environment. However, it’s important to mention that
a field-induced phase transition also occurs for a moderate fine structure
constant ($\kappa=0.5\sim 1$). Specifically, in a strong field scenario
($N_{c}=100$), the system resides in the AF state, whereas in a weak field
($N_{c}=7000$), it transitions to the KD state. Due to the logarithmic
dependence of $u_{z}$ and $u_{\perp}$ on the magnetic field, adjusting the
dielectric environment is a more effective means of modulating the phase
boundary than varying the magnetic field strength.

Here we have taken a common UV cutoff for all the short-range interactions
comparable to the lattice constant. This is appropriate for those that
originate from the Coulomb forces. For those arising from optical phonons it
would be more appropriate to take the UV cutoff to be comparable to the phonon
gap. However, since the phonon gap is smaller than the band-width of the
$p_{z}$ orbitals, accounting for this would mean that the contribution to
$u_{\alpha}$ coefficients from the interactions mediated by phonons (i.e.
$g_{\perp z}^{0},g_{z\perp}^{0}$ as discussed in Eq.(S23)), would have a
smaller enhancement factor due to renormalizations than the one we have used
here, because the phonon interactions would effectively renormalize over a
shorter RG time. However, since the phonon contribution to interaction
coefficients is already smaller than the Coulomb contribution, accounting for
this would not significantly alter our estimates. For example, this would mean
that the coefficient $c_{0}^{z}$ that accompanies the $K$-optical-phonon
potential $g_{\perp z}$ contribution to $u_{\perp}$ would be smaller than the
one listed before (e.g. in Eq.(S45)), because it would have smaller UV cutoff.
However, this contribution to the linear in $\kappa$ enhancement of
$u_{\perp}$ is smaller compared to that coming from the coefficient
$g^{0}_{\perp\perp}$ generated by the Coulomb interactions. Thus we expect
accounting for this would not significantly alter our estimates.

###  II.2 $\hat{V}_{aa}-$ induced Haldane pseudopotential corrections to
$u_{\alpha}(m)$

Figure S5: Feynman diagrams representing the quadratic-in-$g_{\alpha\mu}^{0}$
corrections to two-body valley-dependent Haldane pseudopotentials in the
zeroth LL. The dashed lines stand for the contact sublattice-and valley-
dependent interactions.

The non-vanishing valley-dependent interactions generated $\hat{V}_{aa}$ can
be represented by Feynman diagrams depicted in Fig. S5a (particle-particle), b
(double-exchange), c (vertex correction), and d (bubble). The first three
types of Feynman diagrams can be evaluated according to the same procedure
outlined in the previous section and we will explain later the details to
compute the bubble diagrams, which cannot be generated by $\hat{V}_{sa}$. The
results show that the main contributions from $\hat{V}_{aa}$ are the
corrections to $m=0$ Haldane pseudopotentials that diverge as $\sqrt{N_{c}}$,

$\displaystyle u_{\perp}^{(2)^{\prime}}(m=0)=$ $\displaystyle\frac{1}{2\pi
l_{B}^{2}\bar{g}}\left[g_{zz}^{0}g_{\perp z}^{0}z_{a}-g_{zz}^{0}g_{\perp
z}^{0}z_{b}+2(g_{zz}^{0}g_{\perp\perp}^{0}+g_{z\perp}^{0}g_{\perp
z}^{0})y_{b}-4g_{z\perp}^{0}g_{\perp\perp}^{0}x_{b}-g_{zz}^{0}g_{\perp
z}^{0}z_{c}-2g_{z\perp}^{0}g_{\perp z}^{0}y_{c}-(g_{\perp
z}^{0})^{2}z_{d}\right],$ (S47) $\displaystyle u_{z}^{(2)^{\prime}}(m=0)=$
$\displaystyle\frac{1}{2\pi l_{B}^{2}\bar{g}}\left[(g_{\perp
z}^{0})^{2}z_{a}-(g_{\perp z}^{0})^{2}z_{b}+4g_{\perp
z}^{0}g_{\perp\perp}^{0}y_{b}-4(g_{\perp\perp}^{0})^{2}x_{b}+g_{zz}^{0}(g_{zz}^{0}-2g_{\perp
z}^{0})z_{c}+2g_{zz}^{0}(g_{z\perp}^{0}-2g_{\perp\perp}^{0})y_{c}-(g_{zz}^{0})^{2}z_{d}\right],$
(S48)

where we have kept $g_{0\mu}=g_{\alpha 0}=0$ for simplicity, $\bar{g}=\hbar
v_{F}l_{B}/\sqrt{N_{c}}$, and the coefficients $z_{a}\approx 0.11,z_{b}\approx
y_{b}\approx x_{b}\approx 0.1,z_{c}=0.2,y_{c}=0.017,$ and $z_{d}=0.45$. The
subscripts of these coefficients represent the associated Feynman diagrams.
For the bubble diagram,

$\displaystyle\frac{u_{\alpha,d}^{(2)\prime}(m)}{2}$
$\displaystyle=-4\times\frac{\left(g_{\alpha z}^{0}\right)^{2}}{2\pi
l_{B}^{2}}\int\frac{d^{2}q}{(2\pi)^{2}}\prescript{}{r}{\langle
m|}e^{i\bm{q}\cdot\bm{\eta}_{ij}}|m\rangle_{r}\sideset{}{{}^{\prime}}{\sum}_{n,n^{\prime}=-\infty}^{\infty}\frac{\Theta(\epsilon_{n})-\Theta(\epsilon_{n^{\prime}})}{\epsilon_{n}-\epsilon_{n^{\prime}}}|F_{n^{\prime},n}^{z}(\bm{q})F_{0,0}^{z}(\bm{q})|^{2}$
$\displaystyle=-\frac{1}{2\pi l_{B}^{2}}\frac{2\left(g_{\alpha
z}^{0}\right)^{2}}{\hbar
v_{F}l_{B}}\int\frac{d^{2}\bm{q}}{(2\pi)^{2}}e^{-\frac{3\bm{q}^{2}}{2}}L_{m}(\bm{q}^{2})\sideset{}{{}^{\prime}}{\sum}_{n,n^{\prime}=0}^{N_{c}}\frac{1}{\sqrt{2|n|}+\sqrt{2|n^{\prime}|}}\frac{j!}{l!}\left(\frac{\bm{q}^{2}}{2}\right)^{l-j}\left(L_{j}^{l-j}(\frac{\bm{q}^{2}}{2})+L_{j-1}^{l-j}(\frac{\bm{q}^{2}}{2})\right)^{2}.$
(S49)

where $j=\min(n,n^{\prime})$, $l=\max(n,n^{\prime})$, and the summation
excludes the $n=n^{\prime}=0$ case. Our numerical calculations suggest that
$u_{\alpha,d}^{(2)^{\prime}}(m)\approx-z_{d}\frac{\left(g_{\alpha
z}^{0}\right)^{2}}{2\pi l_{B}^{2}\bar{g}}$ for large $N_{c}$ with $z_{d}$
converging to $\sim 0.45$. $u_{\alpha,d}^{(2)^{\prime}}(m>0)$ is at least one
order of magnitude smaller than $u_{\alpha,d}^{(2)^{\prime}}(m=0)$ and thus
the valley-dependent interactions remains short-ranged, although we would like
to mention that
$|u_{\alpha,d}^{(2)^{\prime}}|\gg|u_{\alpha,a/b/c}^{(2)^{\prime}}|$ for $m>0$,
indicating that the effective interactions generated by the bubble diagrams
can have a longer interaction range than the other three channels.

We comment that the square root divergence of
$\delta{u}_{\perp,z}^{\prime}(0)$ causes significant uncertainties on the
estimation of their values. For instance, if we fix $g_{v}N_{c}/2\pi
l_{B}^{2}=$number of graphene unit cell per area ($g_{v}=2$ is the valley
degeneracy), $\bar{g}=\hbar v_{F}a_{0}\sqrt{\sqrt{3}/2\pi}\approx
85\textrm{meV}\cdot\textrm{nm}^{2}$ and, according to the estimated strength
of bare valley-dependent interactions, Eq. (S24), $2\pi
l_{B}^{2}u_{\perp}^{(2)\prime}(m=0)\sim 159\textrm{meV}\cdot\textrm{nm}^{2}$
and $2\pi
l_{B}^{2}u_{z}^{(2)\prime}(m=0)\sim-492\textrm{meV}\cdot\textrm{nm}^{2}$. On
the other hand, if $\hbar v_{F}\sqrt{2N_{c}}/l_{B}=\Lambda\sim 2$eV, the
energy cut off of the continuum model of graphene, $2\pi
l_{B}^{2}u_{\perp}^{(2)\prime}(m=0)\sim 44\textrm{meV}\cdot\textrm{nm}^{2}$
and $2\pi
l_{B}^{2}u_{z}^{(2)^{\prime}}(m=0)\sim-137\textrm{meV}\cdot\textrm{nm}^{2}$.
Despite these uncertainties on the estimates, we can see that
$u_{\alpha}^{(2)^{\prime}}(m=0)$ are overall smaller than
$u_{\alpha}^{(2)}(m=0)$ depicted in Fig.1 for large $\kappa>0.5$, although
their ratio is much larger than the ratio of the bare lattice-scale
interaction and long-range Coulomb interaction,
$\mathcal{O}(g_{\alpha\mu}^{0}\epsilon/e^{2}l_{B})$, due to the $\sqrt{N_{c}}$
divergence in Eqs. (S47) and (S48). More importantly, we emphasize that unlike
$u_{\alpha}^{(2)}(m=0)$ generated by $\hat{V}_{sa}^{(2)}$,
$u_{\alpha}^{(2)^{\prime}}(m=0)$ are insensitive to the long-wave dielectric
screening (i.e., $\kappa$), which justifies neglecting
$u_{\alpha}^{(2)^{\prime}}(m=0)$ in the study of the dielectric screening
effects on the valley-dependent Haldane pseudopotentials in the main text.