PHYSICAL REVIEW B
108
, 195148 (2023)
Electrostatic fate of
N
-layer moiré graphene
Kryštof Kolá
ˇ
r,
1
,
*
Yiran Zhang,
2
,
3
,
4
Stevan Nadj-Perge,
2
,
3
Felix von Oppen,
1
and Cyprian Lewandowski
5
,
6
1
Dahlem Center for Complex Quantum Systems and Fachbereich Physik, Freie Universität Berlin, 14195 Berlin, Germany
2
T. J. Watson Laboratory of Applied Physics, California Institute of Technology, 1200 East California Boulevard,
Pasadena, California 91125, USA
3
Institute for Quantum Information and Matter, California Institute of Technology, Pasadena, California 91125, USA
4
Department of Physics, California Institute of Technology, Pasadena, California 91125, USA
5
National High Magnetic Field Laboratory, Tallahassee, Florida 32310, USA
6
Department of Physics, Florida State University, Tallahassee, Florida 32306, USA
(Received 21 July 2023; revised 28 October 2023; accepted 7 November 2023; published 27 November 2023)
Twisted
N
-layer graphene (TNG) moiré structures have recently been shown to exhibit robust superconduc-
tivity similar to twisted bilayer graphene (TBG). In particular for
N
=
4and
N
=
5, the phase diagram features
a superconducting pocket that extends beyond the nominal full filling of the flat band. These observations
are seemingly at odds with the canonical understanding of the low-energy theory of TNG, wherein the TNG
Hamiltonian consists of one flat-band sector and accompanying dispersive bands. Using a self-consistent
Hartree-Fock treatment, we explain how the phenomenology of TNG can be understood through an interplay of
in-plane Hartree and inhomogeneous layer potentials, which cause a reshuffling of electronic bands. We extend
our understanding beyond the case of
N
=
5 realized in experiment so far. We describe how the Hartree and layer
potentials control the phase diagram for devices with
N
>
5 and tend to preclude exchange-driven correlated
phenomena in this limit. To circumvent these electrostatic constraints, we propose a flat-band paradigm that
could be realized in large-
N
devices by taking advantage of two nearly flat sectors acting together to enhance the
importance of exchange effects.
DOI:
10.1103/PhysRevB.108.195148
I. INTRODUCTION
As theoretically predicted in seminal papers [
1
,
2
],
magic-angle twisted bilayer graphene (TBG) comprises
eight flat bands near charge neutrality, two per each spin and
valley flavor. At specific “magic” twist angles, these bands
become maximally flat. The first experimental realization
of this system [
3
,
4
] revealed an intriguing phase diagram as
a function of filling, featuring superconducting domes and
correlated insulating behavior at integer filling reminiscent
of the cuprates. Further experiments confirmed and extended
these findings [
5
–
13
], but also pointed to a strong device
dependence of the phase diagram. An important feature
that emerged as rather robust is the presence of “resets”
near integer fillings, corresponding to Stoner-like flavor
polarization phase transitions [
14
,
15
].
In searching for other promising moiré systems beyond
TBG, Refs. [
16
,
17
] suggested to stack more than two
graphene layers with alternating twist angles, see Fig.
1(a)
.
For any number of layers, this procedure preserves the moiré
translation symmetry of TBG. The resulting band structure
consists of a set of twisted-bilayer-graphene-like bands, re-
ferred to as sectors, at different effective twist angles as
illustrated in Fig.
1(c)
[
16
–
18
]. These sectors (labeled by
k
) are characterized by different spatial charge distributions
over the layers as sketched in Fig.
1(b)
. At the magic
*
Corresponding author: kolar@zedat.fu-berlin.de
angle of the multilayer system, one sector’s effective an-
gle is just the magic angle of twisted bilayer graphene.
Importantly, this procedure allows for the construction of
flat-band devices at larger physical twist angles than for
TBG, yielding a smaller unit cell and potentially more robust
devices.
The proposal of Refs. [
16
,
17
] was first realized for the
case of three layers [
6
,
19
,
20
], which host a set of flat bands
coexisting with a graphene-like Dirac cone. These twisted
trilayer graphene (TTG) devices exhibited strongly-coupled
superconductivity tunable by an out-of-plane displacement
field. More recently, twisted
N
-layer graphene (TNG) devices
have been successfully fabricated for up to
N
=
5 layers
[
21
,
22
]. As for twisted bilayer graphene, the phase diagrams
of the multilayer systems featured prominent superconducting
domes in the temperature-doping plane. However, while in
TBG superconductivity typically terminates at filling
|
ν
|=
3,
Ref. [
22
] observed that for
N
=
5, superconductivity persists
up to a filling of five electrons per unit cell. Strikingly, such
an extended superconducting pocket is at odds with a picture
of decoupled sectors [
16
–
18
].
The extended superconducting pockets reported in
Refs. [
21
,
22
] motivate us to study TNG systems theoreti-
cally. We employ mean-field theory for analytical estimates
and perform fully self-consistent Hartree-Fock calculations.
While mean-field theory has inherent drawbacks and is an
approximate technique, it has proven remarkably successful
in the study of TBG [
23
–
27
], with certain models of TBG
exhibiting Slater-determinant ground states at integer filling
2469-9950/2023/108(19)/195148(22)
195148-1
©2023 American Physical Society
KRYŠTOF KOLÁ
ˇ
R
et al.
PHYSICAL REVIEW B
108
, 195148 (2023)
FIG. 1. (a) Device schematics. We consider
N
-layer graphene
with alternating twist angles in a double-gated setup. Here
θ
is the
physical twist angle. (b) Schematics of layer charge distribution [see
Eq. (
8
)] for
N
=
11, showing the three sectors with lowest effective
twist angle,
k
=
1: red,
k
=
2: orange,
k
=
3: green. In experiment
to date,
k
=
1 is the flattest, “magic” sector. (c) Single-particle band
structure for
N
=
11. (d) Band filling of the magic sector at different
total gate densities. 25
×
10
12
cm
−
2
is the threshold of dielectric
breakdown in current hBN-based samples [
29
].
[
28
]. For a larger number of layers, the Hartree-Fock approx-
imation accounts well for the screening of classical charge
distributions, which we will argue play a crucial role in the
physics of TNG.
We focus first on understanding and explaining the experi-
mental results of Refs. [
21
,
22
]for
N
=
3
,
4
,
5 layers, and then
apply the developed framework to study the twisted
N
-layer
problem for larger values of
N
, characterizing its electronic
properties. We will argue that both in- and out-of-plane elec-
trostatics play a crucial role in shaping the phase diagrams of
TNG systems, providing a simple picture in terms of sector
shifts in Sec.
II
. For larger layer numbers, we show that elec-
trostatically doping the moiré system requires a larger charge
density on the metallic gates. This behavior arises because it is
necessary for the metallic gates to compensate for the charge
redistribution due to interactions. This effect makes it increas-
ingly prohibitive to electrostatically dope
N
>
5 multilayer
structures into the regime where the magic flat band is opti-
mally filled for superconductivity. This is shown in Fig.
1(d)
,
where we also indicate estimates for the gate charge densities.
Interestingly, we find that while going beyond
N
=
5 layers to
study interaction effects of the
k
=
1 flat band presents little
advantage, focusing on the second-harmonic bands (
k
=
2) in
Fig.
1(b)
for
N
5 overcomes the prohibitive electrostatic
barrier and yields very flat bands conducive to correlation
effects.
Our paper is structured as follows. Section
II
presents
a summary of our results focusing on physical under-
standing and experimental trends. Section
III
outlines the
formal description of the
N
-layer problem and introduces
the Hartree-Fock machinery, emphasizing the similarities and
differences with twisted bilayer graphene (TBG). In Sec.
IV
,
we combine physical understanding with Hartree-Fock calcu-
lations for
N
=
3
,
4
,
5, focusing on explaining experimental
trends. Section
V
discusses the electronic properties of
N
>
5 devices in more detail. We conclude with a summary
and discussion in Sec.
VI
. Readers uninterested in details
of the mathematical description can focus on Secs.
II
,
V
,
and
VI
.
II. PHYSICAL UNDERSTANDING AND SUMMARY
OF MAIN RESULTS
A. Experimental motivation
A key physical effect in twisted alternating-angle graphene
multilayers is the cascade of “resets” close to integer fillings
of the flat bands. The resets already occur at relatively high
temperatures, well above those required for the correlated
superconducting and insulating states, and are deduced from
measurements of the chemical potential [
14
,
15
]aswellasthe
Hall conductivity [
13
,
19
]. The cascade of transitions can be
explained in different ways [
14
,
15
,
26
,
30
–
33
], with Ref. [
14
]
interpreting it as Stoner-like flavor (spin and valley) polar-
ization. Within this picture, flat-band superconductivity is
unlikely to exist when three of the four flavors are fully occu-
pied and time-reversed partners at the Fermi level are absent.
In TBG, this happens beyond
ν
=±
3 (see further discussion
in Sec.
III
regarding intervalley coherent orders). Irrespective
of the detailed theoretical symmetry-breaking mechanism,
this expectation is in line with experimental trends. In TBG, a
cascade transition near
ν
=±
3 typically serves as an upper
filling bound for superconductivity [
6
,
13
,
19
,
34
]. Similarly,
a lower filling bound for superconductivity is the cascade
transition at
ν
=±
2.
Cascade phenomenology has also been reported for TNG
systems with
N
=
3
,
4
,
5 layers [
6
,
19
,
21
,
22
], c.f., Fig.
2(a)
.
The band structure of TNG decomposes into decoupled
sectors of TBG-like and (for
N
odd) monolayer-graphene
(MLG)-like bands. This is illustrated by the example band
structure for
N
=
11 in Fig.
1(c)
, which features five TBG-like
bands and a Dirac cone. When one of the TBG-like sectors
is effectively at the magic angle, the cascade features can be
understood as occurring in the magic sector, with the other
sectors being filled uniformly [
21
,
22
].
Startlingly, as shown in Fig.
2(b)
and reported in
Refs. [
21
,
22
], superconductivity persists to higher total fill-
ings (
ν
total
) in TQG (twisted quadrilayer graphene,
N
=
4)
and TPG (twisted pentalayer graphene,
N
=
5), extending up
195148-2
ELECTROSTATIC FATE OF
N
-LAYER MOIRÉ ...
PHYSICAL REVIEW B
108
, 195148 (2023)
0
2
4
6
8
ν
total
0
2
4
6
8
n
H
[10
12
cm
−
2
]
(a)
N
=3
N
=4
N
=5
0
2
4
6
8
ν
total
0
.
0
0
.
5
1
.
0
1
.
5
2
.
0
T
C
[
K
]
(b)
N
=3
N
=4
N
=5
2
4
6
8
10
12
14
⊥
18
16
14
12
10
8
6
(c)
3
4
5
6
7
8
ν
total
K
M
K
M
M
Γ
-40
-20
0
20
40
60
80
E [meV]
(d)
HFX:
BS
0
2
4
DOS [AU]
K
M
K
M
M
Γ
(e)
XFL:
BS
0
2
4
DOS [AU]
K
M
K
M
M
Γ
(f)
HFL:
BS
0
2
4
DOS [AU]
0
2
4
ν
magic
(g)
(h)
N
=3
N
=4
N
=5
(i)
0
2
4
6
8
ν
total
0
1
ν
flavor
0
2
4
6
8
ν
total
0
2
4
6
8
ν
total
FIG. 2. (a) Experimental data of the Hall density vs total filling showing the the cascade transitions (arrows) for
N
=
3
,
4
,
5, see Ref. [
22
]
for details on the samples and measurements. (b) Corresponding experimental data for
T
C
domes for
N
=
3
,
4
,
5. (c) Colormap of
ν
total
needed
to reach filling
ν
magic
=
3 of the magic sector in the
ε
-
ε
⊥
plane. (d) Interacting structures and densities of states for TPG at
ν
magic
≈
4, including
in-plane Hartree and Fock (HFX) terms from Eq. (
3
). Shown are the
k
=
1 magic sector (red) and
k
=
2 nonmagic TBG-like sector (blue).
Dashed red lines denote the location of the Fermi level. (e) Same as (d) but with layer Hartree potentials and Fock (XFL). (f) Same as (d) but
including all terms, that is, HFL. (g) Total magic-sector filling (top) and flavor-resolved magic filling (bottom), showing the cascade with
in-plane Hartree and Fock (HFX) for
N
=
3
,
4
,
5. (h) Same as (g), but with out-plane Hartree and Fock (XFL) (i) Same as (g), but including
all the terms (HFL).
to
ν
total
=
5forthe
N
=
5 case of TPG. Simultaneously the
cascade “resets” also set in at higher filling fractions.
Assuming that doping of the magic sector (
ν
magic
)isin
the optimal range for superconductivity, i.e., approximately
2–3 electrons per moiré cell, these observations would im-
ply substantial filling of the nonmagic sectors at odds with
a simple band-structure picture. The nonmagic sectors are
strongly dispersive, so that their noninteracting band structure
would predict almost no filling. Specifically, complete filling
of the magic bands (
ν
magic
=
4) would be accompanied by
a filling of less than
≈
0
.
06 electrons per moiré unit cell in
the nonmagic bands for TPG and less than
≈
0
.
02 electrons
per moiré unit cell for TQG. (We measure fillings relative to
charge neutrality.)
To obtain these estimates, we note that the
N
/
2
TBG-like
electronic sectors (
k
=
1
,
2
,...,
N
/
2
)haveeffectivetwist
angles [
16
]
θ
eff
k
=
θ
2 cos
π
k
N
+
1
,
(1)
which differ from the physical twist angle
θ
,Fig.
1(a)
.This
formula reveals the advantage of multilayers—one can obtain
a sector effectively at the magic angle for devices at a larger
twist angle
θ
. This is best exploited by choosing the
k
=
1
sector to lie effectively at the magic angle, which maximizes
the physical twist angle. All the current experiments on mul-
tilayer alternating twist angle systems make this choice, and
we shall also make it our default choice for analysis. However,
we note that for large
N
, the choice
k
magic
=
2 also becomes
feasible. We will return to this possibility in Sec.
V
. Approx-
imating the nonmagic sectors as Dirac cones, their filling is
195148-3
KRYŠTOF KOLÁ
ˇ
R
et al.
PHYSICAL REVIEW B
108
, 195148 (2023)
(see Appendix
A4
)
ν
non-magic
=
k
∈
nonmagic
ν
k
≈
k
∈
nonmagic
A
uc
N
f
c
k
4
π
( ̄
h
v
(
k
)
D
)
2
μ
2
k
.
(2)
Here,
c
k
=
2(
c
k
=
1) if the sector
k
is TBG-like (MLG-like),
μ
k
is the effective chemical potential in sector
k
,
N
f
=
4isthe
number of flavors, and
v
(
k
)
D
is the Dirac velocity in sector
k
.In
the absence of interactions,
μ
k
=
μ
magic
with
μ
magic
the magic
sector Fermi energy. A filled magic sector corresponds to
μ
magic
≈
W
/
2, where
W
is the noninteracting bandwidth. This
bandwidth varies with strain, taking values
W
20 meV.
Even at the upper limit for
W
, we then find
ν
non-magic
0
.
06
for TPG (using
v
(
k
=
2)
D
≈
0
.
35
v
D
). For TQG, the
k
=
2 sector
has an even larger detuning from the magic angle (
θ
eff
k
=
2
=
2
.
9
◦
), so that
v
(
k
=
2)
D
≈
0
.
6
v
D
and
ν
non-magic
0
.
02. Therefore,
the enhanced nonmagic-sector filling [
21
,
22
] is an interaction
effect, motivating our Hartree-Fock study of TNG.
B. Physical understanding
Electron-electron interactions alter the above considera-
tions predominantly through two terms in the Hamiltonian, as
can be seen by examining the mean-field decomposition (see
Secs.
III B
and
III C
for details)
H
MF
=
H
SP
+
H
Hartree
+
H
Fock
+
H
layer
.
(3)
First, interactions represented by the Hartree and Fock mean-
field terms broaden the noninteracting magic bands, promote
the onset of symmetry-breaking order, and, crucially for
our analysis, induce filling-dependent upward shifts of the
quasiparticle energies relative to nonmagic sectors. This
Hartree-dominated shift arises because the electron density of
the TBG-like sectors is spatially inhomogeneous in the 2D
plane, which is associated with a cost in Coulomb energy
[
26
,
35
–
39
]. Importantly, the inhomogeneity is particularly
strong in the magic sector and decreases with detuning from
the magic angle. Second, the contribution H
layer
is new to
N
>
2 layers and arises because the sectors have different
vertical charge distributions across layers [
21
,
22
]asshown
in Fig.
1(b)
. These distributions are given by the layer de-
pendence of the wave functions, taking the form of standing
waves analogous to a particle-in-a-box problem. The sector
with lowest effective twist angle,
k
=
1, corresponds to the
first harmonic, which is singly peaked at the center of the
stack. The
k
=
2 sector is the second harmonic with a doubly-
peaked structure, and so on. The different layer-dependent
charge distributions imply that the sectors have different ener-
gies due to the electric potential produced by the gate charges.
For the devices investigated experimentally (magic sector
k
=
1), both
H
Hartree
and
H
layer
effects enhance the occupa-
tion of the nonmagic sectors relative to the noninteracting
band-structure scenario described above. The first mechanism
postpones the occupation of the magic sector as it is broad-
ened and shifted upward in energy as it is filled. A similar
shift in energy also occurs for the second mechanism. The
potential produced by the gate charges in combination with
the induced charges in TNG has a maximum in the central
layer. (Note that in the absence of a displacement field, the
electric field vanishes at the center by symmetry. Moreover,
TABLE I. Inverse capacitance (
C
−
1
)
k
,
k
for
k
,
k
∈{
1
,
2
}
, evalu-
ated for layer numbers
N
=
4,
N
=
5, and
N
→∞
.
N
(
C
−
1
)
1
,
1
(
C
−
1
)
1
,
2
=
(
C
−
1
)
2
,
1
(
C
−
1
)
2
,
2
4
0.262
0.1
0.0382
5
0.403
0.208
0.125
N
→∞
0.147
N
0.115
N
0.099
N
the potential drops towards the, say, positively charged gate
electrodes above and below the TNG stack.) Due to this po-
tential maximum, the energy is higher for sectors, in which
charge is more localized near the central layer. Thus, this
mechanism also predicts that the magic sector is pushed up
in energy relative to the nonmagic sectors.
We can provide an estimate of this electrostatically induced
band shifting, which will be verified in later sections through
extensive Hartree-Fock calculations. First we assume that in
the presence of interactions, the overall band structure of each
sector remains fixed (i.e., given by the noninteracting band
structure) and only the chemical potential of each sector
μ
k
shifts as
μ
k
=
μ
−
U
k
−
G
k
.
(4)
Here,
U
k
and
G
k
quantify the shifts due to
H
layer
and H
Hartree
,
and
μ
is the chemical potential of the whole system. We take
G
k
=
0 for all sectors except the magic sector (
k
=
1) as it
has the largest in-plane inhomogeneity (see Appendix
A2
). In
the magic sector [
26
,
35
,
36
,
40
],
G
k
∼
e
2
/
(4
πε
ε
0
L
M
), where
L
M
is the moiré period. For TPG, depending on dielectric
constant,
G
k
can be as large as 30 meV, giving a filling of
up to
ν
non-magic
≈
1
.
1 of the nonmagic sectors at full filling of
the magic sector,
ν
magic
=
4. This should be compared to the
noninteracting estimate of
ν
non-magic
≈
0
.
06 given above.
Inclusion of the shift
U
k
induced by
H
layer
can further
increase the filling of
ν
non-magic
.ThetermH
layer
contributes
nontrivially due to imperfect screening of the gate electrodes
by the layers and becomes increasingly important as
N
grows.
The energy shift
U
k
of a sector
k
can be approximated by (see
Appendix
C2
)
U
k
=
e
2
d
l
A
uc
ε
0
ε
⊥
k
(
C
−
1
)
k
,
k
ν
k
(5)
for given sector fillings
ν
k
. Here,
A
uc
is the unit-cell area,
d
l
is
the layer distance, and
ε
⊥
the out-of-plane dielectric constant
of the graphene layers. The matrix
C
in sector space is a
dimensionless capacitance-like matrix, which we tabulate for
TQG, TPG, as well as large
N
in Table
I
(see Appendix
C2
for formulas for arbitrary
N
and derivations).
For TPG we obtain a shift of up to 45 meV, allowing for
a filling of up to
ν
non-magic
≈
2
.
5 for full filling of the magic
sector,
ν
magic
=
4 (see Sec.
III B
for further discussion).
Modifications of the quasiparticle dispersion by H
Fock
would tend to reduce the above estimates. However, we also
highlight that the effects of the layer potential H
layer
and the
Hartree correction H
Hartree
mutually reinforce each other. To
illustrate this, consider
U
1
−
U
2
=
10 meV and
G
1
=
10 meV
and small bandwidth
W
/
2
=
2 meV. Taken separately, each
term would only yield a tiny
ν
non-magic
≈
0
.
07. On the other
195148-4
ELECTROSTATIC FATE OF
N
-LAYER MOIRÉ ...
PHYSICAL REVIEW B
108
, 195148 (2023)
hand, taking
μ
2
=
22 meV in Eq. (
2
) yields a four times larger
ν
non-magic
≈
0
.
3. This highlights the importance of considering
both shift mechanisms simultaneously.
To conclude the qualitative analysis of this section, we
comment on the relative role the three interaction terms in
Eq. (
3
) play as the layer number
N
increases. The inverse
capacitance matrix (
C
−
1
)
k
,
k
is a decreasing function of
k
and
k
. Physically, larger-
k
sectors screen the gate field better,
therefore generating smaller layer potentials. This monotonic
decrease implies that
U
magic
−
U
k
>
0 for any (nonmagic)
k
>
1. Thus, the effective chemical potentials
μ
k
of the nonmagic
sectors increase, enhancing their occupations. Secondly, for
fixed
k
and
k
,(
C
−
1
)
k
,
k
scales linearly with the vertical ex-
tent (as the inverse capacitance of a parallel-plate capacitor)
and thus with the number of layers
N
. This suggests that
the layer potential grows in importance with
N
, eventually
dominating over other contributions for large
N
. Indeed, other
contributions to the mean-field Hamiltonian do not grow with
the number of layers. This suggests that the layer potentials
become dominant at large
N
and doping of the central
k
=
1
sector by gating will be preempted by dielectric breakdown
[
29
], as shown in Fig.
1(d)
. We return to this analysis using
Hartree-Fock calculations in subsequent sections.
III. MODEL
In this section, we introduce the noninteracting model,
specify the interaction, and discuss the mean-field decoupling.
While we largely follow standard procedures for the mean-
field description of moiré graphene [
24
–
27
,
41
–
43
], we allow
for layer dependence of the interaction and include the layer
potential term that is usually ignored.
A. Twisted graphene multilayers
We consider
N
-layer alternating angle twisted graphene.
Focusing on the K-valley, the single-particle Bistritzer-
MacDonald Hamiltonian reads [
16
]
H
K
sp
=
⎛
⎜
⎜
⎜
⎜
⎜
⎜
⎜
⎝
h
−
θ/
2
(
k
)
T
†
(
r
)0
···
0
T
(
r
)
h
θ/
2
(
k
)
T
(
r
)
0
T
†
(
r
)
h
−
θ/
2
(
k
)
.
.
.
.
.
.
0
h
(
−
1)
N
θ/
2
(
k
)
⎞
⎟
⎟
⎟
⎟
⎟
⎟
⎟
⎠
,
(6)
where
h
θ/
2
(
k
)
=−
i
̄
h
v
D
(
σ
·
k
)
e
i
θσ
z
denotes the Dirac Hamil-
tonians of the layers (our numerics neglects the rotation
of the Dirac terms) and
T
(
r
)
=
2
j
=
0
T
j
e
i
q
j
·
r
is the in-
terlayer hopping with
T
j
=
w
AA
σ
0
+
w
AB
[
σ
x
cos(2
π
j
/
3)
+
σ
y
sin(2
π
j
/
3)] and
q
j
=
(
O
3
)
j
(
K
2
−
K
1
)
=
2
|
K
|
sin(
θ/
2)
(
O
3
)
j
[0
,
−
1] with
K
i
the Dirac-point positions in layer
i
and
O
3
the matrix of a counterclockwise 120
◦
rotation. Neglecting
possible layer dependence, we account for lattice relaxation
by choosing
w
AA
=
80 meV
,
w
AB
=
110 meV [
18
]. Disper-
sion and Bloch wave functions of the
K
valley follow by
time-reversal symmetry. The model of Eq. (
6
)isamini-
mal description of
N
-layer systems, neglecting relative layer
displacements [
16
,
44
], next-nearest-layer hoppings [
16
], pe-
riodic strain [
45
], and layer dependence of lattice corrugation
[
18
]. While these additional ingredients modify the quantita-
tive details of the electronic spectrum, they do not alter the
two key features, namely the inhomogeneous charge distribu-
tion and the inhomogeneous distribution of electronic sectors
across layers. Both ingredients are crucial to capture the effect
of interactions on the properties of the
N
-layered structure.
The single-particle Hamiltonian
H
sp
transforms into block-
diagonal form under a basis transformation
V
TNG
in layer
space [
16
]. For an even number
N
of layers, there are
N
/
2
=
N
/
2
blocks—or sectors. These blocks describe bands anal-
ogous to twisted bilayer graphene at twist angle
θ
with
interlayer hoppings rescaled by a coefficient
k
. We can
equivalently think of the sectors as corresponding to TBG
with unscaled hoppings, but an effective twist angle
θ
eff
k
=
θ/
k
.
(7)
In this picture, the sector Hamiltonian is multiplied by an
overall scale factor
k
.For
N
odd, in addition to the
N
/
2
TBG-like sectors, there is an additional sector, in which
the band derives from the underlying graphene Dirac cone
folded into the moiré Brillouin zone (BZ). We will denote
this sector as the monolayer-graphene (MLG)-like sector [see
Fig.
1(c)
]. We will choose the physical angle
θ
such that
there is one TBG-like sector—termed magic sector—at the
magic angle,
θ
eff
k
=
θ
magic
≈
1
.
1
◦
. In experiments to date, this
would be the
k
=
1 sector, but in Sec.
V
we also consider
the possibility
k
magic
=
2. We refer to all other sectors as
the nonmagic sectors, including the MLG-like sector for
N
odd [
16
]. Technically, the sector decomposition emerges by
solving an effective
N
/
2
-site open tight-binding chain on
the even layers, with
N
/
2
solutions (see Appendix
A1
for
a pedagogical derivation). The solution for the odd layers
proceeds analogously. The resulting weight distribution for
sector
k
is
W
(
k
)
l
=
2
N
+
1
sin
2
π
kl
N
+
1
,
(8)
as plotted in Fig.
1(b)
, with corresponding eigenvalues
k
=
2 cos
π
k
N
+
1
.
(9)
Combined with Eq. (
7
), this gives Eq. (
1
). With an increasing
number of layers, there is a continuum of twist angles [
16
],
with the largest density of twist angles close to the mini-
mal
θ
eff
k
(attained for
k
=
1). If the physical twist angle
θ
is such that the lowest effective angle sector is magic, there
195148-5
KRYŠTOF KOLÁ
ˇ
R
et al.
PHYSICAL REVIEW B
108
, 195148 (2023)
will thus be other sectors very close to the magic angle.
Moreover, by slightly decreasing the physical twist angle, one
can alternatively tune the larger effective angles to be magic
(see Sec.
V
). The weight distribution in Eq. (
8
) quantifies
the charge distributions across layers for the various sectors,
see Fig.
1(b)
. As discussed above, this is important for the
electrostatic properties of the problem.
B. Coulomb interactions
We assume a symmetric double-gated setup [see Fig.
1(a)
]
as typically employed in experiment. We work at gate charge
densities
en
/
2 per gate, so that
−
en
is the charge density in
TNG. We include Coulomb interactions through
Hint
=
1
2
d
r
d
r
V
(
r
−
r
):
ρ
(
r
)
ρ
(
r
):
,
(10)
where the density
ρ
(
r
) includes free charges in both the
graphene system and on the gates with the positive back-
ground subtracted (:
...
: denotes normal ordering). The
integration ranges over the full 3D space. Integrating out the
electronic degrees of freedom of the metallic gates, one arrives
at an effective screened interaction for the
N
layers for a
fixed electron density
n
(see Appendix
B1
for details). The
resulting interaction Hamiltonian takes the form
H
int
=
1
2
A
q
=
0
i
,
j
V
ij
(
q
):
ρ
i
,
q
ρ
j
,
−
q
:
+
N
−
1
i
=
1
A
ε
⊥
ε
0
d
l
(
E
⊥
i
,
i
+
1
)
2
2
.
(11)
Here,
ρ
i
,
q
is the electron density in layer
i
at in-plane mo-
mentum
q
,
E
⊥
i
,
i
+
1
denotes the uniform component of the
perpendicular electric field between layers
i
and
i
+
1,
A
is
the system area,
d
l
is the interlayer distance, and
V
ij
(
q
)is
the double-gate-screened layer-dependent Coulomb interac-
tion derived in Appendix
B2
. We allow the dielectric constant
of the
q
=
0term(
ε
⊥
) to differ from the dielectric constant
entering
V
ij
(
q
)(
ε
). Physically, the out-of-plane interaction
reflects the out-of-plane response of graphene, while the
q
=
0 component is governed by the dielectric properties of the
substrate. For graphene layers,
ε
⊥
has been estimated to be
around 2 [
46
,
47
], while
ε
is around 5 for hBN substrates
[
28
,
36
,
48
,
49
]. Larger values, accounting for remote band
screening, have also been investigated [
23
,
25
,
36
]. We treat the
dielectric constants as parameters. Without the second term,
Eq. (
11
) is the standard in-plane Coulomb interaction of a 2D
system with screening due to metallic gates. The second term
is not usually included, but is important for multilayer systems
as discussed in Sec.
II
.
C. Mean-field decoupling
We perform our numerical calculations by restricting the
full Hilbert space to a finite number of
N
active
bands with
N
flavor
spin
/
valley flavors and solving the mean-field Hartree-Fock
equations. Specific details of the numerical simulation are
provided in Appendix
E
. We search for the
N
active
×
N
active
density matrix [
P
f
(
k
)]
αβ
=
c
†
f
,
k
,α
c
f
,
k
,β
.Here
c
f
,
k
,β
anni-
hilates a flavor-
f
electron in the single-particle band
β
at
momentum
k
. The single-particle bands fall into sectors
k
∈
{
1
,...,
n
o
}
. We keep
N
remote
remote bands, which generate ad-
ditional Hartree and Fock interaction terms. In projecting onto
a finite set of bands, we are assuming frozen fully filled bands
below and empty bands above this set. To avoid overcounting
of interactions already present in monolayer graphene and
thus included in the BM model [
23
,
24
,
43
], we subtract a
mean-field Hamiltonian corresponding to a reference density
matrix
P
0
f
(
k
). This is implemented in the mean-field equa-
tions by replacing every
P
f
(
k
) with
δ
P
f
(
k
)
=
P
f
(
k
)
−
P
0
f
(
k
)
.
(12)
We choose the subtraction scheme [
23
,
24
,
27
,
43
] in which
P
0
f
(
k
) is the ground density matrix at charge neutrality with
the interlayer hoppings switched off. For bands far below
the charge-neutrality point, interlayer hoppings are ineffective
and this density matrix approximates that of fully filled TNG
bands. It therefore cancels with the remote-band-interaction
term to a good approximation [
24
,
28
], justifying retaining
only a finite number
N
remote
of remote bands.
For the in-plane term, the mean-field decoupling extends
the usual procedure detailed in previous studies [
16
,
23
–
26
,
41
,
42
] to include the layer dependence of
V
i
,
j
(
q
). The
resulting Hartree term reads
H
Hartree
=
1
A
i
,
j
G
ρ
i
,
G
V
i
,
j
(
G
)
ρ
j
,
−
G
,
(13)
where we introduce the projected layer density operator,
ρ
i
,
G
=
f
k
c
†
f
,
k
fi
G
(
k
)
c
f
,
k
, and denote the mean-field density
operator (with the appropriate subtraction) as
ρ
j
,
−
G
=
f
k
c
†
f
,
k
fj
−
G
(
k
)
c
f
,
k
=
f
k
tr
δ
P
T
f
(
k
)
fj
−
G
(
k
)
.
(14)
Here, the trace runs over the space of active bands. Similarly,
the Fock term reads
H
Fock
=−
1
A
f
i
,
j
q
,
k
V
ij
(
q
)
×
c
†
f
,
k
fi
q
(
k
)
δ
P
T
f
(
k
+
q
)
fj
−
q
(
k
+
q
)
c
f
,
k
,
(15)
where in contrast to the Hartree term, each flavor interacts
only with itself.
At the mean-field level, the charge distribution across the
layers enters the Hamiltonian through
H
layer
=−
e
l
ρ
l
,
0
V
l
,
(16)
where
V
l
is the potential and
ρ
l
,
0
the electron number (i.e., the
q
=
0 Fourier component of the electron density
ρ
l
,
q
) of layer
l
.ThetermH
layer
contributes nontrivially due to imperfect
screening of the gate electrodes by the layers and becomes
increasingly important as
N
grows. Appendix
B3
details a
formal derivation of H
layer
in Eq. (
16
) by decoupling the out-
of-plane term in Eq. (
11
). The difference of layer potentials
V
i
+
1
−
V
i
=−
d
l
E
⊥
i
,
i
+
1
(17)
195148-6