of 27
Enhanced superconductivity in spin–orbit
proximitized bilayer graphene
In the format provided by the
authors and unedited
Nature | www.nature.com/nature
Supplementary information
https://doi.org/10.1038/s41586-022-05446-x
Supplementary Information:
Enhanced Superconductivity in Spin-Orbit Proximitized Bilayer Graphene
Yiran Zhang, Robert Polski, Alex Thomson,
́
Etienne Lantagne-Hurtubise, Cyprian Lewandowski,
Haoxin Zhou, Kenji Watanabe, Takashi Taniguchi, Jason Alicea, and Stevan Nadj-Perge
1
Supplementary Figures
SI Fig. 1
|
n
-
D
phase diagram of device D4 including electron doping.
R
xx
versus doping
density
n
and displacement field
D
including electron doping for device D4. A node-like resistive
feature appears around
n
= 0
and
D
= 0
(marked by a black arrow; zoom-in data on the right),
showing a small offset in
D
field (
D/ε
0
0
.
025
±
0
.
005
V/nm).
2
SI Fig. 2
|
Alignment between BLG and WSe
2
for D1. a
c
, optical image of WSe
2
(
a
), BLG
(
b
) and alignment of BLG-WSe
2
(
c
) for device D1. Coloured dashed lines indicate the possible
crystal edges with zigzag or armchair orientation. Black dashed lines in
b
and
c
indicate one edge
of BLG. The scale bar in each panel corresponds to
20
μ
m.
Theoretical Analysis
1 Continuum model band structure of bilayer graphene
We consider the low-energy continuum model commonly used to describe Bernal-stacked bilayer
graphene (BLG)
S1
, under a perpendicular displacement field
D
which generates a potential differ-
ence
u
=
d
D/ε
BLG
between the top and bottom layers. Here
d
= 0
.
33
nm is the interlayer
distance and
ε
BLG
4
.
3
is the relative permittivity of bilayer graphene. A continuum approxima-
tion of the band structure returns a Hamiltonian of the form
H
0
=
X
ξ
=
±
X
k
ψ
ξ
(
k
)
h
0
(
k
)
ψ
ξ
(
k
)
, h
0
(
k
) =
u/
2
v
0
Π
v
4
Π
v
3
Π
v
0
Π
+
u/
2
γ
1
v
4
Π
v
4
Π
γ
1
u/
2
v
0
Π
v
3
Π
v
4
Π
v
0
Π
u/
2
(1)
where
Π = (
ξk
x
+
ik
y
)
and
v
i
3
a
2
γ
i
. Here,
ξ
=
±
1
indicates the valley that has been expanded
about:
K
,
K
= (
ξ
4
π/
3
a,
0)
with
a
= 0
.
246
nm the lattice constant of monolayer graphene. The
4
×
4
matrix
h
ξ
(
k
)
is expressed in the sublattice/layer basis corresponding to creation/annihilation
operators of the form
ψ
ξ
(
k
) = (
ψ
ξ,A
1
(
k
)
ξ,B
1
(
k
)
ξ,A
2
(
k
)
ξ,B
2
(
k
))
T
, where
A
/
B
indicate the
sublattice,
1
,
2
indicate the layer, and the momentum
k
is measured relative to
K
ξ
(indices denot-
ing the spin degrees of freedom have been suppressed). It will sometimes be convenient below to
express the Hamiltonian in terms of the spinors
ψ
(
k
) = (
ψ
+
(
k
)
(
k
))
T
.
The common values quoted for the five parameters entering the continuum model in Eq. (1) are
3
γ
0
= 2
.
61
eV (intralayer nearest-neighbor tunneling),
γ
1
= 361
meV (leading interlayer tunnel-
ing),
γ
3
= 283
meV (also known as trigonal warping term),
γ
4
= 138
meV, and
= 15
meV
(potential difference between dimer and non-dimer sites)
S2
.
A TMD monolayer adjacent to the graphene, such as is the case here with WSe
2
, is known to
induce SOC via virtual tunnelling
S3–S5
:
H
SOC
=
X
ξ
=
±
X
k
ψ
ξ
(
k
)
h
SOC
ψ
ξ
(
k
)
, h
SOC
(
k
) =
P
1

λ
I
2
ξs
z
+
λ
R
2
ξσ
x
s
y
σ
y
s
x
)

