of 8
Nature | www.nature.com |
1
Article
Superconductivity in metallic twisted
bilayer graphene stabilized by WSe
2
Harpreet Singh Arora, Robert Polski, Yiran Zhang, Alex Thomson, Youngjoon Choi, Hyunjin
Kim, Zhong Lin, Ilham Zaky Wilson, Xiaodong Xu, Jiun-Haw Chu, Kenji Watanabe, Takashi
Taniguchi, Jason Alicea & Stevan Nadj-Perge
In the format provided by the
authors and unedited
Supplementary information
https://doi.org/10.1038/s41586-020-2473-8
Nature | www.nature.com/nature
Supplementary Information: Superconductivity in metal-
lic twisted bilayer graphene stabilized by WSe
2
Harpreet Singh Arora
, Robert Polski
, Yiran Zhang
, Alex Thomson, Youngjoon Choi, Hyunjin
Kim, Zhong Lin, Ilham Zaky Wilson, Xiaodong Xu, Jiun-Haw Chu, Kenji Watanabe, Takashi
Taniguchi, Jason Alicea and Stevan Nadj-Perge
I Theoretical Modeling: Spin-orbit coupling in twisted bilayer graphene
Monolayer graphene with induced spin-orbit coupling:
We begin by describing the induced
spin-orbit interaction felt by monolayer graphene adjacent to a transition metal dichalcogenide
(TMD). In the absence of spin-orbit coupling, the low-energy Hamiltonian for the monolayer is
H
0
=
̄
hv
0
k
Ψ
(
k
)
(
k
x
σ
x
τ
z
+
k
y
σ
y
)
Ψ(
k
)
,
(1)
where
Ψ(
k
)
is an eight component spinor with sublattice, valley, and spin indices. The Pauli ma-
trices
σ
x,y,z
act on sublattice indices of the spinor, while
τ
x,y,z
act on the valley indices. The Fermi
velocity is
v
0
=
3
ta/
(2 ̄
h
)
10
6
m/s where
t
= 2
.
61
eV is the nearest-neighbour tunnelling be-
tween carbon atoms and
a
= 0
.
246
nm is the lattice constant
55
. The proximate TMD may induce a
combination of Ising, Rashba, and Kane-Mele spin-orbit terms, which may be included by taking
H
0
H
0
+
H
SO
, where
56, 57
H
SO
=
k
Ψ
(
k
)
(
λ
I
2
τ
z
s
z
+
λ
R
2
(
τ
z
σ
x
s
y
σ
y
s
x
) +
λ
KM
2
τ
z
σ
z
s
z
)
Ψ(
k
)
.
(2)
Here,
s
x,y,z
act on the spin indices. The parameters
λ
I
,
λ
R
, and
λ
KM
respectively quantify
the strength of the Ising, Rashba, and Kane-Mele
58
spin orbit couplings. Their values, as ex-
tracted from density functional theory calculations
56, 59
and various experiments
60–64
, vary widely:
λ
I
1
5
meV, and
λ
R
1
15
meV. While DFT does not predict a substantial
λ
KM
, recent
experiments
65
attributed Kane-Mele SOI as the dominant source of spin relaxation close to the
CNP.
We note that a mass term,
z
/
2
, may also be generated due to inversion symmetry
breaking by the TMD. However, DFT predicts reasonably small values
56, 59
for the mass relative
to Ising and Rashba SOI. Given this expectation, and our focus on the effects of SOI, we hereafter
neglect this contribution.
Twisted bilayer graphene without spin-orbit coupling:
We briefly turn away from the question
of spin orbit to provide a short overview of the continuum model
66, 67
of twisted bilayer graphene
1
with twist angle
θ
. It suffices to specialise to states proximate to the
+
K
valley (
i.e.
,
τ
z
= +1
).
We let
ψ
t
(
ψ
b
) represent the electron annihilation operator of top (bottom) layer. The model Hamil-
tonian may be expressed as
H
cont
=
H
t
+
H
b
+
H
tun
.
(3)
The first two terms on the left-hand side respectively denote the Dirac Hamiltonian of the top and
bottom layers in the absence of interlayer tunnelling:
H
t/b
=
k
ψ
t/b
(
k
)
h
t/b
(
k
)
ψ
t/b
(
k
)
,
(4)
where
h
t
(
k
) =
̄
hv
0
e
iθσ
z
/
4
k
·
σ
e
iθσ
z
/
4
,
h
b
(
k
) =
̄
hv
0
e
iθσ
z
/
4
k
·
σ
e
iθσ
z
/
4
.
(5)
As in Eq. (1),
v
0
is the Fermi velocity of graphene. The layers tunnel through
H
tun
=
`
=1
,
2
,
3
k
ψ
t
(
k
)
T
`
ψ
b
(
k
+
q
`
) +
h.c.,
(6)
where
q
`
=
k
θ
(
sin
[
2
π
3
(
`
1)
]
ˆ
x
+ cos
[
2
π
3
(
`
1)
]
ˆ
y
)
,
k
θ
=
4
π
3
a
2 sin(
θ/
2)
,
(7)
and
T
`
=
w
0
+
w
1
(
e
2
π
(
`
1)
i/
3
σ
+
+
e
2
πi
(
`
1)
/
3
σ
)
.
(8)
2
Twisted bilayer graphene with induced spin-orbit coupling:
Suppose now that a TMD resides
adjacent to one of the graphene monolayers that compose twisted bilayer graphene. We assume
(in correspondence to exprimental realization) that the TMD touches the top layer. The primary
modification to the continuum model presented in Eq. (3) then occurs in
H
t
:
H
t
=
k
ψ
t
(
k
) (
h
t
(
k
) +
h
t,
SO
)
ψ
t
(
k
)
.
(9)
Here,
h
t,
SO
represents the appropriately rotated projection of
H
SO
(Eq. (2)) onto the
+
K
valley:
h
t,
SO
=
e
iθσ
z
/
4
(
λ
I
2
s
z
+
λ
R
2
(
σ
x
s
y
σ
y
s
x
) +
λ
KM
2
σ
z
s
z
)
e
iθσ
z
/
4
.
(10)
Given the interfacial nature of the induced spin-orbit effect, we do not expect substantial modifi-
cations to the Hamiltonian governing the bottom monolayer (when the TMD is only placed on one
side).
II Landau-level spectrum near charge neutrality
We begin by very briefly outlining prior theoretical treatments of the Landau-level spectrum in
BN-encapsulated TBG—for which SOI is not expected to play a role. Specializing to the exper-
iments presented in this paper, we then demonstrate that SOI is not only generically expected to
induce further splittings, but that the simplest model of TBG with SOI can return the observed LL
spectrum (Fig. 3) for a reasonable set of parameters. Despite their importance, interactions are
neglected in the discussion below for the sake of simplicity.
Boron nitride encapsulated TBG:
The Landau-level spectrum observed near the CNP in TBG
encapsulated in BN sufficiently far from the magic angle is usually 8-fold symmetric
68, 69
. This
result is in agreement with the band structure one obtains from the continuum model. In the
absence of both SOI and a magnetic field, the system possesses two Dirac cones per
K
-valley per
spin at the
+
κ
and
κ
points in the moir
́
e Brillouin zone—thus implying eight degenerate Fermi
surfaces upon either electron or hole doping (see Fig. 4 and Extended Data Fig. 10). Equation (1)
serves as a good low-energy description of the system close to the CNP provided a few changes are
made: (i) the Pauli matrices
σ
x,y,z
now act on the band indices of the operators
Ψ
(
τ
x,y,z
and
s
x,y,z
still act on the
K
and spin flavour indices, respectively); (ii) there is an additional flavour index
given by the moir
́
e valley,
±
κ
; and (iii) the Fermi velocity is substantially renormalized from its
value in monolayer graphene.
3
Experimentally, however, this simple band-structure picture appears insufficient close to the
magic angle: the 8-fold degeneracy is lifted to a 4-fold degeneracy
70, 71
. The latter degeneracy
is instead obtained when a
single
Fermi surface is present per spin per (
K
-)valley upon doping
away from the CNP. Such a scenario has been shown to occur for reasonable parameters close
to the magic angle, where the band structure is very sensitive to small changes in the angle and
tunnelling parameters
(
w
0
,w
1
)
72
. Alternatively, in the presence of strain, the Dirac cones are no
longer pinned to
±
κ
and move closer together—again resulting in a single Fermi surface per spin
and
K
-valley upon doping away from charge neutrality
73, 74
. Each Fermi surface yields a sequence
±
1
,
±
2
,
±
3
,...
, which becomes
±
4
,
±
8
,
±
12
,...
upon accounting for the spin and valley degrees
of freedom. The definitive explanation for the observed reduction in degeneracy, however, remains
unclear at present.
Landau levels for TBG with SOI:
Spin-orbit coupling breaks the
SU(2)
spin symmetry and
provides another mechanism for reducing the Landau-level degeneracy. Starting from the sce-
narios highlighted in the previous section, the usual
±
4
,
±
8
,
±
12
,...
sequence observed near the
magic angle would then naturally further split down to
±
2
,
±
4
,
±
6
,...
in the presence of SOI.
We stress, however, that effects such as strain that can return fourfold degenerate Fermi surfaces
close to charge neutrality are not included in the band structures shown in Fig. 4 and Extended
Data Fig. 10. In these figures, two Fermi surfaces enclosing
κ
per spin per
K
-valley are present
at charge neutrality, as in the scenario discussed above for twisted bilayer graphene away from
the magic angle. We show here that the experimentally observed twofold symmetry of Landau
levels may nevertheless be obtained directly in the limit that the low-energy theory is described by
eight Dirac cones with SOI. That is, SOI can break the the eighfold degeneracy down to a twofold
degeneracy without have to take any other effects (e.g., strain) into account.
We first assume that we are close enough to charge neutrality that the Dirac description
of TBG remains valid. As outlined in the previous section, an appropriately repurposed Eq. (1)
provides a good starting point. SOI may then be included in a manner directly analogous to the
approach we took above with monolayer graphene (
i.e.
,
H
SO
defined in Eq. (2)). We note, however,
that the spin-orbit parameters (
λ
I
,
λ
R
,
λ
KM
) are not prohibited from taking different values at the
±
κ
Dirac cones (the parameters at
±
K
are related by time reversal). Nevertheless, we assume for
simplicity that any differences in the effective parameters between
κ
-valleys is unimportant.
Without loss of generality, we therefore begin by considering the Dirac cone at
+
κ
with
valley flavour
+
K
(which, according to the assumption made above, is equally valid for the Dirac
cones at
κ
with valley flavour
K
). The first-quantized Hamiltonian takes the form
H
=
̄
hv
F
(
i∂
x
σ
x
+
i∂
y
σ
y
) +
̃
λ
I
2
s
z
+
̃
λ
R
2
(
σ
x
s
y
σ
y
s
x
) +
̃
λ
KM
2
σ
z
s
z
.
(11)
This effective description may be derived in perturbation theory in powers of the interlayer tun-
nelling (
w
1
/
2
/
( ̄
hv
0
k
θ
)
, to be precise) in a manner analogous to Ref. 67’s original derivation of the
magic angle. We have added a tilde to the SOI parameters as compared with previous sections to
4
emphasise that
̃
λ
I
,
̃
λ
R
, and
̃
λ
KM
describe the
effective
SOI relevant to the moir
́
e Dirac cones, and
that they are not expected to be the same as the SOI induced by the TMD on monolayer graphene.
We have again neglected a possible mass term.
The magnetic field may now be included in a straightforward manner by taking
i
i
+
e
A
, where
e
is the electronic charge and
A
is the vector potential corresponding to a
magnetic field of strength
B
(the Zeeman splitting is negligible at the energy scales considered).
We work in the Landau gauge, setting
A
=
B
(
y,
0)
and write the wavefunction as
Φ(
x,y
) =
e
ikx
φ
(
y
)
,
φ
(
y
) = (
φ
1
(
y
)
1
(
y
)
2
(
y
)
2
(
y
))
T
where
1
,
2
denote the band indices (acted on by
σ
x,y,z
) and
,
denote the spin (acted on by
s
x,y,z
). The eigenvalue equation may be expressed as
55
(
y
) =
̃
λ
+
/
2
0
̄
c
A
0
0
̃
λ
+
/
2
i
̃
λ
R
̄
c
A
̄
c
A
i
̃
λ
R
̃
λ
/
2
0
0
̄
c
A
0
̃
λ
/
2
φ
(
y
)
(12)
where
̃
λ
±
= (
̃
λ
KM
±
̃
λ
I
)
/
2
,
E
is the energy to be obtained, and
ω
c
=
2
v
F
/`
B
is the cyclotron
frequency, with
`
B
=
̄
h/eB
the magnetic length. The operator
A
=
`
B
y
+
`
B
k
y/`
B
is an annihilation operator, while
A
=
`
B
y
+
`
B
k
y/`
B
is a creation operator. They satisfy
[
A,A
] = 1
. It follows that the functions
φ
(
y
)
are superpositions of solutions to the 1
d
simple
harmonic oscillator. Further, each solution
φ
(
y
)
with energy
E
corresponds to an extensively large
set of degenerate states labelled by
k
. With this information, it is straightforward to show that there
are generically four distinct solutions satisfying
0 =
E
̃
λ
,
(13)
0 =
(
E
+
̃
λ
+
2
)(
E
+
̃
λ
2
)(
E
̃
λ
2
)
̃
λ
2
R
(
E
̃
λ
2
)
ω
2
c
(
E
+
̃
λ
2
)
.
We denote these four solutions by
E
(+)
1
,
μ
= 1
4
. The remaining energies are obtained as
solutions to
0 =
(
(
n
1)
ω
2
c
(
E
̃
λ
+
2
)(
E
+
̃
λ
2
))(
2
c
(
E
+
̃
λ
+
2
)(
E
̃
λ
2
))
(14)
λ
2
R
(
E
̃
λ
+
2
)(
E
̃
λ
2
)
for
n
2
. We similarly label these energies by
E
(+)
n,μ
,
μ
= 1
4
. It is then straightforward to show
that the Landau level spectrum corresponding to a Dirac cone of the opposite chirality (
τ
z
=
1
)
may be obtained directly by solving Eqs. (13) and (14) but with
̃
λ
I
→ −
̃
λ
I
(
i.e.
̃
λ
±
̃
λ
). We
denote these energies
E
(
)
n,μ
.
We present in Extended Data Fig. 9b and Extended Data Fig. 9d solutions to Eqs. (13)
and (14) for
n
30
. Since the density of a Dirac cone (without a magnetic field) is proportional to
5
the energy squared, we plot the magnetic field against the energy
squared
. The parameters chosen
are
(
̃
λ
I
,
̃
λ
R
,
̃
λ
KM
) = (3
,
4
,
0)
meV and
(
̃
λ
I
,
̃
λ
R
,
̃
λ
KM
) = (1
.
5
,
2
.
5
,
2)
meV with
v
F
10
5
m/s, as
appropriate for
θ
= 0
.
79
and
θ
= 0
.
87
.
In the physical sample, the Landau levels are not actually expected to be infinitely degenerate,
but rather possess some finite width. We model this effect through a phenomenological broadening
of the density of states
D
(
E,B
)
:
D
(
E,B
) = 2
·
eB
2
π
n
=1
4
μ
=1
η
=
±
1
2
π
Γ
exp
[
(
E
(
η
)
n,μ
E
)
2
2
]
.
(15)
The coefficient out front,
2
·
eB/
(2
π
)
, gives the degeneracy of each Landau level—
eB/
(2
π
)
times
the number of Dirac cones within each of the
K
valleys. The sum over
η
=
±
then accounts for
the two valley flavours
±
K
. We present colour density plots of this function in Extended Data
Fig. 9a and Extended Data Fig. 9c as a function of the magnetic field and
E
2
using the same SOI
parameters as above with broadenings
Γ = 0
.
22
meV and
Γ = 0
.
15
meV, respectively. The sum
over
n
is evaluated up to
n
= 30
. Notably, this simple model is able to reproduce several key
features of the Landau level sequence shown in Fig. 3. For instance, at the higher fields shown, a
sequence 2, 4, 6, 8, 10,. . . is observed. By contrast, at lower fields, we instead find 4, 6, 10, 14,. . . .
Supplementary References:
55. Neto, A. C., Guinea, F., Peres, N. M., Novoselov, K. S. & Geim, A. K. The electronic proper-
ties of graphene.
Rev. of Mod. Phys.
81
, 109–162 (2009).
56. Gmitra, M. & Fabian, J. Graphene on transition-metal dichalcogenides: A platform for prox-
imity spin-orbit physics and optospintronics.
Phys. Rev. B
92
, 155403 (2015).
57. Zaletel, M. P. & Khoo, J. Y. The gate-tunable strong and fragile topology of multilayer-
graphene on a transition metal dichalcogenide. Preprint at https://arxiv.org/abs/1901.01294
(2019).
58. Kane, C. L. & Mele, E. J. Quantum spin Hall effect in graphene.
Phys. Rev. Lett.
95
, 226801
(2005).
59. Gmitra, M., Kochan, D., H
̈
ogl, P. & Fabian, J. Trivial and inverted Dirac bands and the
emergence of quantum spin Hall states in graphene on transition-metal dichalcogenides.
Phys.
Rev. B
93
, 155104 (2016).
60. Zihlmann, S.
et al.
Large spin relaxation anisotropy and valley-Zeeman spin-orbit coupling in
WSe
2
/graphene/
h
-BN heterostructures.
Phys. Rev. B
97
, 075434 (2018).
61. Island, J. O.
et al.
Spin-orbit-driven band inversion in bilayer graphene by the van der Waals
proximity effect.
Nature
571
, 85–89 (2019).
6
62. Wang, D.
et al.
Quantum Hall effect measurement of spin-orbit coupling strengths in ultraclean
bilayer graphene/WSe
2
heterostructures.
Nano Lett.
19
, 7028–7034 (2019).
63. Wang, Z.
et al.
Origin and magnitude of ‘designer’ spin-orbit interaction in graphene on
semiconducting transition metal dichalcogenides.
Phys. Rev. X
6
, 041020 (2016).
64. Cummings, A. W., Garcia, J. H., Fabian, J. & Roche, S. Giant spin lifetime anisotropy in
graphene induced by proximity effects.
Phys. Rev. Lett.
119
, 206601 (2017).
65. Wakamura, T.
et al.
Spin-orbit interaction induced in graphene by transition metal dichalco-
genides.
Phys. Rev. B
99
, 245402 (2019).
66. Sichau, J.
et al.
Resonance microwave measurements of an intrinsic spin-orbit coupling gap
in graphene: A Possible indication of a topological state.
Phys. Rev. Lett.
122
, 046403 (2019).
67. Bistritzer, R. & MacDonald, A. H. Moir
́
e bands in twisted double-layer graphene.
Proc. Natl.
Acad. Sci. USA
108
, 12233–12237 (2011).
68. Sanchez-Yamagishi, J. D.
et al.
Quantum Hall effect, screening, and layer-polarized insulating
states in twisted bilayer graphene.
Phys. Rev. Lett.
108
, 076601 (2012).
69. Chung, T.-F., Xu, Y. & Chen, Y. P. Transport measurements in twisted bilayer graphene:
Electron-phonon coupling and Landau level crossing.
Phys. Rev. B
98
, 035425 (2018).
70. Cao, Y.
et al.
Correlated insulator behaviour at half-filling in magic-angle graphene superlat-
tices.
Nature
556
, 80–84 (2018).
71. Yankowitz, M.
et al.
Tuning superconductivity in twisted bilayer graphene.
Science
363
,
1059–1064 (2019).
72. Hejazi, K., Liu, C. & Balents, L. Landau levels in twisted bilayer graphene and semiclassical
orbits.
Phys. Rev. B
100
, 035115 (2019).
73. Bi, Z., Yuan, N. F. Q. & Fu, L. Designing flat bands by strain.
Phys. Rev. B
100
, 035448
(2019).
74. Zhang, Y.-H., Po, H. C. & Senthil, T. Landau level degeneracy in twisted bilayer graphene:
Role of symmetry breaking.
Phys. Rev. B
100
, 125104 (2019).
7