# Global phase diagram of charge neutral graphene in the quantum Hall regime
for generic interactions

Suman Jyoti De<EMAIL_ADDRESS>Harish-Chandra Research Institute, A CI
of Homi Bhabha National Institute, Chhatnag Road, Jhunsi, Prayagraj 211019,
India  Ankur Das<EMAIL_ADDRESS>Department of Condensed Matter
Physics, Weizmann Institute of Science, Rehovot, 76100 Israel  Sumathi Rao
<EMAIL_ADDRESS>International Centre for Theoretical Sciences (ICTS-
TIFR), Shivakote, Hesaraghatta Hobli, Bangalore 560089, India  Ribhu K. Kaul
<EMAIL_ADDRESS>Department of Physics, The Pennsylvania State University,
University Park, PA 16802, USA  Ganpathy Murthy<EMAIL_ADDRESS>Department
of Physics & Astronomy, University of Kentucky, Lexington, KY 40506, USA

###### Abstract

Monolayer graphene at charge neutrality in a quantizing magnetic field is a
quantum Hall ferromagnet. Due to the spin and valley (near) degeneracies,
there is a plethora of possible ground states. Previous theoretical work,
based on a stringent ultra short-range assumption on the symmetry-allowed
interactions, predicts a phase diagram with distinct regions of spin-
polarized, canted antiferromagnetic, inter-valley coherent, and charge density
wave order. While early experiments suggested that the system was in the
canted antiferromagnetic phase at a perpendicular field, recent scanning
tunneling studies universally find Kekulé bond order, and sometimes also
charge density wave order. Recently, it was found that if one relaxes the
stringent assumption mentioned above, a phase with coexisting canted
antiferromagnetic and Kekulé order exists in the region of the phase diagram
believed to correspond to real samples. In this work, starting from the
continuum limit appropriate for experiments, we present the complete phase
diagram of $\nu=0$ graphene in the Hartree-Fock approximation, using generic
symmetry-allowed interactions, assuming translation invariant ground states up
to an intervalley coherence. Allowing for a sublattice potential (valley
Zeeman coupling), we find numerous phases with different types of coexisting
order. We conclude with a discussion of the physical signatures of the various
states.

##  I Introduction

The quantum Hall effects (QHE) [1, 2], discovered four decades ago [3] in
semiconductor heterostructures, embody many phenomena observed there for the
first time, but later found in many systems. The QHE represents the first and
simplest topological insulator [4, 5], as a consequence of which the electric
and thermal Hall conductances are quantized. The bulk is insulating; charge
and heat are carried by edge modes [6] which are robust against disorder. Due
to the quantization of kinetic energy into discrete values, Landau levels are
also the first example of truly flat bulk bands. As a result, the bulk physics
is controlled entirely by electron-electron interactions in a partially filled
Landau level. Notably, this leads to the fractional QHE (FQHE) states [7],
which host excitations with fractional charge and statistics [8]. Due to spin
or other internal degeneracies (such as valley or layer), ground states at
some integer fillings are also controlled by interactions. Typically,
interactions lead to ferromagnetism, as exemplified by the single layer
$\nu=1$ spin ferromagnet [9] or the bilayer $\nu=1$ state in $GaAs$ quantum
wells [10, 11, 12]. Such quantum Hall ferromagnets also have interesting
topological charged excitations such as skyrmions [9] or merons [11, 12].

Graphene [13, 14, 15, 16] is a single layer of Carbon atoms arranged in a
honeycomb lattice with two sites ($A$ and $B$) in each unit cell. Near charge
neutrality, low-energy electrons in graphene occur in two valleys at the two
inequivalent zone corners ${\mathbf{K}}$ and ${\mathbf{K}}^{\prime}$, and obey
a Dirac equation in each valley. In a quantizing perpendicular magnetic field
$B$, the Dirac spectrum produces particle-hole symmetric Landau levels
$n=0,\pm 1,\pm 2\cdots$ with energy $E_{\pm n}\propto\pm\sqrt{B|n|}$ in each
valley (ignoring the Zeeman splitting). Each Landau level is (almost) four-
fold degenerate, with the four states being labelled by spin and valley. The
$n=0$ Landau level is special; states in each valley are localized on one
sublattice. Furthermore, since the states are comprised of a superposition of
particle-like and hole-like states of the $B=0$ problem, there is valley
mixing at the edge, leading to one edge mode with a particle-like dispersion
and one with a hole-like dispersion [17].

It was recognized very early [18, 19, 20, 21]that partial filling of the $n=0$
manifold of Landau levels (called the zero-LLs or ZLLs) would lead to a rich
set of quantum Hall ferromagnets at integer fillings, the most interesting of
which is the $\nu=0$ case, when two of the four ZLLs are filled. The
particular superpositions of the spin and valley LLs that are occupied
determine the nature of the ground state.

In the non-interacting limit, the orbital part of the kinetic energy has an
$SU(4)$ spin/valley symmetry in the continuum limit. The Zeeman coupling
$E_{Z}$ splits the spin $\uparrow$ and spin $\downarrow$ Landau levels. The
ground state is fully spin-polarized [22], with the ${\mathbf{K}}\uparrow$,
and ${\mathbf{K}}^{\prime}\uparrow$ Landau levels occupied. The chemical
potential lies between the $\uparrow$-spin and $\downarrow$-spin sets of ZLLs.
Due to the nature of the edge states of the ZLLs mentioned above [17], at each
edge, the chemical potential intersects a particle-like level and a hole-like
level of opposite spin. Since the system is symmetric under $U(1)$ spin
rotations around the total $B$ field, these two modes cannot back-scatter,
even if potential disorder is present. Therefore the system should be in a
quantum spin Hall (QSH) phase [23, 4, 22].

Now let us add electron-electron interactions. The dominant long-range part of
the Coulomb interaction does not discriminate between spin and valley indices,
being $SU(4)$ symmetric. The Coulomb interaction changes the single-particle
spectrum at the edge, but leaves the two gapless, opposite-spin, charge-
carrying edge modes intact, preserving the QSH phase [24]. In transport, such
a state should show a two-terminal conductance of $2e^{2}/h$.

Initial experiments [25, 26, 27, 28] saw a trivial insulating state at $\nu=0$
without any protected edges. A seminal experiment measured the two-terminal
conductance in tilted field [29], which allowed independent tuning of the
Zeeman coupling. It was found that at large $E_{Z}$, the two-terminal
conductance does indeed tend asymptotically to $2e^{2}/h$. Below a critical
Zeeman energy $E_{Z}^{*}$, the system remains a trivial insulator. The gap at
the edge vanishes continuously as one approaches $E_{Z}\to E_{Z}^{*}$,
indicating a second-order phase transition.

The fact that the ground state at purely perpendicular field is not a QSH
state means that interactions beyond the $SU(4)$-symmetric Coulomb
interactions must play an important role [18, 19, 20, 21]. The ground state
must be chosen by $SU(4)$-anisotropic residual interactions, arising from
lattice-scale couplings. Low-energy effective symmetries inherited from the
$B=0$ problem can be used to deduce a $U(1)$ valley symmetry at the four-Fermi
level (reduced to a $Z_{3}$ symmetry when higher-Fermi interactions are
included) [18]. One can classify the four-Fermi anisotropic residual
interactions into two types; an Ising-like coupling in the valley space
$v_{z}({\mathbf{q}})$ and an $xy$-like coupling in the valley space
$v_{xy}({\mathbf{q}})$. See Section II.1 for the full definition of
$v_{z}({\mathbf{q}})$, $v_{xy}({\mathbf{q}})$ and the interaction Hamiltonian.

Since the residual interactions arise from lattice-scale couplings, and the
magnetic length $\ell=\sqrt{\frac{\hbar}{eB}}$ is two orders of magnitude
larger than the lattice spacing ($\ell\gg a$), the bare interactions are
ultra-short-range (USR). The first step in obtaining the effective
interactions in the ZLL manifold is to project the bare interactions to this
manifold. Operationally, this assumption implies that $v_{\mu}$ are
independent of ${\mathbf{q}}$ in the ZLL. It should be noted that if one
projects USR bare interactions to a $N\neq 0$ manifold in graphene, the
effective interactions will not be USR [30].

Using the USR assumption for the interactions and building upon previous work
[18, 19, 20, 21], Kharitonov [31] found the phase diagram in the Hartree-Fock
approximation for $\nu=0$ graphene. There are four phases: A fully polarized
phase $F$, an antiferromagnetic phase $AF$ (which becomes a canted AF, or CAF
in the presence of the Zeeman coupling), a charge density wave (CDW) phase,
and a phase with inter-valley coherence (IVC), sometimes also called a Kekulé
distorted or bond-ordered (KD/BO) phase. All the phase boundaries are first-
order, except for the CAF to F transition, which is second-order. Upon the
addition of a valley Zeeman coupling, the CDW and bond order coexist, leading
to a partially sublattice polarized (PSP) phase [32], but the transition
between the PSP and the CAF phases remains first-order.

This picture indeed reproduces the phenomenology of $\nu=0$ graphene in tilted
field [29], assuming that the anisotropic couplings are such that the system
(in perpendicular field) is in the CAF phase [31]. Subsequent magnon
transmission experiments [33, 34, 35, 36] through a $\nu=0$ region surrounded
by ferromagnetic $\nu=1$ regions also confirm that coherently propagating
magnetic excitations are present in the $\nu=0$ state. In the CAF phase the
natural candidate is the gapless Goldstone mode associated with the
spontaneous symmetry breaking of the $U(1)$ spin symmetry.

More recently, three scanning tunneling studies [37, 38, 39] on $\nu=0$
graphene perpendicular $B$ field have introduced more complexity into this
picture. While they cannot directly confirm or rule out CAF order, all three
see Kekulé bond order, and two of the three [38, 39] see CDW order as well.
Partial alignment of the graphene layer with the hexagonal Boron Nitride (HBN)
substrate is known [40, 41, 42, 43] to produce a sublattice potential/valley
Zeeman coupling $E_{V}$, which can lead to CDW order. This origin of the CDW
order was confirmed in one of the experiments [38], while it remains unclear
in the other experiment [39].

It is possible that the samples used in the scanning tunneling experiments are
in a different phase from those used in transport experiments. However, the
most parsimonious explanation is that CAF and bond order coexist in all
samples, the CAF being undetected in STM experiments and the bond order being
undetected in transport.

Since there is no coexistence between CAF and BO in Kharitonov’s phase diagram
[31], one is led to re-examine the assumption of ultra-short-range
interactions. From the point of view of an effective model in the ZLLs, there
is no reason to assume any particular form for the symmetry-allowed
interactions, because integrating out high energy Landau levels will lead to
renormalizations in the form of the interactions [44, 45, 46, 47, 48]. In a
recent work, three of us [49] used this reasoning, generalizing the symmetry-
allowed interactions $v_{z}({\mathbf{q}}),\ v_{xy}({\mathbf{q}})$ to be
arbitrary functions of ${\mathbf{q}}$. This might seem to introduce an
infinite number of new couplings. Remarkably, in the continuum Hartree-Fock
(HF) approximation, assuming translation symmetry is preserved up to an inter-
valley coherence, only two independent numbers per type of coupling suffice
[49], namely the Hartree and the Fock couplings, defined as

$g_{\mu,H}=\frac{v_{\mu}({\mathbf{q}}=0)}{2\pi\ell^{2}};\ \
g_{\mu,F}=\int\frac{d^{2}q}{(2\pi)^{2}}v_{\mu}({\mathbf{q}})e^{-q^{2}\ell^{2}/2}$
(1)

where $\mu=z,xy$. For ultra-short-range (USR) interactions
$g_{\mu,H}=g_{\mu,F}$. Indeed, it was found [49] that given certain
inequalities between the Hartree and Fock couplings, there was a region of
couplings which showed generic coexistence between the CAF and BO orders. We
will refer to this as the B/CAF phase.

It should be noted that non-USR interactions have implicitly been introduced
earlier by Goerbig and collaborators in the context of effective nonlinear
sigma models for $\nu=\pm 1$ in graphene [50, 51]. In $\nu=-1$ for example,
only a single ZLL (some linear combination of the four possibilities) is
occupied. As in any ferromagnet with single occupancy, the many-body wave
function is completely antisymmetric in space. Thus, USR interactions cannot
contribute to the energy of the state and the entire physics is controlled by
the non-USR couplings [52]. The case of $\nu=1$ is related to that of $\nu=-1$
by particle-hole symmetry. More recently, the connection of the parameters of
the effective theory to the non-USR nature of the microscopic anisotropic
interactions was made explicit by Atteia and Goerbig [52], once again in the
context of $\nu=\pm 1$.

There is a different line of reasoning coming from Bernal-stacked bilayer
graphene (BLG) in the quantum Hall regime [53] which also leads to the non-USR
condition $g_{\mu,H}\neq g_{\mu,F}$. Assuming solely nearest neighbor
hoppings, the $B=0$ dispersion has quadratic band touchings at ${\mathbf{K}}$
and ${\mathbf{K}}^{\prime}$. However, upon including the symmetry-allowed
trigonal warping (a hopping between non-Bernal stacked sites in different
layers) the quadratic band touching reconstructs into four Dirac cones [54].
The inclusion of trigonal warping has a profound effect [54] on the eight-fold
(nearly) degenerate manifold of states near charge neutrality in a quantizing
$B$. The upshot is that the the symmetry-allowed interactions, when projected
into the low-energy manifold, now acquire structure on the scale of $\ell$ and
no longer satisfy the USR condition $g_{\mu,H}=g_{\mu,F}$. This fact is
crucial in obtaining phases which show coexistence between different kinds of
order in BLG at $\nu=0$ [53]. There is a deep analogy between MLG and BLG, as
we will see; the states that we will uncover in the full phase diagram of MLG
are identical to a subset of states found earlier in BLG [53].

The purpose of this paper is to find the complete HF phase diagram of
monolayer graphene in the continuum approximation, assuming translation
invariance up to an inter-valley coherence. Our previous work [49] was
motivated by the STM experiments [37, 38, 39], and confined to values of
couplings thought to apply to real samples. Furthermore, the valley Zeeman
coupling was ignored. We will explore the full phase diagram, in the presence
of nonzero $E_{Z},\ E_{V}$ for all possible
$g_{z,H},g_{z,F},g_{xy,H},g_{xy,F}$.

The majority of the results in the main text are for the case when the Hartree
and Fock parts of a given coupling have the same sign:

$\frac{g_{z,F}}{g_{z,H}}>0;\ \ \ \ \frac{g_{xy,F}}{g_{xy,H}}>0$ (2)

This seems natural for weak LL-mixing, when the renormalizations from
integrating out the higher energy states are expected to be small compared to
the bare values of the couplings. However, for strong LL-mixing, one may well
have situations when the Hartree and Fock parts of a given coupling have
opposite signs. We will present some interesting results in this case as well.

To give a brief preview of our results. We find three coexistence phases: (i)
The coexistence phase occurring near the BO/CAF boundary in the USR model,
which was found earlier [49]. This phase, which we call the B/CAF phase, also
has a spin-valley entangled order even at $E_{V}=0$, which we label as SVEY
(we will explain the notation shortly in Section II.2). (ii) A phase occuring
near the CDW/FM boundary in the USR model. This phase displays the coexistence
between CDW and FM order, mediated by a spin-valley entangled order we label
SVE+ (explained in Section II.2). (iii) When $\frac{g_{z,F}}{g_{z,H}}<0$ we
find a phase where FM and SVEX/SVEY order (explained in Section II.2) coexist,
without any other order being present at $E_{V}=0$. We call this the FSVE
phase. When $E_{V}>0$ many of these phases acquire a CDW order parameter, but
remain largely unchanged otherwise. In sum, the full phase diagram of
monolayer graphene for generic interactions is much richer than was previously
believed.

The plan of the paper is as follows: In Section II we will briefly review the
previous work on the effective model for MLG in the continuum approximation.
We will generalize the interactions to be non-USR, discuss the HF
approximation and find the ground state energy. Also in Section II, we present
a parameterization of translation-invariant $\nu=0$ states [55, 51, 56] whose
energy depends on four angles. It turns out that the states that have been
found in the USR limit [31] can be characterized in terms of a single angle.
Instabilities of these states, which can be computed analytically, will signal
the occurrence of more complex states with coexisting order parameters. In
general, the actual ground state in any region of coexistence has to be found
numerically. In Section III we present our results; since there are six
independent tuning parameters, we will present many two-dimensional sections
through the space of coupling constants. Each two-dimensional section will
satisfy different inequalities between the Hartree and Fock couplings. We end
with our summary, conclusions, and open questions in Section IV. The
appendices contain the details of our calculations, analytical expressions for
various instabilities, and sample results for strong Landau-level mixing.

##  II Model Hamiltonian, parameters, and methods