,
(2)
where the Pauli matrices
σ
i
and
s
i
,
i
=
x,y,z
, respectively act on sublattice and spin degrees of
freedom. The operator
P
1
projects onto the top graphene sheet, i.e., only the sites A1 and B1:
P
1
= diag(
1
2
×
2
,
0
2
×
2
)
in the layer/sublattice basis used to express
h
0
(
k
)
in (1). The parameters
λ
I
and
λ
R
quantify the strength of the Ising (also called “valley-Zeeman”) and Rashba SOC. Ab
initio-type numerics and experimental estimates find a range of values
λ
I
0
5
meV and
λ
R
0
15
meV for the SOC parameters
S3, S4, S6–S11
, which are also predicted to be strongly
twist-angle dependent
S8–S10
.
In the absence of SOC and an applied displacement field with
v
3
=
v
4
= 0
, two bands touch
quadratically at charge neutrality. Two remaining bands are at significantly higher and lower en-
ergies; their wavefunction are dominated by the “dimer sites,” i.e., the A2 and B1 which sit im-
mediately on top of one another in the bilayer and hybridize strongly through the onsite tunnelling
parameter
γ
1
. Trigonal warping introduced by the
v
3
,
v
4
associated hoppings in Eq. (1) splits the
quadratic band touching at charge neutrality into four distinct Dirac cones separated by van Hove
singularities (VHS): one Dirac cone remains at
k
= 0
, while the other three are located at
C
3
-
related momenta slightly away from the Dirac point. Turning on a displacement field
D
, a gap
opens at charge neutrality and the VHSs move apart in energy. Further, by flattening the band bot-
tom, the applied
D
field also amplifies divergence of the DOS close to the VHS. The low-energy
states near
K
and
K
become strongly layer- and sublattice-polarized; e.g., on
A
1
sites for the va-
lence band and
B
2
sites for the conduction band, or vice versa for the other sign of
D
. That is, the
low-energy wavefunctions near charge neutrality and under a large
D
field are strongly localized
on the “non-dimer sites” of BLG.
The layer- and sublattice polarization of the low-energy wavefunctions near the
K
,
K
points
has important consequences for SOC induced by the TMD. Indeed, Rashba SOC does not act
effectively in the low-energy theory because it is off-diagonal in the sublattice degree of freedom.
It therefore induces a splitting only at second order in degenerate perturbation theory, with
λ
eff
R
(
λ
2
R
v
0
k
)
2
/
(
γ
2
1
u
)
with
u
the interlayer potential
S12
(neglecting effects of further perturbations such
as trigonal warping, which we discuss below). By contrast, the Ising SOC acts effectively in the
subspace of sublattice- and layer-polarized wavefunctions.
4
The normalized frequencies one expects from quantum oscillations for the non-interacting theory
are shown in Extended Data Fig. 6a as a function of hole doping for
D/ε
0
= 1
V/nm. The red
lines correspond to the spin-orbit-free case, whereas the blue lines are computed in the presence of
SOC; we additionally plot the Fermi surfaces for the spin-orbit coupled band structure in the insets
at a few representative fillings. The Ising coupling is set to the experimentally extracted value
λ
I
= 0
.
7
meV, whereas for the Rashba coupling we select
λ
R
= 3
meV, which is within the upper
bound consistent with the experimental resolution. For the filling range shown, five different sets of
Fermi surface topologies are present for the case with SOC. For
|
n
|
<
0
.
5
×
10
11
cm
2
, an
FP(6)
+
state is realized in which each valley contributed three equally sized (but offset in momentum
space) Fermi surfaces. Assuming
λ
I
>
0
, given that the in-plane mixing by Rashba is relatively
small, the states that make up these pockets largely have valley and spin quantum numbers
K
,
and
K,
. With further doping, another six equally sized Fermi surfaces appear corresponding to
the spin degrees of freedom pushed down in energy by SOC (
K,
and
K
,
). This
FP(6
,
6)
+
state persists with doping until
|
n
| ∼
7
.
4
×
10
11
cm
2
, at which point the system reaches a van
Hove singularity and thus the Fermi surface structure changes. The three majority pockets merge,
resulting in two degenerate hole pockets with a small electron pocket at their centre. With further
doping the electron-like pocket of the majority Fermi surfaces vanishes (
|
n
| ∼
9
×
10
11
cm
2
).
Subsequently, the minority pockets reach the same VHS (
|
n
| ∼
9
.
8
×
10
11
cm
2
) leading to the
formation of the small electron pocket.
A comparison of Extended Data Fig. 6a with the quantum oscillations data in Fig. 3c makes it
clear that the non-interacting theory is insufficient. In particular, SOC does explicitly “polarize”
the band—they are energetically split—so that the
FP(6
,
6)
+
phase is realized; this state persists
up to
|
n
| ∼
7
.
4
×
10
11
cm
2
at
D
= 1
V/nm. By contrast, in the experiment the
FP(6
,
6)
+
phase
terminates at around
|
n
| ∼
5
×
10
11
cm
2
, where it is replaced by the
FP(2
,
2)
+
state. Hence,
while the Ising and Rashba parameters (obtained from measurements at zero
D
and doping) do
arguably polarize the bands in the non-interacting limit, the resulting splitting is not large enough
to account for the observed phase diagram.
2 Interactions
The resistance data as a function of displacement field and doping clearly demonstrate that the
non-interacting band structure implied by
H
0
+
H
SOC
in the previous section cannot fully describe
the system. Instead, given the large density of states close to charge neutrality in the presence of
large displacement fields, a series of polarized phases are observed, which are naturally explained
as a consequence of the Coulomb interaction.
The Coulomb interaction is given by
H
C
=
1
2
A
X
k
,
k
,
q
V
(
q
)
ψ
α
(
k
)
ψ
β
(
k
)
ψ
β
(
k
q
)
ψ
α
(
k
+
q
)
.
(3)
5
Here the indices
α
and
β
sum over valley, layer, spin, and sublattice degrees of freedom and
A
=
A
u
.
c
.
N
site
is the total area of the sample, with
A
u
.
c
.
denoting the unit cell area,
A
u
.
c
.
=
3
a
2
/
2
, and
N
site
denoting the total number of sites. The unscreened Coulomb potential is
V
(
q
) =
e
2
/
(2
ε
r
ε
0
|
q
|
)
, where
ε
r
4
.
3
is the dielectric constant for hBN-screened graphene. Instead of
considering this model, we look at a far simpler model in which the interaction is fully local:
V
(
q
) =
A
u
.
c
.
U
C
. We can roughly estimate
U
C
1
A
u
.
c
.
e
2
4
πε
r
ε
0
d,
(4)
where we have substituted
2
π/q
F
d
with
d
the inter-particle spacing. A density of
n
=
5
×
10
11
cm
2
roughly translates to an inter-particle spacing of
d
= 15
nm, which in turn implies
U
C
100
eV. This estimate should be taken as an upper bound since it does not include the effects
of screening. Accordingly, more reasonable results are obtained by allowed the effective Coulomb
interaction strength to take smaller values. In particular, we often select
U
C
= 35
eV in accord
with earlier calculations of Bernal stacked systems
S13, S14
.
We emphasize that
U
C
should
not
be thought of as the setting the “energy scale” of the problem.
Instead, the interactions naturally scale with the density. In particular, if we let
ν
f
denote the
number of electrons per unit cell,
ν
f
=
A
u
.
c
.
·
n
=
3
a
2
/
(2
d
2
)
, then the energy per electron is
ε
C
=
U
C
electron
e
2
4
πε
r
ε
0
1
d
,
(5)
which is precisely what we would have found with a real space description. In this case, we find
ε
C
20
meV for
d
= 15
nm.
Even with the long-range Coulomb form,
V
(
q
)
1
/
|
q
|
, the interaction presented in Eq. (3) is
not fully general. Instead, it was derived by taking the zero momentum portion of the density. In
effect, the density can be expanded in terms of the continuum model operators as
ρ
(
r
)
ρ
++
(
r
) +
ρ
−−
(
r
) +
e
i
K
·
r
ρ
+
(
r
) +
e
i
K
·
r
ρ
+
(
r
)
where
ρ
ξξ
(
r
) =
ψ
ξ
(
r
)
ψ
ξ
(
r
)
with
ψ
ξ
(
r
)
the real space
version of the annihilation operator defined in SI, section 1. Equation (3) only includes
ρ
++
and
ρ
−−
, which accounts for the long-range part of the Coulomb interaction. The Hund’s term includes
the remaining two piece of the density carrying momentum
K
,
K
and thus necessarily has a
minimum momentum transfer of
K
. Its magnitude can therefore be characterized by
V
(
K
) =
e
2
/
(4
πε
r
ε
0
|
K
|
)
. Translating this scale into the relevant energy scale like in (5), we find
ε
H
e
2
4
πε
r
ε
0
1
a
d
a
ε
C
,
(6)
where factors of order unity have been neglected.
6
3 Symmetries
We begin by discussing the flavour symmetries of the Hamiltonian in the absence of SOC. It
immediately follows that the system is invariant under the usual
SU(2)
s
spin rotation symmetry:
ψ
(
k
)
e
n
·
s
/
2
ψ
(
k
)
, where
n
is an arbitrary unit 3-vector. The system similarly preserves the
familiar
U(1)
c
phase rotation symmetry associated with charge conservation,
ψ
(
k
)
e
ψ
(
k
)
.
These two standard symmetries are further augmented in bilayer graphene by the preservation of
particle number individually within each valley, which follows from the so-called
U(1)
v
valley
symmetry; its action takes the form
ψ
(
k
)
e
iθτ
z
ψ
(
k
)
, where
τ
z
is a Pauli matrix acting on the
valley indices of
ψ
(
k
)
. In essence, the
U(1)
v
valley symmetry is a manifestation of the low energy
scales at work: extrinsic scattering between states originating from valley
K
to those originating
from valley
K
are necessarily short range and thus precluded by the high quality of the sample.
Further inspection of the Hamiltonian
H
0
+
H
C
reveals that the physical symmetry group,
U(1)
c
×
U(1)
v
×
SU(2)
s
, is in fact a subgroup of a much larger
effective
symmetry operative at the dominant
energy scales of the system. In particular, the Hamiltonian is invariant under independent
spin
rotations within each valley:
ψ
(
k
)
P
+
e
+
n
+
·
s
/
2
+
P
e
n
·
s
/
2