We choose the primitive translation vectors for graphene as
${\mathbf{a}}_{1}=a{\hat{e}}_{x}$,
${\mathbf{a}}_{2}=a(\frac{{\hat{e}}_{x}}{2}+\frac{\sqrt{3}{\hat{e}}_{y}}{2})$,
with the general Bravais lattice site
${\mathbf{R}}=n_{1}{\mathbf{a}}_{1}+n_{2}{\mathbf{a}}_{2}$. The noninteracting
Hamiltonian of graphene at zero magnetic field, suppressing the spin index for
the moment, is

$\displaystyle H_{0}$
$\displaystyle=-t\sum\limits_{{\mathbf{R}},j}c^{\dagger}_{A{{\mathbf{R}}}}c_{B{{\mathbf{R}}}+{\bf
d}_{j}}+h.c.$ (3)

where $t$ is the nearest-neighbor hopping matrix element,
$c_{A{\mathbf{R}}},c_{B{\mathbf{R}}}$ destroy electrons at the $A$ and $B$
sublattice sites of the Bravais site ${\mathbf{R}}$, the sum on $j=1,2,3$ with
${\bf d_{1}}=0,{\bf d_{2}}={\mathbf{a}}_{1}-{\mathbf{a}}_{2},{\bf
d_{3}}=-{\mathbf{a}}_{2}$ and $h.c.$ stands for hermitian conjugate. Note that
there is no spin-orbit coupling in the Hamiltonian of Eq. 3. First-principles
[57, 58] and tight-binding calculations [59, 60] show that the spin-orbit
coupling in graphene is of the order of tens of $\mu eV$, smaller than any
other energy scale in the problem. We will set the spin-orbit coupling to zero
here and henceforth. Fourier transforming Eq. 3, we obtain the Bloch
Hamiltonian at wave-vector ${\mathbf{k}}$ as a matrix in the sublattice space

$H({\mathbf{k}})=-t\left(\begin{array}[]{cc}0&f({\mathbf{k}})\\\
f^{*}({\mathbf{k}})&0\end{array}\right)$ (4)

where
$f({\mathbf{k}})=1+e^{i{\mathbf{k}}\cdot({\mathbf{a}}_{1}-{\mathbf{a}}_{2})}+e^{-i{\mathbf{k}}\cdot{\mathbf{a}}_{2}}$.
It is easily checked that $f({\mathbf{k}})$ vanishes at the two inequivalent
zone corners (valleys)
${\mathbf{K}}={\hat{e}}_{x}\frac{4\pi}{3a}=-{\mathbf{K}}^{\prime}$ indicating
Dirac crossings. The low-energy effective Hamiltonian in the ${\mathbf{K}}$
valley (${\mathbf{k}}={\mathbf{K}}+{\mathbf{p}}$, $|{\mathbf{p}}|\ll 2\pi/a$)
can be obtained by expansion.

$H_{{\mathbf{K}}}({\mathbf{p}})=\frac{ta\sqrt{3}}{2}\left(\begin{array}[]{cc}0&p_{x}-ip_{y}\\\
p_{x}+ip_{y}&0\end{array}\right)$ (5)

The low-energy effective Hamiltonian at the ${\mathbf{K}}^{\prime}$ valley can
be obtained by the identity
$H_{{\mathbf{K}}^{\prime}}({\mathbf{p}})=(H_{{\mathbf{K}}}(-{\mathbf{p}}))^{*}$.

In the continuum limit, we turn on a weak perpendicular magnetic field
$B_{\perp}$ by allowing the Hamiltonian to act on slowly varying envelope
functions. Operationally, this involves promoting
$p_{i}\to-i\partial_{i}\to-i\partial_{i}+eA_{i}({\mathbf{r}})$, where the
electron’s charge is $-e$, and ${\mathbf{A}}({\mathbf{r}})$ is the vector
potential satisfying $\nabla\times{\mathbf{A}}={\hat{e}}_{z}B_{\perp}$. In
order for the continuum limit to be justified, the magnetic length has to be
much larger than the lattice spacing; $\ell=\sqrt{\frac{\hbar}{eB_{\perp}}}\gg
a$. This is extremely well-satisfied for realistic fields.

Next, one chooses Landau gauge ${\mathbf{A}}=B_{\perp}x{\hat{e}}_{y}$, and
imposes periodic boundary conditions in the $y$-direction with a length
$L_{y}$. Let us define the Landau level wavefunctions as

$\langle
x,y|n,k\rangle=\frac{e^{iky}}{\sqrt{L_{y}}}\Phi_{n}\left(\frac{x-k\ell^{2}}{\ell}\right)$
(6)

where $\Phi_{n}$ are the normalized harmonic oscillator wavefunctions. Note
that, here and henceforth, $k$ is a one-dimensional guiding center label, and
not a two-dimensional momentum. Now it is straightforward to see that the ZLL
states are $\left(0,|n=0,k\rangle\right)^{T}$ in the ${\mathbf{K}}$ valley,
and $\left(|n=0,k\rangle,0\right)^{T}$ in the ${\mathbf{K}}^{\prime}$ valley.
Thus, in the ZLLs, valley and sublattice are locked together. Now we are ready
to present our model Hamiltonian.

###  II.1 Hamiltonian and the Hartree-Fock Approximation

In what follows, we will index the fermion operators with a valley index
$\alpha,\beta$, which can be ${\mathbf{K}}\equiv 0$ or
${\mathbf{K}}^{\prime}\equiv 1$, and a spin index $s=\uparrow\equiv 0$ or
$s=\downarrow\equiv 1$. In this notation, the non-interacting Hamiltonian of
the ZLLs is

$H_{1b}=-\sum\limits_{\alpha,s,k}\left(E_{Z}(-1)^{s}+E_{V}(-1)^{\alpha}\right)c^{\dagger}_{\alpha,s,k}c_{\alpha,s,k}$
(7)

which introduces the Zeeman energy $E_{Z}$ and the valley Zeeman/sublattice
potential $E_{V}$.

Turning to interactions, Alicea and Fisher [18] noted that, in addition to the
$SU(4)$-symmetric Coulomb interaction, two other types of low-energy effective
interactions were allowed by $SU(2)$ spin-rotation symmetry and momentum
conservation in the $B=0$ problem. We recall that the spin-orbit coupling is
negligible [59, 57, 58, 60] and has been neglected. A $U(1)$ symmetry in the
valley space (separate conservation of the number of electrons in each valley)
emerges when restricting oneself to four-Fermi interactions. Upon including
six-Fermi terms this is reduced to a $Z_{3}$ symmetry because
$3({\mathbf{K}}-{\mathbf{K}}^{\prime})$ is a reciprocal lattice vector. The
full interaction Hamiltonian for the ZLLs in monolayer graphene, in the Landau
gauge basis discussed earlier, is

$\displaystyle H=$ $\displaystyle H_{1b}+H_{int}$ (8) $\displaystyle H_{int}=$
$\displaystyle H_{Coul}+H_{Z}+H_{xy}$ (9) $\displaystyle H_{Z}=$
$\displaystyle\frac{1}{2L_{x}L_{y}}\sum\limits_{k,k^{\prime},{\mathbf{q}}}v_{z}({\mathbf{q}})e^{-iq_{x}(k-k^{\prime}-q_{y})\ell^{2}}e^{-(q\ell)^{2}/2}:c^{\dagger}_{k-q_{y}}\tau_{Z}c_{k}c^{\dagger}_{k^{\prime}+q_{y}}\tau_{Z}c_{k^{\prime}}:$
(10) $\displaystyle H_{xy}=$
$\displaystyle\frac{1}{2L_{x}L_{y}}\sum\limits_{k,k^{\prime},{\mathbf{q}}}v_{xy}({\mathbf{q}})e^{-iq_{x}(k-k^{\prime}-q_{y})\ell^{2}}e^{-(q\ell)^{2}/2}\left(:c^{\dagger}_{k-q_{y}}\tau_{x}c_{k}c^{\dagger}_{k^{\prime}+q_{y}}\tau_{x}c_{k^{\prime}}:+:c^{\dagger}_{k-q_{y}}\tau_{y}c_{k}c^{\dagger}_{k^{\prime}+q_{y}}\tau_{y}c_{k^{\prime}}:\right)$
(11)

We have used a shorthand notation where the spin/valley indices are suppressed
and summed, and $\tau_{i}$ is a Pauli matrix in the valley space (it acts as
the identity in the spin space). Furthermore, since the Coulomb interaction
has no role to play in selecting the ground state, we drop it henceforth. In
the generic case, when $E_{Z},\ E_{V}>0$, the Hamiltonian has a $U(1)_{s}$
spin-rotation symmetry generated by total $\sigma_{z}$, a $U(1)_{v}$ valley-
rotation symmetry generated by total $\tau_{z}$ (which is also the difference
between the number of electrons in the ${\mathbf{K}}$ and
${\mathbf{K}}^{\prime}$ valleys), and an entangled spin-valley $U(1)_{sv}$
symmetry generated by total $\tau_{z}\sigma_{z}$. In the fine-tuned case
$E_{Z}=0$ (not realizable in experimental samples) the spin-rotation symmetry
is enhanced to $SU(2)_{s}$. In the fine-tuned case $E_{V}=0$, the valley
symmetry is enhanced to a $U(1)_{v}\otimes Z_{2v}$, where the $Z_{2v}$
represents the symmetry exchanging the two valleys.

In the HF approximation, one looks for the single Slater determinant that has
the right electron filling and minimizes the energy. Such a state,
symbolically written as $|HF\rangle$, can be completely characterized by the
set of 1-body expectation values.

$\Delta_{\alpha\beta}^{ss^{\prime}}(k,k^{\prime})=\langle
HF|c^{\dagger}_{\alpha,s,k}c_{\beta,s^{\prime},k^{\prime}}|HF\rangle$ (12)

We will restrict the space of HF states to those obeying translation
invariance, up to an intervalley coherence. This means that the $\Delta$
becomes diagonal in $k$ and independent of it.

$\Delta_{\alpha\beta}^{ss^{\prime}}(k,k^{\prime})=\delta_{kk^{\prime}}\Delta_{\alpha\beta}^{ss^{\prime}}$
(13)

Let us make the idea of translation invariance up to an intervalley coherence
more explicit. The $\Delta$ we have assumed allows nonzero averages of the
form $\langle
HF|c^{\dagger}_{{\mathbf{K}},s}c_{{\mathbf{K}}^{\prime}s^{\prime}}|HF\rangle$.
Clearly these break lattice translations since ${\mathbf{K}}$ and
${\mathbf{K}}^{\prime}$ are not identical up to a reciprocal lattice vector.
Allowing such averages introduces a new set of reciprocal lattice vectors
which are ${\mathbf{K}}-{\mathbf{K}}^{\prime}$ and all their rotated versions,
and leads to translation invariance with an enlarged unit cell of size
$\sqrt{3}\times\sqrt{3}$ as compared to the original. Indeed, this is exactly
what is seen STM experiments [37, 38, 39]. The ansatz of Eq. 13 makes sure
that there is no translation symmetry breaking beyond the minimal one implied
by intervalley coherence.

In general, the matrix $\Delta$ is the projector on to the linear space of the
occupied states. Given that two orthogonal linear combinations of the four
ZLLs (call them $|f_{1}\rangle$ and $|f_{2}\rangle$) are occupied at $\nu=0$
we can write

$\Delta=|f_{1}\rangle\langle f_{1}|+|f_{2}\rangle\langle f_{2}|$ (14)

We can now express the HF energy of the Hamiltonian of Eq. 11 per guiding
center in terms of $\Delta$, with $N_{\phi}=\frac{L_{x}L_{y}}{2\pi\ell^{2}}$,
as

$\displaystyle\frac{E_{HF}}{N_{\phi}}=$ $\displaystyle-
E_{Z}Tr\left[\sigma_{Z}\Delta\right]-E_{V}Tr\left[\tau_{Z}\Delta\right]+\frac{g_{z,H}}{2}\left(Tr[\tau_{Z}\Delta]\right)^{2}-\frac{g_{z,F}}{2}Tr\left[\tau_{Z}\Delta\tau_{Z}\Delta\right]$
$\displaystyle+\frac{g_{xy,H}}{2}\left\\{\left(Tr[\tau_{x}\Delta]\right)^{2}+\left(Tr[\tau_{y}\Delta]\right)^{2}\right\\}-\frac{g_{xy,F}}{2}\left(Tr[\tau_{x}\Delta\tau_{x}\Delta]+Tr[\tau_{y}\Delta\tau_{y}\Delta]\right)$
(15)

Note that $g_{\mu,h}$ and $g_{\mu,F}$ are defined as in Eq. 1.

###  II.2 Ansatz for States, Instabilities, and Order Parameters

We will start with an efficient parameterization [55, 51, 56] for the two
orthogonal occupied states $|f_{1}\rangle$ and $|f_{2}\rangle$. This
parameterization has been used not only for uniform states but also for
describing skyrmions [56].

$\displaystyle|f_{1}\rangle=$
$\displaystyle\cos{\frac{\alpha_{1}}{2}}|{\mathbf{n}}\rangle\otimes|{\mathbf{s}}\rangle+e^{i\beta_{1}}\sin{\frac{\alpha_{1}}{2}}|-{\mathbf{n}}\rangle\otimes|-{\mathbf{s}}\rangle$
(16) $\displaystyle|f_{2}\rangle=$
$\displaystyle\cos{\frac{\alpha_{2}}{2}}|{\mathbf{n}}\rangle\otimes|-{\mathbf{s}}\rangle+e^{i\beta_{2}}\sin{\frac{\alpha_{2}}{2}}|-{\mathbf{n}}\rangle\otimes|{\mathbf{s}}\rangle$
(17)

where
${\mathbf{n}}=\sin{\theta_{p}}\cos{\phi_{p}}{\hat{e}}_{x}+\sin{\theta_{p}}\sin{\phi_{p}}{\hat{e}}_{y}+\cos{\theta_{p}}{\hat{e}}_{z}$,
and
${\mathbf{s}}=\sin{\theta_{s}}\cos{\phi_{s}}{\hat{e}}_{x}+\sin{\theta_{s}}\sin{\phi_{s}}{\hat{e}}_{y}+\cos{\theta_{s}}{\hat{e}}_{z}$
are unit vectors indicating the directions of the state on the valley and spin
Bloch spheres respectively. The spinors $|{\mathbf{n}}\rangle$ and
$|{\mathbf{s}}\rangle$ are defined in the standard way

$|{\mathbf{n}}\rangle=\left(\begin{array}[]{c}\cos{\frac{\theta_{p}}{2}}\\\
e^{i\phi_{p}}\sin{\frac{\theta_{p}}{2}}\end{array}\right);\ \
|{\mathbf{s}}\rangle=\left(\begin{array}[]{c}\cos{\frac{\theta_{s}}{2}}\\\
e^{i\phi_{s}}\sin{\frac{\theta_{s}}{2}}\end{array}\right)$ (18)

In going from $|{\mathbf{n}}\rangle\to|-{\mathbf{n}}\rangle$ one substitutes
$\theta_{p}\to\pi-\theta_{p}$ and $\phi_{p}\to\phi_{p}+\pi$, and likewise for
${\mathbf{s}}$.

Given this ansatz, which depends on eight angles, we compute the HF energy.

$\displaystyle\text{E}_{HF}=$ $\displaystyle-
E_{Z}\cos\theta_{s}\left[\cos\alpha_{1}-\cos\alpha_{2}\right]-E_{V}\cos\theta_{p}[\cos\alpha_{1}+\cos\alpha_{2}]+\frac{g_{z,H}}{2}\cos^{2}\theta_{p}\left[\cos\alpha_{1}+\cos\alpha_{2}\right]^{2}$
$\displaystyle-\frac{g_{z,F}}{16}\left[4\cos^{2}\theta_{p}\left(\cos(2\alpha_{1})+\cos(2\alpha_{2})\right)-8\sin^{2}\theta_{p}\left(\cos\alpha_{1}\cos\alpha_{2}-\sin\alpha_{1}\sin\alpha_{2}\cos(\beta_{1}+\beta_{2})\right)+8\right]$
$\displaystyle+\frac{g_{xy,H}}{2}\sin^{2}\theta_{p}\left[\cos\alpha_{1}+\cos\alpha_{2}\right]^{2}-\frac{g_{xy,F}}{16}\bigg{[}8\sin^{2}\theta_{p}\left(\cos\alpha_{1}\cos\alpha_{2}-\sin\alpha_{1}\sin\alpha_{2}\cos(\beta_{1}+\beta_{2})\right)$
$\displaystyle+4\sin^{2}\theta_{p}\left(\cos(2\alpha_{1})+\cos(2\alpha_{2})\right)-16\left(\cos\alpha_{1}\cos\alpha_{2}-1\right)\bigg{]}$
(19)