ψ
(
k
)
, where
n
±
are unit
3-vectors and
P
±
project onto the valley
K
and
K
. Together with the two
U(1)
symmetry groups,
the result is the existence of a
U(2)
K
×
U(2)
K
=
U(1)
c
×
U(1)
v
×
SU(2)
K
×
SU(2)
K
flavour
symmetry. Importantly, this result implies that any degeneracies encoded by the
U(2)
K
×
U(2)
K
effective symmetry are split only at the scale of the Hund’s coupling,
ε
H
(
a/d
)
ε
c
.
The introduction of SOC naturally reduces both the effective and physical symmetry groups. The
Ising term by itself (
λ
R
= 0
) reduces the
SU(2)
s
to
U(1)
z
, the group generating rotations about the
spin-
z
axis, and the remaining physical symmetry group is thus
U(1)
c
×
U(1)
v
×
U(1)
z
. The large
effective symmetry group relevant to scales larger than
ε
H
is similarly diminished by the restriction
that only spin-
z
rotations in either valley leave Hamiltonian unmodified. The result is an effective
symmetry group composed of four different
U(1)
rotations:
U(1)
c
×
U(1)
v
×
U(1)
z,K
×
U(1)
z,K
,
where
U(1)
z,K
(
)
rotates the spin of the valley
K
(
)
fermions about the
z
-axis.
When Rashba spin-orbit coupling is present, with or without Ising SOC, all global, continuous spin
rotations are absent. Both the physical and effective flavour symmetry groups are pared down to
U(1)
c
×
U(1)
v
. The small upper bound for the Rashba coupling
λ
R
imposed by experiment leads
us to largely neglect its symmetry breaking effect.
The Hamiltonian also possesses a number of discrete symmetries, the most important of which is
time reversal symmetry (TRS):
T
:
ψ
(
k
)
x
s
y
ψ
(
k
)
,
i
→−
i.
(7)
Time reversal remains a good symmetry of the system both with and without spin-orbit coupling.
7
4 Mean field approximation
We study the interacting theory using mean field theory. In particular, for each filling
ν
, where
ν
is the number of carriers per unit cell as measured relative to charge neutrality, we find the
Slater determinant ground state
|
φ
ν
that minimize the mean-field ground state energy
E
(
ν
)
MF
[
φ
] =
φ
ν
|
(
H
0
+
H
SOC
+
H
C
)
|
φ
ν
. This procedure is essentially equivalent to replacing the interacting
Hamiltonian with the one-particle mean field Hamiltonian
H
(
ν
)
MF
=
X
k
ψ
(
k
)
h
MF
ψ
(
k
)
,
h
(
ν
)
MF
=
U
C
N
site
X
q
h
P
(
ν
)
(
q
)
tr
P
(
ν
)
(
q
)