There are three noteworthy features of this energy. The first is that it is
independent of $\phi_{p}$ and $\phi_{s}$. This results from the $U(1)_{v}$ and
$U(1)_{s}$ symmetries of the Hamiltonian. Therefore we can set
$\phi_{s}=\phi_{p}=0$ without loss of generality. The second feature is that
the dependence of the energy on $\beta_{1},\beta_{2}$ occurs only in the
interacting part, and only in the combination $\beta_{1}+\beta_{2}$. The
reason the one-body HF energy does not involve $\beta_{1},\beta_{2}$ is that
the averages of $\tau_{z},\ \sigma_{z}$ do not involve $\beta_{i}$.

$\displaystyle\langle f_{i}|\tau_{z}|f_{i}\rangle=$
$\displaystyle\cos{\alpha_{i}}\cos{\theta_{p}}$ (20) $\displaystyle\langle
f_{i}|\sigma_{z}|f_{i}\rangle=$
$\displaystyle-(-1)^{i}\cos{\alpha_{i}}\cos{\theta_{s}}$ (21)

The dependence of the interaction energy solely on $\beta_{1}+\beta_{2}$
arises from the $SU(2)_{spin}$ symmetry of the interactions, which implies
that the $U(1)$ rotation $|{\mathbf{s}}\rangle\rightarrow
e^{i\chi/2}|{\mathbf{s}}\rangle;\ |-{\mathbf{s}}\rangle\rightarrow
e^{-i\chi/2}|-{\mathbf{s}}\rangle$ cannot change the interaction energy.
Ignoring overall phase factors, this rotation has the net effect
$\beta_{1}\rightarrow\beta_{1}-\chi;\ \beta_{2}\rightarrow\beta_{2}+\chi$.
This demonstrates that the energy can only depend on $\beta_{1}+\beta_{2}$.
Thirdly, the dependence on $\beta_{1}+\beta_{2}$ occurs via the term
$\cos(\beta_{1}+\beta_{2})$, which appears linearly. Depending on the sign of
its coefficient, the minimum energy will occur at
$\cos(\beta_{1}+\beta_{2})=\pm 1$.

The bottom line is that the minimum of the HF energy for uniform states can be
found in a subspace in which $|f_{1}\rangle,|f_{2}\rangle$ can both be chosen
real.

Below, we will call states which have Kekulé/BO, and/or CDW order B/CO states
(because they have both bond order and/or charge order). The states originally
found by Kharitonov [31] can be represented in terms of the above angles as
follows (details in Appendix A),

$\displaystyle|FM\rangle=|\alpha_{1}=0,\alpha_{2}=\pi,\theta_{p}=\pi/2,\theta_{s}=0,\beta_{1}=\beta_{2}=\pi\rangle$
(22)
$\displaystyle|CAF\rangle=|\alpha_{1}=\theta_{CAF},\alpha_{2}=\pi-\alpha_{1},\theta_{p}=\pi/2,\theta_{s}=0,\beta_{1}=\beta_{2}=\pi\rangle;\
g_{xy,F}<0$ (23)
$\displaystyle|B/CO\rangle=|\alpha_{1}=\alpha_{2}=0,\theta_{p}=\theta_{B/CO},\theta_{s}=0,\beta_{1}=\beta_{2}=0\rangle$
(24)

where

$\displaystyle\theta_{CAF}=\cos^{-1}\left(\frac{E_{Z}}{2|g_{xy,F}|}\right);\ \
g_{xy,F}<0$ (25)
$\displaystyle\theta_{B/CO}=\cos^{-1}\left(\frac{E_{V}}{g_{V}}\right);\ \
g_{V}>E_{V}$ (26) $\displaystyle\theta_{B/CO}=0;\ \ g_{V}<E_{V}$ (27)
$\displaystyle g_{V}=2g_{z,H}-g_{z,F}-2g_{xy,H}+g_{xy,F}$ (28)

If $\theta_{B/CO}=0$ the state is a pure CDW, if $0<\theta_{B/CO}<\pi/2$ it
has coexisting BO and CDW order, and if $\theta_{B/CO}=\pi/2$ the system is in
a pure Kekulé state. The reason we call these states “simple” is that they can
all be described by at most a single nontrivial angle, which can be
analytically computed as a function of the couplings. Generic states may
depend on more than one nontrivial angle, in which case it is not possible to
solve for the angles or the ground state energy analytically.

The ground state energies of the simple states are

$\displaystyle E_{\text{FM}}=$ $\displaystyle-2(E_{Z}+g_{xy,F})-g_{z,F}$
$\displaystyle E_{\text{CAF}}=$
$\displaystyle\frac{E_{Z}^{2}}{2g_{xy,F}}-g_{z,F};\ 0<E_{Z}<-2g_{xy,F}$
$\displaystyle E_{\text{B/CO}}=$
$\displaystyle-\frac{E_{V}^{2}}{g_{V}}-g_{xy,F}+2g_{xy,H};\ 0<E_{V}<g_{V}$
$\displaystyle E_{\text{CDW}}=$ $\displaystyle-2E_{V}-g_{z,F}+2g_{z,H}$ (29)

Our strategy is to examine the stability of these “simple” states by finding
the eigenvalues of the Hessian matrix of second derivatives of the energy
functional with respect to the four angles
$\alpha_{1},\alpha_{2},\theta_{p},\theta_{s}$.

${\cal E}_{ij}=\frac{\partial^{2}E_{HF}}{\partial\chi_{i}\partial\chi_{j}}$
(30)

where $\chi_{i}$ represents all four angles. For the “simple” states one can
compute the entire Hessian matrix analytically, and also obtain the
eigenvalues analytically.

All eigenvalues being positive means the state is stable to arbitrary small
deformations. As the coupling constants are varied, a formerly positive
eigenvalue may vanish, indicating an instability of the state in question.
This allows us to map out the regions of stability of the “simple” states in
our coupling constant space.

It can happen that when some of the angles take particular values, the
projector on to the occupied subspace becomes independent of certain other
angles. This occurs in the FM and B/CO phases. Consequently, certain rows and
columns ${\cal E}_{ij}$ vanish, which means that one eigenvalue always
vanishes in that state independent of the coupling constants. In such cases,
the instability is marked by the vanishing of an eigenvalue that does depend
on coupling constants. Once the region of possible coexistence has been found,
we use numerical self-consistent Hartree-Fock to obtain the ground state and
confirm the coexistence predicted by the Hessian.

In preparation for showing the results, let us list all the order parameters
which we will encounter and the symmetries they break. We have chosen a
parameterization in which the projector matrix of the occupied states is real.
Thus, out of all possible hermitian matrices that represent order parameters,
only real matrices will have nonzero expectation values.

$\displaystyle FM=$ $\displaystyle\langle\sigma_{z}\rangle/2$ $\displaystyle
CAF=$ $\displaystyle\langle\tau_{z}\sigma_{x}\rangle/2$ $\displaystyle BO=$
$\displaystyle\langle\tau_{x}\rangle/2$ $\displaystyle CDW=$
$\displaystyle\langle\tau_{z}\rangle/2$ $\displaystyle SVEX=$
$\displaystyle\langle\tau_{x}\sigma_{x}\rangle/2$ $\displaystyle SVEY=$
$\displaystyle\langle\tau_{y}\sigma_{y}\rangle/2$ $\displaystyle SVE\pm=$
$\displaystyle\langle\tau_{x}\sigma_{x}\pm\tau_{y}\sigma_{y}\rangle/2$ (31)

The SVE (spin-valley entangled) type of order parameters are so called because
they break the spin and valley symmetries simultaneously in an entangled way.

Let us examine the symmetries broken by the various order parameters. In the
fine-tuned case $E_{Z}=0$ (not realizable in experiment) the FM order
parameter spontaneously breaks the $SU(2)_{s}$ symmetry. In the fine-tuned
case $E_{V}=0$ (which is potentially realizable in experiment) the CDW order
spontaneously breaks the $Z_{2v}$ symmetry. In the generic case $E_{Z},\
E_{V}\neq 0$ the FM and CDW order parameters do not break any symmetries of
the Hamiltonian. The CAF order parameter breaks $U(1)_{s}$ and $U(1)_{sv}$,
but preserves $U(1)_{v}$. Bond order breaks $U(1)_{v}$ and $U(1)_{sv}$ but
preserves $U(1)_{s}$. The SVE order parameters break $U(1)_{s}$ and $U(1)_{v}$
but preserve $U(1)_{sv}$. All three $U(1)$s are spontaneously broken in the
B/CAF phase.

We emphasize that while the $U(1)_{s}$ symmetry is protected by the assumed
vanishing of the spin-orbit coupling in graphene (Eq. 3), there is no such
protection for $U(1)_{v}$ or $U(1)_{sv}$. As mentioned in the introduction,
since $3({\mathbf{K}}-{\mathbf{K}}^{\prime})$ is a reciprocal lattice vector,
any $U(1)$ related to the valley space will be broken down to a $Z_{3}$ upon
including six-Fermi and higher interactions. Thus, there are no Goldstone
modes associated with the spontaneous breaking of the $U(1)_{v}$ symmetry.
However, a spontaneous breaking of the $U(1)_{s}$ and/or the $U(1)_{sv}$
symmetries will lead to a phase with a Goldstone mode.

##  III Results

There are six coupling constants in our Hamiltonian,
$E_{Z},E_{V},g_{z,H},g_{z,F},g_{xy,H},g_{xy,F}$. Since the full six-
dimensional phase diagram is impossible to visualize, we will be forced to
take two-dimensional cuts.

We will often draw a correspondence between MLG and Bernal-stacked BLG [61,
62]. In BLG, in addition to spin and valley, an orbital index $n=0,1$ also
appears [54]. For BLG states that have no orbital mixing, and are orbitally
symmetric, there is a one-to-one correspondence with states in MLG.
Specifially, the role of $E_{V}$ in MLG is played by the perpendicular
electric field $D$ in BLG [53].

While RG arguments tell us that generic effective interactions must have
nonzero range, one does not know precisely how the LL-mixing and the
intergration of high-energy states affects the Hartree and Fock parts of the
couplings in the ZLL manifold. For weak LL-mixing, a natural assumption is
that the sign of $g_{\mu,H},g_{\mu,F}$ are the same, but the magnitudes could
be different. The majority of the results we present in the main text assume
this to be true. For strong LL-mixing, it is conceivable that the effective
$g_{\mu,H}$ and $g_{\mu,F}$ have opposite signs. For completeness we have
analyzed this case as well, though most of the details are relegated to the
appendices.

In Section III.1, we will present results $E_{V}=0$, restricting the Hartree
and Fock parts of the couplings to have the same sign. This will allow us to
examine how the original Kharitonov phase diagram [31] changes when we relax
the USR assumption. We confirm the existence of the B/CAF phase found earlier
[49] for $g_{xy,F}/g_{xy,H}>1$. This phase also has the order parameter SVEY
(Eq. 31). Thus, this phase breaks $U(1)_{s}$, $U(1)_{v}$, and $U(1)_{sv}$
spontaneously. Another coexistence phase occurs near the boundary between the
FM and CDW phases when $g_{z,F}/g_{z,H}>1$. We call this the SVE+ phase
because the corresponding order parameter (Eq. 31) is nonzero in this phase.
This phase breaks $U(1)_{s}$ and $U(1)_{v}$ but preserves $U(1)_{sv}$. The
bilayer analog of this phase was found earlier in Bernal-stacked BLG [53],
although at nonzero $E_{V}$.

Next, in Section III.2 we will keep the restriction that the Hartree and Fock
parts of the couplings have the same sign, but turn on $E_{V}$. Already, for
USR interactions, one finds a phase with coexisting CDW and Kekulé order [56],
which we call the B/CO phase. Guided by the intuition that new phases are most
likely to appear near the transition lines of the original USR model, we will
examine these carefully. When $0>g_{xy,H}>g_{xy,F}$, the B/CAF phase shrinks
as $E_{V}$ increases (all other couplings remaining constant), and vanishes
for large enough $E_{V}$. For any $E_{V}>0$, both SVEX and SVEY order
parameters are nonzero in the B/CAF phase. However, the symmetries that are
spontaneously broken do not change in any of the phases at small $E_{V}$.

Going to larger $E_{V}$ is even more interesting: in certain cases, even if
coexistence does not occur for $E_{V}=0$, it can occur for intermediate
$E_{V}$, and disappear for large $E_{V}$. For example, for
$g_{z,F}=0.75g_{z,H},\ \ g_{xy,F}=0.75g_{xy,H}$, both the B/CAF and SVE+
phases are absent at $E_{V}=0$. However, both phases are present for a range
of intermediate $E_{V}$. This is consistent with earlier results for BLG [53],
which also found coexistence phases for intermediate values of the
perpendicular electric field $D$, which plays the same role in BLG as $E_{V}$
plays in MLG.

From the point of view of experiment, the most interesting phases is the B/CAF
phase, which occurs for $0>g_{xy,H}>g_{xy,F}$. At any nonzero $E_{V}$, the
B/CAF phase evolves into one with coexisting CAF, Kekulé, and CDW order. An
analogous phase breaking all three $U(1)$ symmetries was found earlier in BLG
[53]. (Of course, the $U(1)_{v}$ symmetry will be reduced to a $Z_{3}$ upon
the inclusion of higher-Fermi interactions, and thus does not lead to a
Goldstone mode). Such a phase would be consistent with observations of magnon
transmission [33, 35, 36] through the $\nu=0$ state, which implies some kind
of magnetic order, and the STM experiments [37, 38, 39] which show Kekulé and
CDW order.

Finally, in Section III.3 we will present some results for the signs of
$g_{\mu,H}$ and $g_{\mu,F}$ being different for either or both of $g_{z}$ and
$g_{xy}$. As may be expected, the topology of the phase diagram can change
considerably in such cases. However, a new type of coexistence phase appears,
which is not present when the signs of the H and F parts of both couplings are
identical. This phase, which we call the FSVE phase, has coexisting FM and
SVEY order, without any other order being present at $E_{V}=0$. Such a state
breaks $U(1)_{s}$ and $U(1)_{v}$ while preserving $U(1)_{sv}$. As usual, many
phases acquire nonzero CDW order when $E_{V}>0$. A more detailed set of
results for the H and F parts of $g_{xy/z}$ being of opposite signs are
presented in the appendices.

###  III.1 Vanishing Valley Zeeman Coupling with $g_{\mu.F}/g_{\mu,H}>0$

Throughout this section we will assume that we are “close” to the USR model in
the sense that the Hartree and Fock couplings of a given type have the same
sign ($g_{\mu,F}/g_{\mu,H}>0$), and that there is perfect sublattice symmetry
$E_{V}=0$. In order to enable comparisons with the USR model, we show the
Kharitonov phase diagram [31] in Fig. 1. Here we choose $E_{Z}=1.0$ (in
arbitrary units) as a fixed parameter. All the lines between phases are first-
order transitions, except for transition between the CAF and FM phases, which
is second-order. Our convention here and henceforth is that solid lines
represent first-order phase transitions while dashed lines represent
continuous phase transitions.

Figure 1: Parameters are
$E_{Z}=1.0,g_{z,F}=g_{z,H}=g_{z},g_{xy,F}=g_{xy,H}=g_{xy},E_{V}=0.0$,this is
the Kharitonov’s Ultra short range limit.Here as one can see their is no
coexistence phase.

When one relaxes the USR assumption, it turns out the inequalities
$g_{xy,F}/g_{xy,H}\lessgtr 1$ and $g_{z,F}/g_{z,H}\lessgtr 1$ play a crucial
role in determining whether coexistence occurs at $E_{V}=0$ (we will see later
that nonzero $E_{V}$ overcomes this limitation). Briefly, coexistence between
Kekulé and CAF order occurs only when $g_{xy,F}/g_{xy,H}>1$, while coexistence
between CDW and FM order occurs only when $g_{z,F}/g_{z,H}>1$. This is
illustrated in Figs. 2 and 3 for specific values of other couplings, and turns
out to be generic as long as $E_{V}=0$.

Figure 2: Instability lines for the BO and CAF phases as a function of
$g_{xy,H}$ and $g_{z,F}$ for $E_{Z}=1.0,g_{z,H}=2.0,g_{xy,F}=-1.0,E_{V}=0.0$.
The BO phase is unstable above the black line, while the CAF phase is unstable
below the green line. The red line is where the energies of the BO and CAF
phases cross. The vertical dashed line represents $g_{xy,H}=g_{xy,F}$. One can
see that coexistence between BO and CAF state only occurs when
$0>g_{xy,H}>g_{xy,F}$. The coexistence phase is shaded blue.  Figure 3:
Instability lines for the CDW and FM phases as a function of $g_{xy,F}$ and
$g_{z,H}$ for $E_{Z}=1.0,E_{V}=0.0,g_{z,F}=-2.7$. The FM phase is unstable
below the red line, while the CDW phase is unstable above the blue dotted
line. The dashed horizontal line represents $g_{z,H}=g_{z,F}$. The coexistence
between CDW and FM (the SVE+ phase, shaded green) occurs only for
$0>g_{z,H}>g_{z,F}$.

In what follows, we will fix the ratio of $g_{\mu,F}$ to $g_{\mu,H}$,
$\mu=xy,\ z$, and plot the phase diagrams with the axes being $g_{xy,H}$ and
$g_{z,H}$. This should merely be thought of as taking a certain two-
dimensional section of the full space of coupling constants, and does not
represent any physical assumption about the proportionality between
$g_{\mu,F}$ and $g_{\mu,H}$, say, as $B$ varies.

Fig. 4 shows the phase diagram for the case $g_{xy,F}=1.25g_{xy,H}$,
$g_{z,F}=1.25g_{z,H}$. As found earlier [49], there is a region of coexistence
between Kekulé and CAF order near the first-order line of the USR model. This
phase, which we call B/CAF, also has SVEY order. In addition, there is a
region of coexistence between CDW and FM order (the SVE+ phase) in the
neighborhood of the phase transition between the CDW and FM phases in the USR
case. The transitions bordering the coexistence regions of the phase diagram
are second-order, represented by dashed lines. As in Fig. 1, $E_{Z}=1$.

Figure 4: Phase diagram for
$E_{Z}=1.0,g_{z,F}=1.25g_{z,H},g_{xy,F}=1.25g_{xy,H},E_{V}=0.0$. The ratios
$g_{\mu,F}/g_{\mu,H}>1$ for both types of couplings. There are two different
coexistence phases: The B/CAF coexistence phase (shaded blue) was already
found previously [49], and also has SVEY order. The first-order transition
between the FM and CDW phases for USR interactions has broadened into the SVE+
phase (shaded green) with coexistence between CDW and FM order parameters, so
called because it has a spontaneous spin-valley entangled order parameter.

Fig. 5 shows the phase diagram for the case $g_{z,F}/g_{z,H}<1,\
g_{xy,F}/g_{xy,H}>1$. As mentioned above, this satisfies the condition for
coexistence between CAF and BO order (the B/CAF phase), but fails to meet the
condition for the existence of the SVE phase between the CDW and FM phases.
The B/CAF phase with its attendant SVEY order is still present, though reduced
in extent. Similarly, Fig. 6 has $g_{z,F}/g_{z,H}>1,\ g_{xy,F}/g_{xy,H}<1$,
allowing the SVE+ phase to exist but forbidding coexistence between BO and CAF
order. Finally, Fig. 7 has $g_{z,F}/g_{z,H}<1,\ g_{xy,F}/g_{xy,H}<1$,
disallowing any coexistence. The reason for the difference in detail between
the Kharitonov phase diagram Fig. 1 and Fig. 7 is that the boundaries between
phases sometimes depend on the Hartree coupling and sometimes on the Fock
coupling, as seen by the ground state energies in Eq. 29.

Figure 5: Phase diagram for
$E_{Z}=1.0,g_{z,F}=0.75g_{z,H},g_{xy,F}=1.25g_{xy,H},E_{V}=0.0$. The ratios of
couplings ensure $g_{xy,F}/g_{xy,H}>1$, but $g_{z,F}/g_{z,H}<1$. The B/CAF
coexistence region (shaded blue) still occurs but there is no coexistence near
the CDW/FM phase boundary, which remains first-order.  Figure 6: Phase diagram
for $E_{Z}=1.0,g_{z,F}=1.25g_{z,H},g_{xy,F}=0.75g_{xy,H},E_{V}=0.0$. Now the
conditions are no longer met for the B/CAF phase to occur. However, the SVE+
phase (shaded green) does occur near the CDW/FM phase boundary, which is now
split into two second-order lines.  Figure 7: Phase diagram for
$E_{Z}=1.0,g_{z,F}=0.75g_{z,H},g_{xy,F}=0.75g_{xy,H},E_{V}=0.0$. There is no
coexistence.

It is of experimental interest to ask how the system evolves with increasing
Zeeman coupling when all interaction parameters are fixed, and $E_{V}=0$. This
corresponds to applying a tilted field to the system [29], keeping the
perpendicular component of $B$ constant. The evolution of the order parameters
depends on the values chosen for $g_{z,H/F},g_{xy,H/F}$. Clearly, if the
system is in the FM phase already at $E_{Z}=0$, there will be no further
evolution with increasing $E_{Z}$. Similarly, if the system is in the CAF
phase at $E_{Z}=0$, it can only evolve into the FM with increasing $E_{Z}$.
Since these possibilities have been thoroughly explored in the past
literature, we will ignore them in favor of more interesting possibilities.

Figure 8: Evolution of order parameters vs $E_{Z}$ for
$g_{z,H}=-1.2,g_{z,F}=1.25g_{z,H},g_{xy,H}=-0.5,g_{xy,F}=1.25g_{xy,H},E_{V}=0.0$.
The ratios of the H and F couplings corresponds to Fig. 4. The system is in
CDW phase at $E_{Z}=0$. As $E_{Z}$ increases there is a second-order
transition to a phase where FM and CDW order coexist. This coexistence is
mediated by the presence of the SVE+ order parameter, which spontaneously
breaks the $U(1)_{s}$ and $U(1)_{v}$ symmetries. For large enough $E_{Z}$ the
system goes into FM phase. Figure 9: Evolution of order parameters vs $E_{Z}$
for
$g_{z,H}=-0.5,g_{z,F}=1.25g_{z,H},g_{xy,H}=-1.0,g_{xy,F}=1.25g_{xy,H},E_{V}=0.0$.The
ratios of the H and F couplings corresponds to the phase diagram of Fig. 4.
The system is in BO phase at $E_{Z}=0$. As $E_{Z}$ increases the system
undergoes a second-order transition to the intermediate B/CAF phase, where BO,
CAF, and SVEY order coexist. As $E_{Z}$ increases further, there is another
second-order transition to the pure CAF phase, and eventually a third second-
order transition to the FM phase at large $E_{Z}$.

Let us first consider the case $g_{z,F}/g_{z,H}>1,\ g_{xy,F}/g_{xy,H}>1$,
which corresponds to the phase diagram of Fig. 4. Fig. 8 shows the evolution
of order parameters with $E_{Z}$ when the system is in the CDW phase at
vanishing Zeeman coupling. In addition to the FM and CDW order parameters, we
also show the SVE+ order parameter. As $E_{Z}$ increases, the system undergoes
a phase transition from the CDW phase into the SVE+ phase (coexisting FM, CDW,
and SVE+ order), and then into the FM phase. Alternatively, as shown in Fig. 9
the system could start in the BO phase at $E_{Z}=0$. In this case the system
first goes into the B/CAF phase (which also has SVEY order), then into the
pure CAF phase, and finally into the FM phase at large $E_{Z}$.

Coming next to the case of $g_{z,F}/g_{z,H}<1,\ g_{xy,F}/g_{xy,H}>1$, which
corresponds to the phase diagram of Fig. 5, Fig. 10 shows the evolution of
order parameters as a function of $E_{Z}$ for $g_{z,H}=-0.2,g_{xy,H}=-1.2$.
The system starts in the BO phase at vanishing Zeeman coupling, undergoes a
second-order transition to the B/CAF coexistence phase, undergoes yet another
second-order transition to the pure CAF phase, and finally goes into the FM
phase at large $E_{Z}$.

Figure 10: Evolution of order parameters vs $E_{Z}$ for
$g_{z,H}=-0.2,g_{z,F}=0.75g_{z,H},g_{xy,H}=-1.2,g_{xy,F}=1.25g_{xy,H},E_{V}=0.0$.
The ratios of the H and F couplings correspond to the phase diagram of Fig. 5.
The system is in the BO phase at $E_{Z}=0$. As $E_{Z}$ increases there is a
second-order transition to the B/CAF phase (with BO, CAF, FM, and SVEY order).
For larger $E_{Z}$, there is another second-order transition to the pure CAF
phase. Finally, the system goes through another second-order transition into
the FM phase at very large $E_{Z}$.

Finally, we consider the case $g_{z,F}/g_{z,H}<1,\ g_{xy,F}/g_{xy,H}>1$, which
corresponds to the phase diagram of Fig. 6. Fig. 11 shows the evolution of the
order parameter as a function of $E_{Z}$ for $g_{z,H}=-1.5,\ g_{xy,H}=-0.5$.
The system starts in the CDW phase at vanishing $E_{Z}$, makes a second-order
phase transition to the SVE+ phase at intermediate $E_{Z}$, and finally goes
into the FM phase via another second-order phase transition at large $E_{Z}$.

Figure 11: Evolution of order parameters vs $E_{Z}$ for
$g_{z,H}=-1.5,g_{z,F}=1.25g_{z,H},g_{xy,H}=-0.5,g_{xy,F}=0.75g_{xy,H},E_{V}=0.0$.
The ratios of the H and F couplings correspond to the phase diagram of Fig. 6.
The system is in the CDW phase at $E_{Z}=0$. As $E_{Z}$ increases, the system
enters the SVE+ phase (coexisting CDW, FM, and SVE+ order) via a second-order
transition. For larger $E_{Z}$ there is another second-order transition to the
FM phase.

###  III.2 Nonvanishing Valley Zeeman Coupling and $g_{\mu,F}/g_{\mu,H}>0$

In this subsection we consider how the phase diagrams change when $E_{V}>0$.
We will still stay “close” to the USR model, assuming $g_{\mu,F}/g_{\mu,H}>0$.
The most obvious change is already present in the USR limit: The BO phase is
replaced by the B/CO phase, in which both bond order and CDW order coexist.
This is analogous to the replacement of the antiferromagnetic phase at
$E_{Z}=0$ by the CAF phase at arbitrarily small $E_{Z}$. Concurrently, the
first-order phase transition between the CDW and BO phases in the USR phase
diagram Fig. 1 is converted into a second-order transition. These are all
previously known results [56, 53].

The introduction of small $E_{V}$ does not modify the SVE+ phase (whenever it
occurs). While $E_{V}>0$ does not change the symmetries broken in the B/CAF
phase, it does change the phase quantitatively. Let us first look at Fig. 12,
which shows the order parameters at $E_{V}=0$ along a horizontal cut across
Fig. 4 at $g_{z,H}=-0.45$. Recall that the B/CAF phase, at $E_{V}=0$, has BO,
CAF, FM, and SVEY order. The SVEY order parameter changes discontinuously at
the first-order phase transition with the SVE+ phase. All these features
change when one adds a tiny $E_{V}$. In Fig. 13 we show the evolution of order
parameters along the same horizontal cut ($g_{zH}=-0.45$) when $E_{V}=0.01$.
It can be seen that the SVE+ phase does not change character qualitatively.
However, the B/CAF phase now changes significantly. Firstly, it acquires a
nonzero CDW order parameter, because the B/CO phase also has CDW order. More
importantly, it now has both SVEX and SVEY order, and undergoes a second-order
phase transition to the SVE+ phase. Thus, even a very tiny $E_{V}$ can affect
certain phases and phase transitions significantly.

Figure 12: Evolution of order parameters at $E_{V}=0$ as a function of
$g_{xy,H}$ for $g_{z,H}=-0.45,g_{z,F}=1.25g_{z,H},g_{xy,F}=1.25g_{xy,H}$. The
ratios of the H and F couplings correspond to the phase diagram of Fig. 4. It
can be seen that the system undergoes a first-order transition from the B/CAF
to the F/CDW phase. Figure 13: Evolution of order parameters as a function of
$g_{xy,H}$ for $g_{z,H}=-0.45,g_{z,F}=1.25g_{z,H},g_{xy,F}=1.25g_{xy,H}$ for
vanishingly small $E_{V}=0.01$. The ratios of the H and F couplings correspond
to the phase diagram of Fig. 4. The first-order transition between the B/CAF
and SVE+ phases has now become second-order.

Next we turn a moderate value of $E_{V}=0.25$. Fig. 14 shows the phase diagram
for the ratio of Hartree and Fock parts of the couplings being
$g_{z,F}/g_{z,H}=g_{xy,F}/g_{xy,H}=1.25$, the same as in Fig. 4. The topology
of the phase diagram is identical to that of Fig. 4, with the solitary change
that the first-order transition between the B/CAF and SVE+ phases has been
replaced by a second-order transition. The B/CAF region has shrinks, while the
SVE+ region expands.

Figure 14: For this figure the coupling ratios of $g_{xy}$ and $g_{z}$ are
same as in Fig. 4 and we have considered a small valley Zeeman $E_{V}=0.25$.By
comparing this figure with Fig. 4 one can see that for this choice of coupling
ratios, the presence of finite $E_{V}$ doesn’t qualitatively change the nature
of the phase diagram. Figure 15: For this figure the coupling ratios of
$g_{xy}$ and $g_{z}$ are same as in Fig. 5 and we have considered a valley
Zeeman field $E_{V}=0.5$.By comparing this figure with Fig. 5 one can see that
for this choice of coupling ratios, the presence of finite $E_{V}$ reduce the
coexistence phase area of B/CAF. And here CDW and CAF phases, CDW and FM
phases are separated by first order lines. Figure 16: Phase diagram at
$E_{V}=0.5$ for the coupling ratios $g_{z,F}/g_{z,H}=1.25,\
g_{xy,F}/g_{xy,H}=0.75$. Comparing this with the phase diagram for the same
ratios at $E_{V}=0$ (Fig. 6) one can see that for this choice of coupling
ratios, $E_{V}$ has suppressed the B/CAF phase completely and shrunk the SVE+
phase.

Fig. 15 shows the phase diagram at $E_{V}=0.5$ when the ratios are
$g_{z,F}=0.75g_{z,H};\ \ g_{xy,F}=1.25g_{xy,H}$. As in the case of the
previous figure, some of the phase boundaries move, but the topology remains
the same as at $E_{V}=0$. The same is true for $g_{z,F}/g_{z,H}=1.25;\ \
g_{xy,F}/g_{xy,H}=0.75$, shown in Fig. 16 and $g_{\mu,F}/g_{\mu,H}=0.75$ shown
in Fig. 17.

Figure 17: Phase diagram at $E_{V}=0.5$ for the coupling ratios
$g_{\mu,F}/g_{\mu,H}=0.75$. There is no qualitative change as compared to the
corresponding phase diagram at $E_{V}=0$ (Fig. 7)

More interesting phenomena occur if one increases $E_{V}$ to an even larger
value. Recall that for $g_{z,F}/g_{z,H}=g_{xy,F}/g_{xy,H}=0.75$, there is no
coexistence anywhere in the phase diagram (Fig. 7) at $E_{V}=0$, or at
$E_{V}=0.5$ (Fig. 17). The phase diagram for this ratio of the H and F parts
of the couplings, at $E_{V}=1.4$, is shown in Fig. 18. It can be seen that a
sufficiently large $E_{V}$ can induce coexistence, even when it does not occur
at $E_{V}=0$. A similar phenomenon is seen in previous work in BLG [53].

Figure 18: Phase diagram for $E_{Z}=1$, and $E_{V}=1.4$. The ratios of the
Hartree and Fock couplings are $g_{z,F}=0.75g_{z,H};\ \ g_{xy,F}=0.75g_{xy,H}$
as in Fig. 7 and Fig. 17. In contrast to those cases where there is no
coexistence anywhere in the phase diagram, we do obtain the coexistence phases
B/CAF (blue shaded region) and SVE+ (green shaded region).  Figure 19: The
coupling constants are fixed such that at $E_{V}=0$ there is no B/cAF phase.
We choose $E_{Z}=1,\ g_{z,H}=3.4,\ g_{xy,H}=-2.5,\ g_{xy,F}=-2.1$. Varying the
other two parameters, $g_{z,F},\ E_{V}$, we see that there is an intermediate
regime in both in which the B/CAF phase does appear.

Fig. 19 shows a slightly different way of looking at the occurrence of the
B/CAF phase as $E_{V}$ varies. We have fixed the Zeeman coupling at $E_{Z}=1$,
and the $xy$ interactions such that $g_{xy,F}/g_{xy,H}<1$, implying that
coexistence will not occur for $E_{V}=0$. The coupling $g_{z,H}$ is also
fixed. We show the range of parameters in $E_{V},\ g_{z,F}$ where the B/CAF
phase occurs. As can be seen, there is an intermediate range of $E_{V}$ and
$g_{z,F}$ where the B/CAF phase appears.

Figure 20: The parameters are $E_{Z}=1.0,g_{z,F}=0.8,g_{xy,F}=-0.6$. One can
see from this figure also that only when $(|g_{z,H}|>|g_{z,F}|)$ (for both
$g_{z,H},g_{z,F}>0$), their can be a coexistence phase of (FM+CDW) for some
range of $E_{V}$. The vertical dashed line denotes $g_{z,H}=g_{z,F}$.