1
i
,
(8)
where the projector
P
(
k
)
is given by
P
(
ν
)
αβ
(
q
) =
ψ
β
(
q
)
ψ
α
(
q
)
ν
=
ψ
β
(
q
)
ψ
α
(
q
)
⟩−⟨
ψ
β
(
q
)
ψ
α
(
q
)
CNP
.
(9)
Here, the values of the correlation functions
⟨·⟩
are in turn obtained by diagonalizing
H
0
+
H
SOC
+
H
(
ν
)
MF
, and the subscript “CNP” indicates that the expectation value is being taken with respect to
the charge neutrality point. Self-consistency is attained when the mean field term
H
(
ν
)
MF
used to
calculate
P
(
ν
)
is in turn defined via Eq. (8). It can be shown that the ground state of this self-
consistent Hamiltonian is a local minimum of the mean field energy functional
E
(
ν
)
MF
.
We solve for
P
(
ν
)
through the following procedure. We select an initial value
H
init
.
MF
and then
iterate between Eqs. (8) and (9) until self-consistency is reached. Crucially, states possessing less
symmetry than the initial Hamiltonian
H
0
+
H
SOC
+
H
init
.
MF
are inaccessible. For instance, if the
initial mean field Hamiltonian is invariant under the
U(1)
v
symmetry, the final wavefunction
|
φ
ν
(and the corresponding
P
(
ν
)
) must also be invariant under the
U(1)
v
symmetry and hence so must
H
(
ν
)
MF
. As noted, the symmetries allow us to separate the problem into those that preserve the
U(1)
v
and those that break it.
In principle, the result should be the minimal energy state that respects the same symmetries as
H
init
.
MF
and the non-interacting terms,
H
0
+
H
SOC
. However, in practice, the algorithm sketched
above sometimes finds itself trapped in local minima, unable to attain the true ground state within
that symmetry class. This happenstance is particularly common when there are many nearly de-
generate ground states, which, as we describe in the following section, is the case here. We have
not rigorously explored the phase diagram to ensure that all of the solutions presented below repre-
sent true symmetry-class ground states largely because the simplicity of the model makes it more
appropriate for a qualitative study of trends, as opposed to a quantitative one. There is therefore
little reason to ignore low energy states in favour of what, according to this model, is the “true”
ground state. In fact, the phenomenological arguments we make below in SI, section 7 imply that
a different ground state is realized than suggested by our simulations.
We are primarily interested in what happens upon hole doping the system in the presence of a
positive displacement field, and we therefore specialize to this scenario; our discussion can readily
8
be translated to the case with opposite
D
-field sign as well as with electron doping. We further
note that provided the
D
-induced gap at charge neutrality is sufficiently large, we do not expect
H
C
to induce significant mixing between the four sets of (effectively) degenerate spin-valley bands
defined by
H
0
+
H
SOC
, allowing us to restrict our focus entirely to the active bands of interest.
While the mean field Hamiltonian
h
(
ν
)
MF
is independent of momentum, it nevertheless acts on the
original 16 degrees of freedom as opposed to the four bands of interest. It is convenient to distill
the resulting Hamiltonian to the only degrees of freedom that remain upon projecting to the bands
of interest. In particular, instead of directly discussing
h
(
ν
)
MF
, we focus instead on
h
(
ν
)
MF
=
X
a,i
=0
,x,y,z
(
a,i
)
̸
=(0
,
0)
t
ai
τ
a
s
i
,
t
ai
=
1
4
tr
h
(
ν
)
MF
τ
a
s
i