Fig. 20 shows the fate of the SVE+ phase at large $E_{V}$. Once again, we
choose $g_{z,F}/g_{z,H}<1$, which implies that the SVE+ phase does not occur
for $E_{V}=0$. We see that for large enough $E_{V}$, the SVE+ phase is
stabilized.

It is also interesting to consider the evolution of the order parameters as
$E_{Z}$ increases in Fig. 18. If one starts in the B/CO phase, the expected
sequence of phases as $E_{Z}$ increases for small or vanishing $E_{V}$,
provided the B/CAF phase occurs, is B/CO$\ \rightarrow\ $B/CAF$\ \rightarrow\
$CAF$\ \rightarrow\ $FM. However, for $E_{Z}=1.4$ the sequence is different:
BO$\ \rightarrow\ $B/CAF$\ \rightarrow\ $SVE+$\ \rightarrow\ $FM. This is
shown in Fig. 21.

Figure 21: Evolution of order parameters vs $E_{Z}$ for
$g_{z,H}=0.9,g_{z,F}=0.75g_{z,H},g_{xy,H}=-0.8,g_{xy,F}=0.75g_{xy,H},E_{V}=1.4$.
The ratios of the H and F couplings correspond to the phase diagram of Fig.
18. The system is in the B/CO phase at $E_{Z}=0$. As $E_{Z}$ increases there
is a second-order transition to the B/CAF phase (with BO, CDW, CAF, FM, SVEX
and SVEY order). For larger $E_{Z}$, there is another second-order transition
to the SVE+ phase(with FM and CDW order). The system finally goes into the FM
phase at very large $E_{Z}$.

###  III.3 Hartree and Fock parts of $g_{z/xy}$ having opposite signs

It is conceivable that for strong LL-mixing, the renormalizations of the
interactions could be large enough to make the signs of the Hartree and Fock
parts of $g_{z/xy}$ opposite. For completeness we present some phase diagrams
for this type of situation in this subsection.

First we consider the case $g_{z,F}/g_{z,H}=-1$, but $g_{xy,F}/g_{xy,H}>1$.
The most obvious change is in the topology of the phase diagram. The CAF phase
is completely surrounded by other phases. In addition, there are two
coexistence phases. The blue shaded region is the B/CAF phase familiar from
the previous subsections. It has BO, CAF, and SVEY order. The brown shaded
region represents a new type of coexistence phase which does not occur when
$g_{\mu,F}/g_{\mu,H}>0$. This phase has FM order coexisting with SVEX/SVEY
order (all the ground states generated from SVEX by $U(1)_{s}\otimes U(1)_{v}$
are degenerate), with no other order parameters being present. This state
spontaneously breaks the $U(1)_{s}$ and $U(1)_{v}$ symmetries, but preserves
$U(1)_{sv}$. Fig. 23 shows the order parameters along a horizontal cut in the
phase diagram of Fig. 22 at $g_{z,H}=1$. We will call this the FSVE phase.

Figure 22: Phase diagram for
$E_{Z}=1.0,E_{V}=0.0,g_{z,F}=-1.0*g_{z,H},g_{xy,F}=1.25*g_{xy,H}$. The CAF now
occurs in a finite region of the parameter space. There are two coexistence
phases. The blue shaded region is the B/CAF phase with BO, CAF, and SVEY
order. The region shaded brown shows a new type of coexistence between FM and
SVEX/SVEY order.  Figure 23: Nonzero order parameters along a horizontal cut
of the phase diagram of Fig. 22 at $g_{z,H}=0.8$. All other coupling constants
are the same as in Fig. 22. One starts in the BO phase at large negative
$g_{xy,H}$. There is a second-order transition into the B/CAF phase with BO,
CAF, and SVEY order. Next, there is another second-order transition into the
FSVE phase showing the coexistence of FM and SVEX/SVEY order, the two being
degenerate. We have chosen to plot
SVE-$=\langle|\tau_{x}\sigma_{x}-\tau_{y}\sigma_{y}|\rangle/2$, which is
continuous across the transition. Finally, for larger $g_{xy,H}$ the system
goes into the FM phase.

When $E_{V}>0$, the FSVE phase acquires both SVEX and SVEY order in addition
to FM and CDW order (which is natural since $E_{V}>0$). Thus, the FSVE phase
remains different from the B/CAF phase (because there is no CAF order in the
FSVE phase) and the SVE+ phase, which requires a precise equality of the SVEX
and SVEY order parameters. More details are in the appendices.

Figure 24: Phase diagram for the ratios $g_{z,F}/g_{z,H}=1.25>0$ and
$g_{xy,F}/g_{xy,H}=-1<0$. The other couplings are $E_{Z}=1,\ E_{V}=0$. While
the topology has changed, the SVE+ phase that interpolates between the CDW and
FM phases occurs here as well.

As a second example, we show the case when $g_{z,F}/g_{z,H}=1.25>0$ and
$g_{xy,F}/g_{xy,H}=-1<0$, which is shown in Fig.24. As expected for negative
ratios between the H and F parts of the coupling, the phases are drastically
rearranged in the phase diagram. However, the nature of the phases remains the
same as in Section III.1. There is no B/CAF phase, but the SVE+ phase does
occur, interpolating between the CDW and FM phases.

A more complete set of figures for other cases when one or both of the H and F
coupling ratios are negative appears in the appendices.

##  IV Caveats, Conclusions, and Open Questions

Monolayer graphene at charge neutrality in the quantum Hall regime is an
example of quantum Hall ferromagnetism [9, 10, 11, 12]. While earlier
theoretical [31] and experimental work [29, 33] seemed to suggest a simple
canted antiferromagnetic phase at perpendicular magnetic field, recent
scanning tunneling observations [37, 38, 39] clearly show Kekulé and even
charge density wave order.

A major assumption underlying most earlier theoretical work [18, 20, 21, 31]
is that the residual interactions (other than Coulomb), being descended from
lattice-scale couplings, should be ultra-short-range on the scale of the
magnetic length $\ell$. Recently, based on renormalization group ideas, it was
argued [49] that the effective low-energy interactions will naturally acquire
the length scale $\ell$ even if they were ultra-short-range microscopically.
It should be noted that this argument does not depend on whether one considers
the $N=0$ manifold (the ZLLs) or some other manifold of Landau levels of
graphene. For such generic interactions, coexistence between CAF and Kekulé
order in the physically relevant region of the coupling constant space was
shown to occur [49] in a robust regime of couplings.

Our goal in this paper was to present a more complete phase diagram of MLG at
charge neutrality (and more generally, in any half-filled manifold of Landau
levels). As in previous theoretical work, we assume a clean system with four-
Fermi interactions only. Momentum conservation constrains the residual four-
Fermi interactions to have a $U(1)$ valley symmetry (reduced to a $Z_{3}$ upon
including higher-Fermi interactions, which we ignore). This still leaves two
coupling functions $v_{z}({\mathbf{q}}),v_{xy}({\mathbf{q}})$ to be
determined. In the Hartree-Fock approximation, restricting to ground states
respecting translation invariance up to an intervalley coherence, the
energetics can be captured by a Hartree coupling and a Fock coupling for each
of the two types. After the inclusion of the Zeeman and valley Zeeman
couplings $E_{Z},\ E_{V}$, one has six coupling constants to deal with.

Our main physically relevant finding is that coexistence between CAF and B/CO
order (the bond order always coexists with CDW order when $E_{V}>0$) occurs in
a large region of the coupling constant space when $g_{xy,H},g_{xy,F}<0$.
Depending on the particular values of the interaction couplings, the valley
Zeeman coupling may either enhance or suppress coexistence. In fact, in a
certain range of couplings, a re-entrant transition from a pure B/CO through a
coexistence phase back to a pure B/CO is possible upon increasing $E_{V}$. The
fact that coexistence is generic is consistent with the ubiquity of B/CO order
seen in STM observations [37, 38, 39].

We found two other coexistence phases in regions of the phase diagram which
may not be directly relevant to graphene. The first, called the SVE+ phase,
interpolates between the CDW and FM phases, and also has a spin-valley
entangled SVE+ order parameter. The corresponding phase in BLG was found in
previous work [53]. The second one, called the FSVE phase, occurs when
$g_{z,F}$ is allowed to be of the opposite sign as $g_{z,H}$, which could
conceivably occur for strong Landau-level mixing. For generic $E_{V}>0$, the
FSVE phase has FM, CDW, SVEX and SVEY order, but no CAF or bond order.

Our results are complete given our assumptions, but they come with caveats. We
have ignored spontaneous translation symmetry breaking beyond that required by
intervalley coherence, static disorder, and quantum/thermal fluctuations. Let
us consider each in turn. Intervalley coherence is allowed in our approach,
and is a signal of minimal translation symmetry breaking because it implies a
new reciprocal lattice vector ${\mathbf{K}}-{\mathbf{K}}^{\prime}$. In
combination with the reduction of the $U_{v}(1)$ symmetry to $Z_{3}$ upon
including higher-Fermi interactions, this leads to the Kekulé order, which
does break lattice translation symmetry with a three-fold enlargement of the
unit cell. Thus, this type of translation symmetry breaking in implicitly
included in our approach. Since the STM experiments see only this minimal type
of translation symmetry breaking, we believe our assumption of translation
invariance up to an intervalley coherence is not a serious limitation. Next,
static disorder is present in all samples. When B/CO order is present, it is
expected to be pinned by the local value of disorder. In fact, since the B/CO
order couples to static disorder while the CAF order does not, one expects the
region of the coupling constant space where B/CO order is present to increase
as disorder increases. Otherwise, static disorder is expected to have a
quantitative effect on the transport gaps, but leave the nature of the state
unaffected.

Now we turn to quantum fluctuations. The Hartree-Fock approximation has an
excellent track record in describing quantum Hall ferromagnets at zero
temperature [10, 9, 11, 12]. Near a second-order phase boundary, HF will
always predict a mean-field transition. Quantum fluctuations will shift the
phase boundary, and correct the critical behavior to the appropriate
universality class; for example, the transition from the B/CO to the
coexistence phase with the CAF order should have the universality class of the
three-dimensional $XY$ model. If a phase occurs in a very narrow sliver of
coupling constant space, one might envisage quantum fluctuations making it
disappear. However, since all the phases we find are robust, occurring over
substantial ranges of coupling constants, we expect quantum fluctuations to
alter the phase boundaries in detail, but not affect the phase diagram
qualitatively. An important potential exception is the first-order transition
seen in many parameter regimes between the FM and CDW phases. When $E_{Z},\
E_{V}>0$, neither of these phases breaks any symmetry of the Hamiltonian. One
cannot rule out a first-order transition without a change in symmetry between
the two phases: An example is the liquid-gas transition. However, a phase
transition between the FM and CDW phases is not necessary, since both have the
same symmetry. Quantum fluctuations may destroy the first-order line in favor
of a smooth crossover between the CDW and FM. Consider parameter regimes when
the SVE+ phase intervenes between the CDW and FM mean-field phases. Since the
SVE+ phase breaks $U(1)_{s}$ and $U(1)_{v}$ spontaneously, a second-order
phase transition is allowed between it and the FM or CDW regions. Quantum
fluctuations may change the topology of the phase diagram to make the SVE+
phase an island in the middle of the crossover between the CDW-dominated and
FM-dominated regions.

Finally, we turn to $T>0$ and thermal fluctuations. Recall that spontaneous
intervalley coherence is subject to a $Z_{3}$ symmetry, and hence does not
result in a Goldstone mode. Long-range B/CO order is expected to be survive to
a critical temperature $T_{c}>0$. The $U_{spin}(1)$ symmetry is an almost
exact symmetry, being broken only by the tiny spin-orbit coupling ($\approx
10\mu eV$) in graphene [59, 57, 58, 60]. Thus, setting spin-orbit coupling to
zero, any spontaneous breaking of the $U(1)_{s}$ symmetry leads to a gapless
Goldstone mode at $T=0$. For $T>0$ the system is in the universality class of
the two-dimensional $XY$-model. Long-range order is absent at any nonzero $T$,
and there should be a Berezinskii-Kosterlitz-Thouless (BKT) transition [63,
64] at $T_{KT}$, below which there is power-law order. $T_{KT}$ should vanish
as the CAF order parameter vanishes, because the stiffness vanishes as well.
Thus, there is an intermediate temperature regime in the B/CAF phase which is
above $T_{KT}$ but below the $T_{c}$ of the $Z_{3}$ bond order.

Let us turn to experimental signatures of coexistence. STM experiments can
directly measure the B/CO order; thus the key question is how to detect CAF
order. Any order parameter that spontaneously breaks $U(1)_{s}$ will lead to a
gapless Goldstone mode. Thus, the CAF phase and all three coexistence phases
we have found would support a Goldstone mode at $T=0$. Magnon transmission
experiments [33, 35, 36] reveal the presence of magnetic excitations, but
because the magnons are created in ferromagnetic regions with a gap of
$E_{Z}$, such experiments are unable to reveal whether the magnetic
excitations at $\nu=0$ are gapless. In Bernal-stacked bilayer graphene, where
a CAF state is also expected to be present, a very recent experiment has
confirmed the presence of the gapless Goldstone mode [65]. If such an
experiment can be carried out for MLG it would be direct confirmation of the
spontaneous breaking of spin-rotation symmetry $U(1)_{s}$. More broadly, in
the context of transport, a measurement of the bulk thermal conductivity below
$T_{KT}$ should reveal the Goldstone mode. Additionally, the BKT transition
itself should have a signature in $R_{xx}$ [66].

A more indirect way to probe the CAF order parameter is to examine the
detailed structure of spin/valley skyrmions, which can be induced by external
charges. A thorough analysis of skyrmions in the ultra-short-range model of
charge-neutral graphene was carried out very recently [56]. Using this
framework, an analysis of the B/CO texture near a charge defect shows [38]
that it is consistent with theory, assuming that the true ground state is pure
B/CO. A phase with coexistence between B/CO and CAF order will have skyrmions
that differ in detail from those of the pure B/CO phase.

There are two broad open questions. Firstly, given a microscopic model at some
intermediate energy scale much larger than $\hbar\omega_{c}$, how does one
reliably deduce the effective coupling functions in the manifold of the $n=0$
Landau levels? Kharitonov [31], following earlier RG treatments [67, 68],
carried out just such a procedure, under the assumption that the couplings
(other than Coulomb) remained ultra-short-range under RG, implying a finite
number of couplings to renormalize. Based on a general fermionic RG procedure
[69] which includes all low-energy interactions, there has been quite a bit of
previous work attempting to integrate out higher Landau levels perturbatively
[44, 45, 46, 47, 48]. While these works restricted themselves to the Coulomb
interaction, it should be straightforward to extend them to include all
symmetry-allowed interactions.

The second broad question is complementary to the first. Given an experimental
sample, is there a complete set of measurements that could determine the
couplings $g_{z,H},g_{z,F},g_{xy,H},g_{xy,F}$? Given that $E_{Z}$ is
determined by the total field, and $E_{V}$ can be deduced from zero-magnetic-
field gap measurements at charge neutrality, this would fully determine the
effective theory at the mean-field level. One way to approach this is via a
detailed investigation of skyrmions [56]. As long as the size of the skyrmions
is much larger than the magnetic length, a nonlinear sigma model approach is
capable of capturing their structure and energetics. The parameters that enter
the nonlinear sigma model are exactly those that enter the mean-field theory,
with the exception of the stiffness, which is determined by the dominant
Coulomb interaction.

Last, but not least, let us briefly consider fractionally filled states in the
$n=0$ manifold of Landau levels in graphene. For the case of $SU(4)$ Coulomb
interactions plus ultra short range residual interactions, it is possible to
construct variational states with integer and/or fractional fillings in the
different flavors and compute their energies [70, 71]. Determining whether
this construction can be extended to generic residual interactions of
arbitrary range is an important open question.

We look forward to addressing these and other questions in the near future.

###### Acknowledgements.

SJD would like to acknowledge the Infosys funding for final year students. He
also wants to acknowledge ICTS for its hospitality and kind support towards
academic collaboration. AD was supported by the German-Israeli Foundation
Grant No. I-1505-303.10/2019, DFG MI 658/10-2, DFG RO 2247/11-1, DFG EG
96/13-1, and CRC 183 (project C01). AD also thanks the Israel planning and
budgeting committee (PBC) and the Weizmann Institute of Science, the Dean of
Faculty fellowship, and the Koshland Foundation for financial support. S.R.
and G.M. would like to thank the VAJRA scheme of SERB, India for its support.
R.K.K. was supported in part by NSF DMR-2026947. G.M. would like to
acknowledge partial support from the US-Israel Binational Science Foundation
(grant no. 2016130). G.M. and R.K.K. are grateful for the wonderful
environment at the Aspen Center for Physics (NSF grant PHY-1607611).

##  Appendix A Different Phases and Hessian

For charge-neutral graphene, assuming that interactions are ultra-short-range
(USR), there are four phases, namely, the ferromagnet (FM), the canted
antiferromagnet (CAF), the bond-ordered phase (BO), and the charge density
wave (CDW) phase. We call these phases “simple”, because they can all be
described by a single nontrivial angle which is known analytically in terms of
the couplings. As shown in the main text, when one removes the USR restriction
on the interactions, other coexistence phases become possible. Among them is a
spin-valley entangled phase SVE+, which can also be described by a single
angle, and is also “simple”.

A primary tool in our investigation of the phase diagram is a study of the
stability of a given ground state. Recall that each candidate state is
described by filling in two linear combinations of the four spin-valley
degenerate states at each guiding center. We reproduce the equations from the
main text here for convenience, [55, 50, 51, 52]

$\displaystyle|f_{1}\rangle=$
$\displaystyle\cos{\frac{\alpha_{1}}{2}}|{\mathbf{n}}\rangle\otimes|{\mathbf{s}}\rangle+e^{i\beta_{1}}\sin{\frac{\alpha_{1}}{2}}|-{\mathbf{n}}\rangle\otimes|-{\mathbf{s}}\rangle$
(32) $\displaystyle|f_{2}\rangle=$
$\displaystyle\cos{\frac{\alpha_{2}}{2}}|{\mathbf{n}}\rangle\otimes|-{\mathbf{s}}\rangle+e^{i\beta_{2}}\sin{\frac{\alpha_{2}}{2}}|-{\mathbf{n}}\rangle\otimes|{\mathbf{s}}\rangle$
(33)

where
${\mathbf{n}}=\sin{\theta_{p}}\cos{\phi_{p}}{\hat{e}}_{x}+\sin{\theta_{p}}\sin{\phi_{p}}{\hat{e}}_{y}+\cos{\theta_{p}}{\hat{e}}_{z}$,
and
${\mathbf{s}}=\sin{\theta_{s}}\cos{\phi_{s}}{\hat{e}}_{x}+\sin{\theta_{s}}\sin{\phi_{s}}{\hat{e}}_{y}+\cos{\theta_{s}}{\hat{e}}_{z}$
are unit vectors indicating the directions of the state on the valley and spin
Bloch spheres, respectively. The spinors $|{\mathbf{n}}\rangle$ and
$|{\mathbf{s}}\rangle$ are defined in the standard way

$|{\mathbf{n}}\rangle=\left(\begin{array}[]{c}\cos{\frac{\theta_{p}}{2}}\\\
e^{i\phi_{p}}\sin{\frac{\theta_{p}}{2}}\end{array}\right);\\\
|{\mathbf{s}}\rangle=\left(\begin{array}[]{c}\cos{\frac{\theta_{s}}{2}}\\\
e^{i\phi_{s}}\sin{\frac{\theta_{s}}{2}}\end{array}\right)$ (34)

As shown in the main text, the $U(1)_{s}$ symmetry allows us to set
$\phi_{s}=0$, and the $U(1)_{v}$ symmetry allows us to set $\phi_{p}=0$. The
$SU(2)_{s}$ symmetry of the interactions forces the HF energy to depend only
on $\beta_{1}+\beta_{2}$. The energy depends on $\beta_{1}+\beta_{2}$ only via
$\cos(\beta_{1}+\beta_{2})$, which appears linearly in the energy. Therefore,
we can restrict consideration to the two discrete possibilities
$\beta_{1}+\beta_{2}=0,\ \pi$ mod $2\pi$. The HF energy of a given state thus
depends only on the four continuously varying angles
$\alpha_{1},\alpha_{2},\theta_{p},\theta_{s}$.

We now indicate how to determine the values of the four angles for the
“simple” states shown in Section II.2. We begin by examining the expressions
for the order parameters as functions of the four angles:

$\displaystyle FM$
$\displaystyle=\frac{\langle\sigma_{z}\rangle}{2}=\frac{(\cos[\alpha_{1}]-\cos[\alpha_{2}])\cos[\theta_{s}]}{2},$
$\displaystyle CAF$
$\displaystyle=\frac{\langle\tau_{z}\sigma_{x}\rangle}{2}=\frac{1}{2}\Bigg{[}\cos[\theta_{s}]\Big{(}\cos[\beta_{1}]\sin[\alpha_{1}]$
$\displaystyle+\cos[\beta_{2}]\sin[\alpha_{2}]\Big{)}\sin[\theta_{p}]\Bigg{]},$
$\displaystyle BO$
$\displaystyle=\frac{\langle\tau_{x}\rangle}{2}=\frac{(\cos[\alpha_{1}]+\cos[\alpha_{2}])\sin[\theta_{p}]}{2},$
$\displaystyle CDW$
$\displaystyle=\frac{\langle\tau_{z}\rangle}{2}=\frac{(\cos[\alpha_{1}]+\cos[\alpha_{2}])\cos[\theta_{p}]}{2}$
(35)

First we focus on the CAF phase. In this phase the BO and CDW order parameters
should vanish identically, implying that $\alpha_{2}=\pi-\alpha_{1}$.
Examining the FM and CAF order-parameters, we see that we can choose
$\theta_{s}=0$ and $\theta_{p}=\pi/2$ because $\cos{\theta_{s}}$ and
$\cos{\theta_{p}}$ appear as a overall normalization factors which can be set
to one. After imposing these constraints on the angles, the HF energy will
depend only upon the angles $\alpha_{1}$, $\beta_{1}$, $\beta_{2}$, and can be
expressed as

$\displaystyle E_{HF}=$
$\displaystyle\frac{1}{4}\Big{(}-8E_{Z}\cos(\alpha_{1})$
$\displaystyle+2\sin^{2}(\alpha_{1})(g_{xy,F}-g_{z,F})\cos(\beta_{1}+\beta_{2})$
$\displaystyle-\cos(2\alpha_{1})(3g_{xy,F}+g_{z,F})$
$\displaystyle-5g_{xy,F}-3g_{z,F}\Big{)}.$ (36)

The CAF phase occurs for $g_{xy,F}<0$, $g_{z,F}>0$. This restricts
$\cos(\beta_{1}+\beta_{2})=1$, allowing us to choose
$\beta_{1}=\beta_{2}=\pi$. The CAF and FM order-parameters are
$FM=\cos[\alpha_{1}]$ , $CAF=\sin[\alpha_{1}]$. Here $\alpha_{1}$ is the
single non-trivial angle which varies through the CAF phase. It has the
functional dependence
$\alpha_{1}=\cos^{-1}\left[\frac{-E_{Z}}{2g_{xy,F}}\right]$, which we found by
minimizing Eq. 36. Clearly when $E_{Z}>2|g_{xy,F}|$, the angle $\alpha_{1}$
will be fixed to zero and this corresponds to the FM phase. The
parameterization of the FM phase is subsumed in the above.

Now we turn our attention to the B/CO phase which generally occurs for
$E_{V}>0$. Since the CDW and BO phases are restricted versions of the B/CO
phase, their parameterizations are subsumed in that of the B/CO phase. In the
B/CO phase the FM and CAF order-parameters should vanish identically, allowing
us to choose $\alpha_{1}=\alpha_{2}=0$. The B/CO phase is a singlet, which
means that the direction of ${\mathbf{s}}$ can be chosen arbitrarily, allowing
us to fix $\theta_{s}=0$. With this choice of angles the BO and CDW order-
parameters are $BO=\sin[\theta_{p}]$, $CDW=\cos[\theta_{p}]$. Thus, in this
phase the angle $\theta_{p}$ is the non-trivial angle. To find its functional
dependence, we examine the HF energy with the constraints on angles
$\alpha_{1}=\alpha_{2}=\theta_{s}=0$, which is

$\displaystyle E_{HF}=$
$\displaystyle-2E_{V}\cos(\theta_{p})+\frac{1}{2}g_{V}\cos(2\theta_{p})$
$\displaystyle-\frac{g_{xy,F}}{2}+g_{xy,H}-\frac{g_{z,F}}{2}+g_{z,H},$ (37)

where $g_{V}=2g_{z,H}-g_{z,F}-2g_{xy,H}+g_{xy,F}$. Minimizing Eq. 37 leads to
the functional dependence
$\theta_{p}=\cos^{-1}\left[\frac{E_{V}}{g_{V}}\right]$ for the B/CO phase. For
$E_{V}>g_{V}$, the angle $\theta_{p}$ will be fixed to zero, which corresponds
to the pure CDW phase. For $E_{V}=0$, the angle $\theta_{p}=\pi/2$, which
corresponds to the BO phase.

Having described how to fix the angles for the simple states, we now turn to
the instabilities of these states.

Any HF state must be at least a local extremum. Thus, all the first
derivatives of the energy $E_{HF}$ with respect to
$\alpha_{1},\alpha_{2},\theta_{p},\theta_{s}$ must vanish. To look for
instabilities we need to compute the second derivatives of $E_{HF}$ with
respect to the four angles (the Hessian matrix).

${\cal
E}(\alpha_{1},\alpha_{2},\theta_{p},\theta_{s})=\begin{pmatrix}\frac{{\partial^{2}{E_{\text{HF}}}}}{\partial^{2}{\alpha_{1}}}&\frac{{\partial^{2}{E_{\text{HF}}}}}{\partial{\alpha_{1}}\partial{\alpha_{2}}}&\frac{{\partial^{2}{E_{\text{HF}}}}}{\partial{\alpha_{1}}\partial{\theta_{p}}}&\frac{{\partial^{2}{E_{\text{HF}}}}}{\partial{\alpha_{1}}\partial{\theta_{s}}}\\\
\frac{{\partial^{2}{E_{\text{HF}}}}}{\partial{\alpha_{2}}\partial{\alpha_{1}}}&\frac{{\partial^{2}{E_{\text{HF}}}}}{\partial^{2}{\alpha_{2}}}&\frac{{\partial^{2}{E_{\text{HF}}}}}{\partial{\alpha_{2}}\partial{\theta_{p}}}&\frac{{\partial^{2}{E_{\text{HF}}}}}{\partial{\alpha_{2}}\partial{\theta_{s}}}\\\
\frac{{\partial^{2}{E_{\text{HF}}}}}{\partial{\theta_{p}}\partial{\alpha_{1}}}&\frac{{\partial^{2}{E_{\text{HF}}}}}{\partial{\theta_{p}}\partial{\alpha_{2}}}&\frac{{\partial^{2}{E_{\text{HF}}}}}{\partial^{2}{\theta_{p}}}&\frac{{\partial^{2}{E_{\text{HF}}}}}{\partial{\theta_{p}}\partial{\theta_{s}}}\\\
\frac{{\partial^{2}{E_{\text{HF}}}}}{\partial{\theta_{s}}\partial{\alpha_{1}}}&\frac{{\partial^{2}{E_{\text{HF}}}}}{\partial{\theta_{s}}\partial{\alpha_{2}}}&\frac{{\partial^{2}{E_{\text{HF}}}}}{\partial{\theta_{s}}\partial{\theta_{p}}}&\frac{{\partial^{2}{E_{\text{HF}}}}}{\partial^{2}{\theta_{s}}}\end{pmatrix}.$
(38)

The eigenvalues of the Hessian determine the stability of the given state. An
eigenvalue crossing zero signals an instability of the given state. For
“simple” states, one can analytically obtain the Hessian matrix and its
eigenvalues. We will use the Hessian eigenvalues of “simple” states to map out
the gross features of the phase diagram. Finer details of the phase diagrams
are obtained by self-consistent iterative HF.

In the following sections we will present explicit expressions for the Hessian
and its eigenvalues and attendant instabilities in the CAF (Appendix B), the
B/CO (Appendix C), the FM (Appendix D), the CDW (Appendix E), and the SVE
(Appendix F) phases. In Appendix G we present phase diagrams when the Hartree
and Fock parts of either/both of the couplings have opposite signs, which may
be relevant for very strong Landau-level mixing.

##  Appendix B CAF phase

The CAF phase occurs only for $g_{xy,F}<0$ and $E_{Z}<2|g_{xy,F}|$. It is
described by the following values of the angles

$\displaystyle\alpha_{1}=$
$\displaystyle\cos^{-1}\left[\frac{-E_{Z}}{2g_{xy,F}}\right]=\pi-\alpha_{2}$
(39) $\displaystyle\theta_{p}=$ $\displaystyle\frac{\pi}{2},\;\theta_{s}=0$
(40) $\displaystyle\beta_{1}=$ $\displaystyle\beta_{2}=\pi$ (41)

The energy of the CAF phase is

$E_{\text{CAF}}=\frac{E_{Z}^{2}}{2g_{xy,F}}-g_{z,F}.$ (42)

The Hessian matrix for the CAF state has the following block diagonal form:

${\cal E}_{\text{CAF}}=\left(\begin{array}[]{cc}A_{3\times 3}&0_{3\times 1}\\\
0_{1\times 3}&-\frac{E_{Z}^{2}}{g_{xy,F}}\end{array}\right)_{4\times 4},$ (43)

where

$\displaystyle A=$
$\displaystyle\left(\begin{array}[]{ccc}\frac{1}{4}\left(\frac{E_{Z}^{2}(g_{xy,F}-g_{xy,H})}{g_{xy,F}^{2}}-6g_{xy,F}+4g_{xy,H}+2g_{z,F}\right)&\frac{1}{2}(g_{xy,F}+2g_{xy,H}+g_{z,F})-\frac{E_{Z}^{2}(g_{xy,F}+g_{xy,H})}{4g_{xy,F}^{2}}&-E_{V}\sqrt{1-\frac{E_{Z}^{2}}{4g_{xy,F}^{2}}}\\\
\frac{1}{2}(g_{xy,F}+2g_{xy,H}+g_{z,F})-\frac{E_{Z}^{2}(g_{xy,F}+g_{xy,H})}{4g_{xy,F}^{2}}&\frac{1}{4}\left(\frac{E_{Z}^{2}(g_{xy,F}-g_{xy,H})}{g_{xy,F}^{2}}-6g_{xy,F}+4g_{xy,H}+2g_{z,F}\right)&-E_{V}\sqrt{1-\frac{E_{Z}^{2}}{4g_{xy,F}^{2}}}\\\
-E_{V}\sqrt{1-\frac{E_{Z}^{2}}{4g_{xy,F}^{2}}}&-E_{V}\sqrt{1-\frac{E_{Z}^{2}}{4g_{xy,F}^{2}}}&\frac{\left(E_{Z}^{2}-4g_{xy,F}^{2}\right)(g_{xy,F}-g_{z,F})}{2g_{xy,F}^{2}}\\\
\end{array}\right).$ (47)

Thus we find the instability equations of the CAF phase are

$\displaystyle g_{xy,F}$ $\displaystyle=-\frac{E_{Z}}{2}$ (48a) $\displaystyle
E_{V}$
$\displaystyle=\frac{\sqrt{g_{xy,F}-g_{z,F}}\sqrt{\frac{E_{Z}^{2}g_{xy,H}}{2}+g_{xy,F}^{3}-g_{xy,F}^{2}(2g_{xy,H}+g_{z,F})}}{g_{xy,F}}.$
(48b)

The first instability corresponds to the second-order line between the FM and
CAF phases, while the second corresponds to the instability in the B/CAF
phase.

##  Appendix C BO phase

The angles $\beta_{1}$ and $\beta_{2}$ never appear in the expression for the
HF energy of the BO state. Thus, $\beta_{1}$ and $\beta_{2}$ are undetermined
at this point. The other angles are

$\displaystyle\alpha_{1}=$ $\displaystyle\alpha_{2}=\theta_{s}=0$ (49)
$\displaystyle\theta_{p}=$
$\displaystyle\cos^{-1}\left[\frac{E_{V}}{g_{V}}\right],$ (50)

with

$g_{V}=2g_{z,H}-g_{z,F}-2g_{xy,H}+g_{xy,F}.$ (51)

The energy in the BO phase is

$E_{\text{BO}}=-\frac{E_{V}^{2}}{g_{V}}-g_{xy,F}+2g_{xy,H}.$ (52)

The Hessian matrix and its eigenvalues do depend on $\beta_{1}+\beta_{2}$,
with the eigenvalues depending on $\cos[2(\beta_{1}+\beta_{2})]$. The most
severe constraint on the region of stability occurs when
$\cos[2(\beta_{1}+\beta_{2})]=1$. Below, we choose $\beta_{1}=\beta_{2}=0$,
and obtain the form of the Hessian

${\cal E}_{\text{BO}}=\left(\begin{array}[]{cc}\frac{A_{2\times
2}}{2}&0_{2\times 2}\\\ 0_{2\times 2}&B_{2\times 2}\end{array}\right)_{4\times
4}$ (53)

where