,
(10)
where we do not include the constant shift of the chemical potential. Notably,
h
(
ν
)
MF
is still a
16
×
16
matrix, but with any dependence on either the layer or sublattice removed. It follows that
h
(
ν
)
MF
does
not account for some of the spatial dependence that results when one projects onto the four hole
bands close to charge neutrality. These effects, while present in our numerics, are largely irrelevant
for the purpose of understanding the resulting phases.
5 Polarized phases
We are most interested here in the spontaneous breaking of the spin-valley symmetries, resulting in
the polarized states seen in experiment. The propensity for this type of symmetry breaking follows
from the large density of states induced by the displacement field. Interactions make having a large
density of states at the Fermi energy energetically costly. At the expense of the kinetic energy, the
DOS at the Fermi energy and its associated energy cost may be lowered by breaking the flavour
symmetry and alternately increasing and decreasing the filling of certain flavours. The advantage
of this process is roughly encapsulated in the Stoner criterion, which states that polarization occurs
when
V ρ
1
, where
V
is the interaction scale and
ρ
the density of states.
At the mean field level, the polarized phases are characterized by the (simplified) mean field Hamil-
tonian
h
(
ν
)
MF
of Eq. (10). We begin by addressing the phases in the absence of SOC where the effec-
tive
U(2)
K
×
U(2)
K
symmetry remains a good approximation. In the simplest case, only a single
t
ai
in Eq. (10) is non-zero:
h
(
ν
)
MF
=
t
ai
τ
a
s
i
.
(11)
This mean field term pushes two flavours up and two flavours down in energy, resulting in a set of
minority and a set of majority Fermi pockets. In what follows, this type of phase is denoted “singly
polarized.” Such singly polarized phases are not limited by mean field Hamiltonians of the form
Eq. (11), but are more generally induced by any
h
(
ν
)
MF
given by a sum of
anticommuting
matrices
τ
a
s
i
.
9
We group the singly polarized phases resulting from Eq. (11) into two broad categories. First are
the “simple” polarized states that preserve the
U(1)
v
valley symmetry, implying that the mean field
Hamiltonian associated with such states satisfies
[
τ
z
,h
(
ν
)
MF
] = 0
:
̃
h
(
ν
)
MF
τ
a
s
i
where
τ
a
s
i
∈{
s
x,y,z
, τ
z
s
0
,x,y,z
}
.
(12)
Notably, the
h
(
ν
)
MF
above commutes with the non-interacting Hamiltonian
h
0
(
k
)
, and it follows that
its primary effect is to generate a relative shift of the band energies. The action of
U(2)
K
×
U(2)
K
rotates the orders leading to spin polarization (SP) (
h
(
ν
)
MF
s
x,y,z
) and spin-valley polarization
(SVP) (
h
(
ν
)
MF
τ
z
s
x,y,z
) into each other, hence these states must have the same energy (with respect
to the Hamiltonian under consideration currently). Similarly, the action of
U(2)
K
×
U(2)
K
may
also rotate
h
(
ν
)
MF
to a linear combination of these order parameters, for instance to induce spin
polarization along an arbitrary direction. By contrast, the valley polarized (VP) state characterized
by
h
(
ν
)
MF
τ
z
does not break
U(2)
K
×
U(2)
K
(although it does break time reversal), and it is
therefore not necessarily degenerate with the SP and SVP states. However, the density-density
form of the Coulomb interaction renders this distinction meaningless and prevents the system
from distinguishing whether valley
K
or
K
is filled on average. Additional interaction terms—
say arising from short-range interactions—will split this accidental degeneracy. For instance, the
phonon interaction
S15
takes the form
R
r
(
ψ
τ
z
ψ
)
2
and thus both preserves the
U(2)
K
×
U(2)
K
symmetry while clearly distinguishing between VP and SP/SVP states. (The Hund’s term whose
inclusion does decrease the effective symmetry group also distinguishes these two sets of states.)
The second category of states breaks the
U(2)
K
×
U(2)
K
to a
diagonal
subgroup through the
spontaneous generation of inter-valley tunnelling. These “inter-valley coherent” (IVC) ordered
states occur when
[
τ
z
,h
(
ν
)
MF
]
̸
= 0
:
̃
h
(
ν
)
MF
τ
a
s
i
where
τ
a
s
i
∈{
τ
x
s
0
,x,y,z
y
s
0
,x,y,z
}
.
(13)
As with the SP and SVP states, the IVC order parameters may all be mapped to one another
through the action of
U(2)
K
×
U(2)
K
, meaning that they must be degenerate. Unlike the previous
set of states, the IVC mean field Hamiltonian mixes states from different valleys and therefore
significantly alters the form of the band structure.
In addition the singly polarized states, the system may also favour breaking more than a single sym-
metry, resulting in a “multiply polarized” state. In this case,
h
(
ν
)
MF
is a sum of multiple
commuting
τ
a
s
i
matrices. An example of such a mean field term is
h
(
ν
)
MF
=
t
0
τ
z
+
t
1
s
z
+
t
3
τ
z
s
z
.
(14)
When
|
t
1
|
=
|
t
2
|
=
|
t
3
|
, this mean field Hamiltonian pushes one flavour to a higher or lower energy
on average, leaving the remaining three degenerate. More commonly, however, we find that the
coefficients satisfy
|
t
1
|
>
|
t
2
|
=
|
t
3
|
. As with the “singly polarized” state above, the multiply
polarized states may also be categorized depending on whether they break or preserve
U(1)
v
.
10
We calculated the self-consistent mean field Hamiltonians and corresponding ground states in the
absence of SOC for a variety of parameters. For filling ranges that prefer singly polarized states, we
consistently find IVC ordered states to have the lowest energies. In this case, no other symmetries
are broken. Similarly, for filling ranges preferring multiply polarized states, IVC order is typically
also generated, although this time it must be present alongside another symmetry breaking order.
We stress, however, that our model is very crude and is not expected to yield quantitatively accurate
results.
We finally turn to the case of primary interest: bilayer graphene with proximity-induced SOC.
Given the relative smallness of the effective Rashba coupling in the parameter range of interest, we
focus on a system with only Ising SOC; modifications brought by the reintroduction of Rashba are
briefly addressed below. In this case, the
U(2)
K
×
U(2)
K
is reduced to a
U(1)
c
×
U(1)
v
×
U(1)
z,K
×
U(1)
z,K
, where
U(1)
c
and
U(1)
v
denote the charge and valley symmetries while
U(1)
z,K
(
)
rep-
resents spin rotations about the
z
axis in either valley. It naturally follows that the
SU(2)
s
spin
degeneracies present in the absence of SOC are lifted. In fact, because the Ising SOC term resem-
bles a self-generated “spin-valley order”,
τ
z
s
z
, its presence results in a ‘singly polarized’ state even
in the non-interacting limit. (As discussed in SI, section 1 and shown in Extended Data Fig. 6a,
the splitting induced by the bare Ising coupling is substantially smaller than what is required to ex-
plain the phases seen in experiment.) Clearly,
if
the SP or SVP phases were the preferred ground
state without SOC, the mean field Hamiltonian generated in the presence of SOC would have a
clear energetic preference for Ising-like SVP polarized states. Indeed, we consistently find that the
effective Ising SOC is enhanced by the interactions.
The introduction of Rashba SOC breaks the
SU(2)
spin symmetries operative in either valley,
leaving only a
U(1)
c
×
U(1)
v
flavour symmetry. Although it is included below, it has relatively
little qualitative effect on the resulting phase diagrams.
In Extended Data Fig. 6c,e we present the normalized frequencies expected in quantum oscillations
and the corresponding Fermi surfaces obtained through simulations performed with
λ
I
= 0
.
7
meV
and
λ
R
= 3
meV. Before describing these results in detail, we emphasize that c and e were both
simulated using a single choice of
H
init
.
MF
. As described in SI, section 4, although all phases repre-
sented in Extended Data Fig. 6c,e are
low energy
states, it is possible that our algorithm has not
found the
true
ground state of the model. Since the simplicity of the model prevents us from mak-
ing quantitative predictions based on its behaviour, we do not view this as a particularly significant
failing of the simulations. Nevertheless, we have verified by doing multiple runs with different
starting positions that the phases given in Extended Data Fig. 6c,e are not flukes of our specific
choice of
H
init
.
MF
but are instead overall representative of the different low energy states present.
The plot in Extended Data Fig. 6c was obtained in the absence of
U(1)
v
breaking (as explained in
SI, section 4, the presence of
U(1)
v
is enforced by our choice of
H
init
.
MF
). Comparing with Extended
Data Fig. 6a, it is clear even for low dopings, in the
FP(6
,
6)
+
phase, that the splitting between
11