$A=\left(\begin{array}[]{cc}\frac{E_{V}^{2}(g_{z,F}-g_{xy,F})}{(g_{xy,F}-2g_{xy,H}-g_{z,F}+2g_{z,H})^{2}}+2E_{Z}+g_{xy,F}-4g_{xy,H}-g_{z,F}&(g_{xy,F}-g_{z,F})\left(1-\frac{E_{V}^{2}}{(g_{xy,F}-2g_{xy,H}-g_{z,F}+2g_{z,H})^{2}}\right)\\\
(g_{xy,F}-g_{z,F})\left(1-\frac{E_{V}^{2}}{(g_{xy,F}-2g_{xy,H}-g_{z,F}+2g_{z,H})^{2}}\right)&\frac{E_{V}^{2}(g_{z,F}-g_{xy,F})}{(g_{xy,F}-2g_{xy,H}-g_{z,F}+2g_{z,H})^{2}}-2E_{Z}+g_{xy,F}-4g_{xy,H}-g_{z,F}\\\
\end{array}\right),$ (54)

and

$B=\left(\begin{array}[]{cc}-\frac{2E_{V}^{2}}{g_{xy,F}-2g_{xy,H}-g_{z,F}+2g_{z,H}}+2g_{xy,F}-4g_{xy,H}-2g_{z,F}+4g_{z,H}&0\\\
0&0\\\ \end{array}\right).$ (55)

The instability lines in the BO phase are

$\displaystyle g_{z,H}=$
$\displaystyle\frac{1}{2}(E_{V}-g_{xy,F}+2g_{xy,H}+g_{z,F})$ (56a)
$\displaystyle E_{V}=$
$\displaystyle\sqrt{\frac{\left(E_{Z}^{2}-2g_{xy,H}(-g_{xy,F}+2g_{xy,H}+g_{z,F})\right)(g_{xy,F}-2g_{xy,H}-g_{z,F}+2g_{z,H})^{2}}{2g_{xy,H}(g_{xy,F}-g_{z,F})}}.$
(56b)

These correspond to the instability of the CDW phase and the instability of
the B/CAF phase, respectively.

##  Appendix D FM phase

The FM phase is described by the angles

$\displaystyle\alpha_{1}$ $\displaystyle=0,$ $\displaystyle\ \alpha_{2}=\pi$
(57a) $\displaystyle\theta_{p}=$ $\displaystyle\frac{\pi}{2},\ \theta_{s}=0$
(57b) $\displaystyle\beta_{1}=$ $\displaystyle\beta_{2}=\pi$ (57c)

The ground state energy of the FM phase is

$E_{\text{FM}}=-2(E_{Z}+g_{xy,F})-g_{z,F}.$ (58)

The form of the hessian matrix of the FM phase is

${\cal E}_{\text{FM}}=\left(\begin{array}[]{cc}\frac{A_{2\times
2}}{2}&0_{2\times 2}\\\ 0_{2\times 2}&B_{2\times 2}\end{array}\right)_{4\times
4},$ (59)

where

$A=\left(\begin{array}[]{cc}2E_{Z}+3g_{xy,F}+g_{z,F}&g_{z,F}-g_{xy,F}\\\
g_{z,F}-g_{xy,F}&2E_{Z}+3g_{xy,F}+g_{z,F}\\\ \end{array}\right),$ (60)

and

$B=\left(\begin{array}[]{cc}0&0\\\ 0&2E_{Z}\\\ \end{array}\right).$ (61)

The instability lines in the FM phase are

$\displaystyle g_{xy,F}=$ $\displaystyle-\frac{E_{Z}}{2}$ (62a) $\displaystyle
g_{z,F}=$ $\displaystyle-E_{Z}-g_{xy,F}.$ (62b)

The first is the second-order line between the FM and CAF phases, while the
second is the instability towards the CDW phase. Note that this is not the
actual position of the dominant instability, which should depend on $E_{V}$ as
well. The reason is that $E_{V}$ does not appear in the Hessian matrix for the
FM phase. Note that one eigenvalue is always zero, which means one direction
in angle space is flat. While we believe that some higher-order derivative of
the ground state energy must reveal the instability towards the SVE+ phase, we
have not pursued this issue because we can find the instability easily from
the SVE+ side.

##  Appendix E CDW phase

The CDW phase is described by the angles

$\alpha_{1}=\alpha_{2}=\theta_{p}=\theta_{s}=0.$ (63)

The angles $\beta_{1},\ \beta_{2}$ never appear in the expression of the
ground state energy, which is

$E_{\text{CDW}}=-2E_{V}-g_{z,F}+2g_{z,H}.$ (64)

The angles $\beta_{1}$ and $\beta_{2}$ also do not appear in the Hessian
matrix, which has the following form

$\displaystyle{\cal E}_{\text{CDW}}=$
$\displaystyle\left(\begin{array}[]{cccc}E_{V}+E_{Z}-g_{xy,F}+g_{z,F}-2g_{z,H}&0&0&0\\\
0&E_{V}-E_{Z}-g_{xy,F}+g_{z,F}-2g_{z,H}&0&0\\\
0&0&2(E_{V}-g_{xy,F}+2g_{xy,H}+g_{z,F}-2g_{z,H})&0\\\ 0&0&0&0\\\
\end{array}\right).$ (69)

Note that one eigenvalue is always zero. This does not indicate instability
but rather the fact that one of the four continuously varying angles does not
appear in the Hessian. The lines of instability of the CDW phase are

$\displaystyle g_{z,H}=$
$\displaystyle\frac{1}{2}(E_{V}-g_{xy,F}+2g_{xy,H}+g_{z,F})$ (70a)
$\displaystyle g_{z,H}=$
$\displaystyle\frac{1}{2}(E_{V}-E_{Z}-g_{xy,F}+g_{z,F}).$ (70b)

The first is the instability toward the B/CO phase, while the second is the
instability towards the SVE+ phase.

##  Appendix F SVE+ phase

The SVE+ phase occurs near the boundary of the FM and CDW phases. In addition
to the FM and CDW order parameters, this phase also has a nonzero expectation
value of

$SVE+=\langle\left[\tau_{x}\sigma_{x}+\tau_{y}\sigma_{y}\right]/2\rangle$ (71)

Ordering the rows and columns as
${\mathbf{K}}\uparrow,~{}{\mathbf{K}}\downarrow,~{}{\mathbf{K}}^{\prime}\uparrow,~{}{\mathbf{K}}^{\prime}\downarrow$,
the projector matrix has the form

$\Delta_{SVE+}\left(\xi\right)=\left(\begin{array}[]{cccc}1&0&0&0\\\
0&\cos\left(\frac{\xi}{2}\right)^{2}&-\frac{\sin(\xi)}{2}&0\\\
0&-\frac{\sin(\xi)}{2}&\sin\left(\frac{\xi}{2}\right)^{2}&0\\\ 0&0&0&0\\\
\end{array}\right),$ (72)

with

$\xi=\cos^{-1}\left(\frac{E_{Z}-E_{V}+g_{xy,F}+g_{z,H}}{g_{z,F}-g_{z,H}}\right).$
(73)

Clearly, the phase exists only for the cosine argument having a magnitude
smaller than unity. The angles describing this phase are

$\displaystyle\alpha_{1}=$ $\displaystyle 0;\ \alpha_{2}=\xi$ (74a)
$\displaystyle\theta_{p}=$ $\displaystyle\theta_{s}=0$ (74b)
$\displaystyle\beta_{1}=$ $\displaystyle\beta_{2}=\pi$ (74c)

The HF energy is $E_{SVE+}=$

$\frac{E_{V}^{2}-2E_{V}(E_{Z}+g_{xy,F}+g_{z,F})+E_{Z}^{2}+2E_{Z}(g_{xy,F}-g_{z,F})+4g_{z,H}(E_{Z}+g_{xy,F})+g_{xy,F}^{2}+2g_{z,F}(g_{z,H}-g_{xy,F})-g_{z,F}^{2}}{2(g_{z,F}-g_{z,H})}$
(75)

The Hessian matrix is

${\cal E}_{SVE+}=\left(\begin{array}[]{cc}A_{2\times 2}&0_{2\times 2}\\\
0_{2\times 2}&B_{2\times 2}\end{array}\right)_{4\times 4}$ (76)

where,

$A=\left(\begin{array}[]{cc}\frac{E_{V}(g_{xy,F}+g_{z,F})-(E_{Z}+g_{xy,F}+g_{z,F})(g_{xy,F}-g_{z,F}+2g_{z,H})}{g_{z,F}-g_{z,H}}&0\\\
0&\frac{(E_{V}-E_{Z}-g_{xy,F}-g_{z,F})(E_{V}-E_{Z}-g_{xy,F}+g_{z,F}-2g_{z,H})}{g_{z,F}-g_{z,H}}\\\
\end{array}\right)$ (77)

and

$B=\left(\begin{array}[]{cc}C&0\\\ 0&D\end{array}\right)$ (78)

with
$C=-\frac{(E_{V}-E_{Z}-g_{xy,F}-g_{z,F})\left(E_{V}(g_{xy,F}-g_{xy,H})+E_{Z}(-g_{xy,F}+g_{xy,H}+g_{z,F}-g_{z,H})-g_{xy,F}^{2}+g_{xy,F}g_{xy,H}+g_{xy,F}g_{z,F}-2g_{xy,F}g_{z,H}+g_{xy,H}g_{z,F}\right)}{(g_{z,F}-g_{z,H})^{2}}$
and $D=\frac{E_{Z}(E_{V}-E_{Z}-g_{xy,F}+g_{z,F}-2g_{z,H})}{g_{z,F}-g_{z,H}}$.

The instabilities of SVE+ phase are given by

$\displaystyle g_{z,F}=$ $\displaystyle E_{V}-E_{Z}-g_{xy,F}$ (79a)
$\displaystyle g_{z,H}=$
$\displaystyle\frac{1}{2}(E_{V}-E_{Z}-g_{xy,F}+g_{z,F})$ (79b) $\displaystyle
g_{z,H}=$
$\displaystyle\frac{(g_{xy,F}-g_{xy,H})(E_{V}-E_{Z}-g_{xy,F})+g_{z,F}(E_{Z}+g_{xy,F}+g_{xy,H})}{E_{Z}+2g_{xy,F}}$
(79c) $\displaystyle E_{V}=$
$\displaystyle\frac{(E_{Z}+g_{xy,F}+g_{z,F})(g_{xy,F}-g_{z,F}+2g_{z,H})}{g_{xy,F}+g_{z,F}}.$
(79d)

The first corresponds to the instability towards the FM phase, while the
second corresponds to the instability towards the CDW phase. The third is the
instability towards the B/CAF phase, while the fourth seems to not be
relevant.

##  Appendix G Hartree and Fock couplings having opposite signs

When Landau-level mixing is strong, it could happen that the Hartree and Fock
parts of a given coupling are of opposite sign. We will present sample results
for the three cases (i) $g_{z,F}/g_{z,H}<0;\ g_{xy,F}/g_{xy,H}>0$. (ii)
$g_{z,F}/g_{z,H}>0;\ g_{xy,F}/g_{xy,H}<0$. (iii) $g_{z,F}/g_{z,H}<0;\
g_{xy,F}/g_{xy,H}<0$.

We have found that if the ratio F/H of a particular coupling is negative, the
magnitudes of H and F parts of that particular coupling seem to be irrelevant
to the phase diagram. This is in contrast to the case when the ratio F/H is
positive, where it matters a great deal whether the ratio is bigger or smaller
than unity.

Consider first case (i): $g_{z,F}/g_{z,H}<0;\ g_{xy,F}/g_{xy,H}>0$. We set
$g_{z,F}=-g_{z,H}$, and vary the ratio $g_{xy,F}/g_{xy,H}$.

Figure 25: Phase diagram for $g_{z,F}/g_{z,H}=-1$, and
$g_{xy,F}/g_{xy,H}=0.75$, $E_{Z}=1$ and $E_{V}=0.0$. The topology of the phase
diagram changes; in particular, the CAF region no longer extends to infinity,
but is enclosed by the other phases. There is no coexistence anywhere in the
phase diagram.

The phase diagram for $g_{xy,F}/g_{xy,H}<1$ is shown in Fig. 25. The topology
of the phase diagram changes dramatically, but there is no coexistence
anywhere. The phase diagram for $g_{xy,F}/g_{xy,H}>1$ is shown in Fig. 26), in
which there are two coexistence regions. The first one is the B/CAF phase,
with BO, FM, CAF, and SVEY order, which we already encountered when
$g_{z,F}/g_{z,H}>0$. The second phase, which we call the FSVE phase, only
occurs when $g_{z,F}/g_{z,H}<0$.

Figure 26: Phase diagram for $g_{z,F}/g_{z,H}=-1$, and
$g_{xy,F}/g_{xy,H}=1.25$, $E_{Z}=1$ and $E_{V}=0.0$. The blue shaded region is
the B/CAF phase, which has BO, FM, CAF, and SVEY order. The brown shaded
region is a new coexistence phase with FM and either of SVEX or SVEY order.

At $E_{V}=0$ this phase shows the coexistence of FM and either SVEX or SVEY
order (they are degenerate). Fig. 27 shows the order parameters along a
horizontal section of the phase diagram at $g_{z,H}=0.8$.

Figure 27: Nonzero order parameters vs $g_{xy,H}$ for $g_{z,H}=0.8$. All other
coupling constants are the same as in Fig. 26. The B/CAF phase shows BO, FM,
and SVEY (=SVE-) order, while the FSVE phase shows FM and SVEY order. Figure
28: Order parameters vs $g_{xy,H}$ for various values of $E_{V}$, with
$E_{Z}=1.0,~{}g_{z,F}=-1.0,~{}g_{z,H}=0.8,~{}g_{xy,F}=1.25g_{xy,H}$. (a) For a
tiny $E_{V}=0.05$ the FSVE phase splits into two phases. One of them has CDW,
FM, and unequal SVEX and SVEY order parameters, while the other is the
familiar SVE+ phase with CDW, FM, and SVE+ order.(b) As $E_{V}$ increases to
$0.2$, the B/CAF and FSVE phases shrink, while the SVE+ phase expands at their
expense. (c) At $E_{V}=0.5$ the FSVE phase has vanished, leaving behind the
B/CAF and SVE+ phases.

Fig. 28 shows the changes that occur on this section when $E_{V}$ is turned
on. The B/CAF phase changes in a familiar manner, acquiring a CDW order
parameter as well as unequal SVEX and SVEY order parameters in addition to the
BO and CAF order already present. The FSVE phase splits into two phases, both
having some CDW order. The first has FM and unequal SVEX and SVEY order
parameters, while the second is the familiar SVE+ phase.

Let us now go on to case (ii), $g_{z,F}/g_{z,H}>0;\ g_{xy,F}/g_{xy,H}<0$. We
choose $g_{xy,F}=-1.0*g_{xy,H}$. Fig. 29 shows the phase diagram for
$g_{z,F}/g_{z,H}<1$. The topology again changes dramatically, but there are no
coexistence phases. Fig. 30 shows the phase diagram for $g_{z,F}/g_{z,H}>1$.
The familiar SVE+ phase interpolates between the FM and CDW phases.

Fig. 31 shows the order parameters vs $g_{xy,H}$ along a section of Fig. 30 at
$g_{z,H}=-1.3$. One sees the first-order transition from the BO to the SVE+
phase, and the subsequent second-order transition into the CDW phase.

Figure 29: For this phase diagram we have considered
$E_{Z}=1.0,E_{V}=0.0,g_{z,F}=0.75*g_{z,H},g_{xy,F}=-1.0*g_{xy,H}$. One can
clearly see that as here $g_{xy,F}=-1.0*g_{xy,H}$ CAF is the stable phase for
$g_{xy,H}>\frac{E_{Z}}{2}$. Figure 30: Phase diagram for
$E_{Z}=1.0,E_{V}=0.0,g_{z,F}=1.25*g_{z,H},g_{xy,F}=-1.0*g_{xy,H}$. The
topology of the phase diagram is different from the earlier cases, but now the
SVE+ coexistence phase appears between the FM and CDW phases.  Figure 31:
Nonzero order parameters vs $g_{xy,H}$ for $g_{z,H}=-1.3$. All other coupling
constants are chosen as in Fig. 30. The system is in the BO phase at the
extreme left, makes a first-order transition into the SVE+ phase, and finally,
a second-order transition into the CDW phase at the extreme right.

Finally we turn to case (iii), $g_{z,F}/g_{z,H}<0;\ g_{xy,F}/g_{xy,H}<0$. The
various panels of Fig. 32 show the phase diagrams for different F/H ratios.
The common feature is the presence of the FSVE coexistence phase.

Figure 32: Phase diagrams for case (iii) $g_{z,F}/g_{z,H}<0;\
g_{xy,F}/g_{xy,H}<0$. (a)
$E_{Z}=1.0,E_{V}=0.0,g_{z,F}=-0.75*g_{z,H},g_{xy,F}=-0.75*g_{xy,H}$ (b)
$E_{Z}=1.0,~{}E_{V}=0.0,~{}g_{z,F}=-0.75*g_{z,H},~{}g_{xy,F}=-1.25*g_{xy,H}$
(c)
$E_{Z}=1.0,~{}E_{V}=0.0,~{}g_{z,F}=-1.25*g_{z,H},~{}g_{xy,F}=-1.25*g_{xy,H}$
(d)
$E_{Z}=1.0,~{}E_{V}=0.0,~{}g_{z,F}=-1.25*g_{z,H},~{}g_{xy,F}=-0.75*g_{xy,H}$.
The brown shaded region denotes the FSVE phase in all the figures.

## References

  * Prange and Girvin [1990] R. Prange and S. M. Girvin, _The Quantum Hall Effect_ (Springer Verlag, 1990).
  * Das Sarma and Pinczuk [1996] S. Das Sarma and A. Pinczuk, eds., _Perspectives in Quantum Hall Effects_ (Wiley, 1996).
  * Klitzing et al. [1980] K. v. Klitzing, G. Dorda, and M. Pepper, Phys. Rev. Lett. 45, 494 (1980), URL https://link.aps.org/doi/10.1103/PhysRevLett.45.494.
  * Kane and Mele [2005] C. L. Kane and E. J. Mele, Phys. Rev. Lett. 95, 146802 (2005), URL https://link.aps.org/doi/10.1103/PhysRevLett.95.146802.
  * Hasan and Kane [2010] M. Z. Hasan and C. L. Kane, Rev. Mod. Phys. 82, 3045 (2010), URL https://link.aps.org/doi/10.1103/RevModPhys.82.3045.
  * Halperin [1982] B. I. Halperin, Phys. Rev. B 25, 2185 (1982), URL https://link.aps.org/doi/10.1103/PhysRevB.25.2185.
  * Tsui et al. [1982] D. C. Tsui, H. L. Stormer, and A. C. Gossard, Phys. Rev. Lett. 48, 1559 (1982), URL https://link.aps.org/doi/10.1103/PhysRevLett.48.1559.
  * Laughlin [1983] R. B. Laughlin, Phys. Rev. Lett. 50, 1395 (1983), URL https://link.aps.org/doi/10.1103/PhysRevLett.50.1395.
  * Sondhi et al. [1993] S. L. Sondhi, A. Karlhede, S. A. Kivelson, and E. H. Rezayi, Phys. Rev. B 47, 16419 (1993), URL https://link.aps.org/doi/10.1103/PhysRevB.47.16419.
  * Fertig [1989] H. A. Fertig, Phys. Rev. B 40, 1087 (1989), URL https://link.aps.org/doi/10.1103/PhysRevB.40.1087.
  * Yang et al. [1994] K. Yang, K. Moon, L. Zheng, A. H. MacDonald, S. M. Girvin, D. Yoshioka, and S.-C. Zhang, Phys. Rev. Lett. 72, 732 (1994), URL https://link.aps.org/doi/10.1103/PhysRevLett.72.732.
  * Moon et al. [1995] K. Moon, H. Mori, K. Yang, S. M. Girvin, A. H. MacDonald, L. Zheng, D. Yoshioka, and S.-C. Zhang, Phys. Rev. B 51, 5138 (1995), URL https://link.aps.org/doi/10.1103/PhysRevB.51.5138.
  * Berger et al. [2004] C. Berger, Z. Song, T. Li, X. Li, A. Y. Ogbazghi, R. Feng, Z. Dai, A. N. Marchenkov, E. H. Conrad, P. N. First, et al., The Journal of Physical Chemistry B 108, 19912 (2004), ISSN 1520-6106, URL https://doi.org/10.1021/jp040650f.
  * Novoselov [2004] K. S. Novoselov, Science 306, 666 (2004), URL https://doi.org/10.1126/science.1102896.
  * Zhang et al. [2005] Y. Zhang, J. P. Small, M. E. S. Amori, and P. Kim, Phys. Rev. Lett. 94, 176803 (2005), URL https://link.aps.org/doi/10.1103/PhysRevLett.94.176803.
  * Castro Neto et al. [2009] A. H. Castro Neto, F. Guinea, N. M. R. Peres, K. S. Novoselov, and A. K. Geim, Rev. Mod. Phys. 81, 109 (2009), URL https://link.aps.org/doi/10.1103/RevModPhys.81.109.
  * Brey and Fertig [2006a] L. Brey and H. A. Fertig, Phys. Rev. B 73, 195408 (2006a), URL https://link.aps.org/doi/10.1103/PhysRevB.73.195408.
  * Alicea and Fisher [2006] J. Alicea and M. P. A. Fisher, Phys. Rev. B 74, 075422 (2006), URL https://link.aps.org/doi/10.1103/PhysRevB.74.075422.
  * Yang et al. [2006] K. Yang, S. Das Sarma, and A. H. MacDonald, Phys. Rev. B 74, 075423 (2006), URL https://link.aps.org/doi/10.1103/PhysRevB.74.075423.
  * Herbut [2007a] I. F. Herbut, Phys. Rev. B 75, 165411 (2007a), URL https://link.aps.org/doi/10.1103/PhysRevB.75.165411.
  * Herbut [2007b] I. F. Herbut, Phys. Rev. B 76, 085432 (2007b), URL https://link.aps.org/doi/10.1103/PhysRevB.76.085432.
  * Abanin et al. [2006] D. A. Abanin, P. A. Lee, and L. S. Levitov, Phys. Rev. Lett. 96, 176803 (2006), URL https://link.aps.org/doi/10.1103/PhysRevLett.96.176803.
  * Murakami et al. [2004] S. Murakami, N. Nagaosa, and S.-C. Zhang, Phys. Rev. Lett. 93, 156804 (2004), URL https://link.aps.org/doi/10.1103/PhysRevLett.93.156804.
  * Brey and Fertig [2006b] L. Brey and H. A. Fertig, Physical Review B 73 (2006b), ISSN 1550-235X, URL http://dx.doi.org/10.1103/PhysRevB.73.195408.
  * Zhang et al. [2006] Y. Zhang, Z. Jiang, J. P. Small, M. S. Purewal, Y.-W. Tan, M. Fazlollahi, J. D. Chudow, J. A. Jaszczak, H. L. Stormer, and P. Kim, Phys. Rev. Lett. 96, 136806 (2006), URL https://link.aps.org/doi/10.1103/PhysRevLett.96.136806.
  * Jiang et al. [2007] Z. Jiang, Y. Zhang, H. L. Stormer, and P. Kim, Phys. Rev. Lett. 99, 106802 (2007), URL https://link.aps.org/doi/10.1103/PhysRevLett.99.106802.
  * Young et al. [2012] A. F. Young, C. R. Dean, L. Wang, H. Ren, P. Cadden-Zimansky, K. Watanabe, T. Taniguchi, J. Hone, K. L. Shepard, and P. Kim, Nature Physics 8, 550–556 (2012), ISSN 1745-2481, URL http://dx.doi.org/10.1038/nphys2307.
  * Maher et al. [2013] P. Maher, C. R. Dean, A. F. Young, T. Taniguchi, K. Watanabe, K. L. Shepard, J. Hone, and P. Kim, Nature Physics 9, 154–158 (2013), ISSN 1745-2481, URL http://dx.doi.org/10.1038/nphys2528.
  * Young et al. [2014] A. F. Young, J. D. Sanchez-Yamagishi, B. Hunt, S. H. Choi, K. Watanabe, T. Taniguchi, R. C. Ashoori, and P. Jarillo-Herrero, Nature 505, 528 (2014), ISSN 1476-4687, URL https://doi.org/10.1038/nature12800.
  * Stefanidis and Villadiego [2022] N. Stefanidis and I. S. Villadiego, _Theory of broken symmetry quantum hall states in the $n=1$ landau level of graphene_ (2022), URL https://arxiv.org/abs/2210.03752.
  * Kharitonov [2012a] M. Kharitonov, Phys. Rev. B 85, 155439 (2012a), URL https://link.aps.org/doi/10.1103/PhysRevB.85.155439.
  * Zibrov et al. [2018] A. Zibrov, E. Spanton, H. Zhou, C. Kometter, T. Taniguchi, K. Watanabe, and A. Young, Nature Physics 14 (2018).
  * 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–233 (2018), ISSN 1095-9203, URL http://dx.doi.org/10.1126/science.aar4061.
  * Zhou et al. [2019] H. Zhou, H. Polshyn, T. Taniguchi, K. Watanabe, and A. F. Young, Nature Physics 16, 154–158 (2019), ISSN 1745-2481, URL http://dx.doi.org/10.1038/s41567-019-0729-8.
  * Assouline et al. [2021] A. Assouline, M. Jo, P. Brasseur, K. Watanabe, T. Taniguchi, T. Jolicoeur, D. C. Glattli, N. Kumada, P. Roche, F. D. Parmentier, et al., Nature Physics 17, 1369–1374 (2021), ISSN 1745-2481, URL http://dx.doi.org/10.1038/s41567-021-01411-z.
  * 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, Physical Review X 12 (2022), URL https://doi.org/10.1103%2Fphysrevx.12.021060.
  * Li et al. [2019] S.-Y. Li, Y. Zhang, L.-J. Yin, and L. He, Phys. Rev. B 100, 085437 (2019), URL https://link.aps.org/doi/10.1103/PhysRevB.100.085437.
  * 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–326 (2022), ISSN 1095-9203, URL http://dx.doi.org/10.1126/science.abm3770.
  * 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), URL https://doi.org/10.1038%2Fs41586-022-04513-7.
  * Hunt et al. [2013] B. Hunt, J. D. Sanchez-Yamagishi, A. F. Young, M. Yankowitz, B. J. LeRoy, K. Watanabe, T. Taniguchi, P. Moon, M. Koshino, P. Jarillo-Herrero, et al., Science 340, 1427 (2013), URL https://www.science.org/doi/abs/10.1126/science.1237240.
  * Amet et al. [2013] F. Amet, J. R. Williams, K. Watanabe, T. Taniguchi, and D. Goldhaber-Gordon, Phys. Rev. Lett. 110, 216601 (2013), URL https://link.aps.org/doi/10.1103/PhysRevLett.110.216601.
  * Jung et al. [2015] J. Jung, A. M. DaSilva, A. H. MacDonald, and S. Adam, Nature Communications 6, 6308 (2015), ISSN 2041-1723, URL https://doi.org/10.1038/ncomms7308.
  * Jung et al. [2017] J. Jung, E. Laksono, A. M. DaSilva, A. H. MacDonald, M. Mucha-Kruczyński, and S. Adam, Phys. Rev. B 96, 085442 (2017), URL https://link.aps.org/doi/10.1103/PhysRevB.96.085442.
  * Murthy and Shankar [2002] G. Murthy and R. Shankar, Phys. Rev. B 65, 245309 (2002), URL https://link.aps.org/doi/10.1103/PhysRevB.65.245309.
  * Bishara and Nayak [2009] W. Bishara and C. Nayak, Phys. Rev. B 80, 121302 (2009), URL https://link.aps.org/doi/10.1103/PhysRevB.80.121302.
  * Sodemann and MacDonald [2013] I. Sodemann and A. H. MacDonald, Phys. Rev. B 87, 245425 (2013), URL https://link.aps.org/doi/10.1103/PhysRevB.87.245425.
  * Peterson and Nayak [2013] M. R. Peterson and C. Nayak, Phys. Rev. B 87, 245129 (2013), URL https://link.aps.org/doi/10.1103/PhysRevB.87.245129.
  * Peterson and Nayak [2014] M. R. Peterson and C. Nayak, Phys. Rev. Lett. 113, 086401 (2014), URL https://link.aps.org/doi/10.1103/PhysRevLett.113.086401.
  * Das et al. [2022] A. Das, R. K. Kaul, and G. Murthy, Phys. Rev. Lett. 128, 106803 (2022), URL https://link.aps.org/doi/10.1103/PhysRevLett.128.106803.
  * Lian et al. [2016] Y. Lian, A. Rosch, and M. O. Goerbig, Phys. Rev. Lett. 117, 056806 (2016), URL https://link.aps.org/doi/10.1103/PhysRevLett.117.056806.
  * Lian and Goerbig [2017] Y. Lian and M. O. Goerbig, Phys. Rev. B 95, 245428 (2017), URL https://link.aps.org/doi/10.1103/PhysRevB.95.245428.
  * Atteia and Goerbig [2021] J. Atteia and M. O. Goerbig, Phys. Rev. B 103, 195413 (2021), URL https://link.aps.org/doi/10.1103/PhysRevB.103.195413.
  * Murthy et al. [2017] G. Murthy, E. Shimshoni, and H. A. Fertig, Phys. Rev. B 96, 245125 (2017), URL https://link.aps.org/doi/10.1103/PhysRevB.96.245125.
  * McCann and Fal’ko [2006] E. McCann and V. I. Fal’ko, Phys. Rev. Lett. 96, 086805 (2006), URL https://link.aps.org/doi/10.1103/PhysRevLett.96.086805.
  * Douçot et al. [2008] B. Douçot, M. O. Goerbig, P. Lederer, and R. Moessner, Phys. Rev. B 78, 195327 (2008), URL https://link.aps.org/doi/10.1103/PhysRevB.78.195327.
  * Atteia et al. [2021] J. Atteia, Y. Lian, and M. O. Goerbig, Phys. Rev. B 103, 035403 (2021), URL https://link.aps.org/doi/10.1103/PhysRevB.103.035403.
  * Min et al. [2006] H. Min, J. E. Hill, N. A. Sinitsyn, B. R. Sahu, L. Kleinman, and A. H. MacDonald, Phys. Rev. B 74, 165310 (2006), URL https://link.aps.org/doi/10.1103/PhysRevB.74.165310.
  * Yao et al. [2007] Y. Yao, F. Ye, X.-L. Qi, S.-C. Zhang, and Z. Fang, Phys. Rev. B 75, 041401 (2007), URL https://link.aps.org/doi/10.1103/PhysRevB.75.041401.
  * Huertas-Hernando et al. [2006] D. Huertas-Hernando, F. Guinea, and A. Brataas, Phys. Rev. B 74, 155426 (2006), URL https://link.aps.org/doi/10.1103/PhysRevB.74.155426.
  * Zarea and Sandler [2007] M. Zarea and N. Sandler, Phys. Rev. Lett. 99, 256804 (2007), URL https://link.aps.org/doi/10.1103/PhysRevLett.99.256804.
  * Kharitonov [2012b] M. Kharitonov, Phys. Rev. Lett. 109, 046803 (2012b), URL https://link.aps.org/doi/10.1103/PhysRevLett.109.046803.
  * Kharitonov [2012c] M. Kharitonov, Phys. Rev. B 86, 195435 (2012c), URL https://link.aps.org/doi/10.1103/PhysRevB.86.195435.
  * Kosterlitz and Thouless [1973] J. M. Kosterlitz and D. J. Thouless, Journal of Physics C: Solid State Physics 6, 1181 (1973), URL https://doi.org/10.1088/0022-3719/6/7/010.
  * Berezinskii [1972] V. L. Berezinskii, Sov. Phys. JETP 34, 610 (1972), URL http://www.jetp.ras.ru/cgi-bin/e/index/r/61/3/p1144?a=list.
  * Fu et al. [2021] H. Fu, K. Huang, K. Watanabe, T. Taniguchi, and J. Zhu, Phys. Rev. X 11, 021012 (2021), URL https://link.aps.org/doi/10.1103/PhysRevX.11.021012.
  * Champagne et al. [2008] A. R. Champagne, J. P. Eisenstein, L. N. Pfeiffer, and K. W. West, Phys. Rev. Lett. 100, 096801 (2008), URL https://link.aps.org/doi/10.1103/PhysRevLett.100.096801.
  * Aleiner et al. [2007] I. L. Aleiner, D. E. Kharzeev, and A. M. Tsvelik, Phys. Rev. B 76, 195415 (2007), URL https://link.aps.org/doi/10.1103/PhysRevB.76.195415.
  * Basko and Aleiner [2008] D. M. Basko and I. L. Aleiner, Phys. Rev. B 77, 041409 (2008), URL https://link.aps.org/doi/10.1103/PhysRevB.77.041409.
  * Shankar [1994] R. Shankar, Rev. Mod. Phys. 66, 129 (1994), URL https://link.aps.org/doi/10.1103/RevModPhys.66.129.
  * Sodemann and MacDonald [2014] I. Sodemann and A. H. MacDonald, Phys. Rev. Lett. 112, 126804 (2014), URL https://link.aps.org/doi/10.1103/PhysRevLett.112.126804.
  * Hegde and Villadiego [2022] S. S. Hegde and I. S. Villadiego, Physical Review B 105 (2022), URL https://doi.org/10.1103%2Fphysrevb.105.195417.