# Short vs. long range exchange interactions in twisted bilayer graphene

Alejandro Jimeno-Pozo<EMAIL_ADDRESS>Imdea Nanoscience, Faraday
9, 28049 Madrid, Spain  Zachary A. H. Goodwin  Departments of Physics and
Materials and the Thomas Young Centre for Theory and Simulation of Materials,
Imperial College London, South Kensington Campus, London SW7 2AZ, UK  John A.
Paulson School of Engineering and Applied Sciences, Harvard University,
Cambridge, MA 02138, USA  Pierre A. Pantaleón  Imdea Nanoscience, Faraday 9,
28049 Madrid, Spain  Valerio Vitale  Departments of Physics and Materials and
the Thomas Young Centre for Theory and Simulation of Materials, Imperial
College London, South Kensington Campus, London SW7 2AZ, UK  Dipartimento di
Fisica, Università di Trieste, Strada Costiera 11, 34151, Trieste, Italy
Lennart Klebl  I. Institute of Theoretical Physics, University of Hamburg,
Notkestrasse 9, 22607 Hamburg, Germany  Dante M. Kennes  Institute for Theory
of Statistical Physics, RWTH Aachen University, and JARA Fundamentals of
Future Information Technology, 52062 Aachen, Germany  Max Planck Institute for
the Structure and Dynamics of Matter, Center for Free Electron Laser Science,
22761 Hamburg, Germany  Arash Mostofi  Departments of Physics and Materials
and the Thomas Young Centre for Theory and Simulation of Materials, Imperial
College London, South Kensington Campus, London SW7 2AZ, UK  Johannes Lischner
Departments of Physics and Materials and the Thomas Young Centre for Theory
and Simulation of Materials, Imperial College London, South Kensington Campus,
London SW7 2AZ, UK  Francisco Guinea  Imdea Nanoscience, Faraday 9, 28049
Madrid, Spain  Donostia International Physics Center, Paseo Manuel de
Lardizabal 4, 20018 San Sebastian, Spain

(28th August 2024)

###### Abstract

We discuss the effect of long-range interactions within the self-consistent
Hartree-Fock (HF) approximation in comparison to short-range atomic Hubbard
interactions on the band structure of twisted bilayer graphene (TBG) at charge
neutrality for various twist angles. Starting from atomistic calculations, we
determine the quasi-particle band structure of TBG with Hubbard interactions
for various magnetic orderings: modulated anti-ferromagnetic (MAFM), nodal
anti-ferromagnetic (NAFM) and hexagonal anti-ferromagnetic (HAFM). Then, we
develop an approach to incorporate these magnetic orderings along with the HF
potential in the continuum approximation. Away from the magic angle, we
observe a drastic effect of the magnetic order on the band structure of TBG
compared to the influence of the HF potential. Near the magic angle, however,
the HF potential seems to play a major role on the band structure compared to
the magnetic order. These findings suggest that the spin-valley degenerate
broken symmetry state often found in HF calculations of charge neutral TBG
near the magic angle should favour magnetic order, since the atomistic Hubbard
interaction will break this symmetry in favour of spin polarization.

##  I Introduction

Magic-angle twisted bilayer graphene (TBG) has generated tremendous interest
in twinstronics [1, 2, 3] since the discovery of correlated insulating states
and superconductivity in the $\sim$1.1$\degree$ moiré supperlattice [4, 5].
The initial reports [6, 4, 5] indicated strong electron-electron correlations
in TBG which gives rise to unconventional superconductivity [5]. While this is
not unanimously agreed upon [7, 8, 9, 10, 11], TBG has also been found to host
strange metallic behaviour [12, 13], nematic order [14, 15, 16, 17], Dirac
revivals [18, 19], Pomeranchuk effect [20, 21] and Chern insulators [22, 23,
24, 25, 26], amongst other effects and phases [27, 28, 29, 30, 31, 32, 33].

To understand these phases, given that the magic-angle of TBG contains
$\sim$12,000 atoms [34, 35, 36, 37], it is typical for the low-energy
continuum model [38, 39], based on that of Bistritzer and MacDonald [40], to
be utilised. This theory couples states of the Dirac cones of each layer and
valley at different moiré crystal momenta, which causes the onset of flat
bands at $\sim$1.1$\degree$ [8]. The continuum model of TBG can naturally be
extended to include long-ranged Hartree-Fock interactions [41, 42, 43, 44, 45,
46, 47, 48, 49], since it is an expansion in the moiré crystal momenta (which
are very small values, corresponding to large length scales). This interacting
theory has provided some understanding into the phase diagram of TBG [7], in
terms of the superconducting phase, correlated insulating states, Dirac
revivals, pinning of van Hove singularities [41, 45, 48], for example.

Including short-ranged interactions in this continuum model has proven more
difficult, however. Short-ranged interactions can be included in the Wannier
orbital Hamiltonians of the flat bands [50, 51, 52, 53, 54, 55, 56, 57, 58],
which provides a reduced Hamiltonian matrix which can be solved with strongly
correlated methods [59, 60, 52, 61, 62], but this also has long-ranged
interactions. It is more natural, however, to include short-ranged
interactions in atomistic models, such as DFT [34, 35] or tight-binding (TB)
[36, 37], since the atomic-scale information is retained in such approaches.
For example, Klebl and Honerkamp [63] studied the magnetic phase diagram of
TBG based on on-site Hubbard interactions from RPA spin-susceptibility
calculations [64, 65]. Moreover, these atomistic approaches can also handle
long-ranged interactions, such as self-consistent Hartree interactions [66,
67, 68, 69, 70, 71].

A significant limiting factor of self-consistent atomistic approaches for
broken symmetry phases is their computational cost [72, 73, 69]. Some examples
exist in the literature [74, 71, 70, 69], but a full phase diagram - over
twist angle and doping level, amongst other experimental variables - has not
yet been achieved. For example, Stauber and González [70, 71] developed a
theory based on Green’s functions, where only some of the states were
retained. They were able to study long-ranged interactions, and also the
interplay between long and short ranged interactions, but only at
1.16$\degree$ and either at charge neutrality or $-2$ electrons per moiré unit
cell. Moreover, Vahedi et al. [69] investigated several twist angles
(1.08$\degree$, 1.30$\degree$ and 1.47$\degree$) but only focused on charge
neutrality. Usually, this has been overcome from re-scaling the TB parameters
[75, 69, 76, 77] or applying hydrostatic pressure [78, 79, 80], such that flat
bands can be created with unit cells only containing a few hundred to a few
thousand atoms.

Here we develop an approach which can include short-ranged interactions, such
as the on-site Hubbard interaction of the pz orbitals, in the continuum model.
Starting from the RPA spin susceptibility calculations of Klebl and Honerkamp
[63], we perform self-consistent atomistic Hubbard calculations to obtain the
mean-field magnetic order parameters for different ordering tendencies (at a
large twist angle and charge neutrality). We develop analytical forms for the
magnetic ordering in real space, which we are able to include in the continuum
model as an effective a scalar sublattice potential. In order to elucidate the
angle dependence of the interplay between the magnetic orderings and the
Hartree-Fock potential, we perform self-consistent Hartree-Fock calculations
at charge neutrality, to which we later add the effective magnetic potential
at different twist angles. Overall, it is found that the long range
contribution dominates at the magic angle, but away from the magic-angle,
these magnetic orders become more significant. We discuss the competition
between these long and short range exchange interactions in detail, and finish
with a discussion of future directions.

##  II Methods

###  II.1 Atomistic Calculations

We study commensurate moiré unit cells of TBG [36], starting from AA stacked
bilayers and rotating the top layer anticlockwise about an axis perpendicular
to the graphene sheets that passes through a carbon atom in each layer. The
moiré lattice vectors of the commensurate structures are
$\textbf{R}_{1}=n\textbf{a}_{1}+m\textbf{a}_{2}$ and
$\textbf{R}_{2}=-m\textbf{a}_{1}+(n+m)\textbf{a}_{2}$, where $n$ and $m$ are
integers which define the commensurate TBG structure, and
$\textbf{a}_{1}=(\sqrt{3}/2,-1/2)a_{0}$ and
$\textbf{a}_{2}=(\sqrt{3}/2,1/2)a_{0}$ are the lattice vectors of graphene,
where $a_{0}=2.46\leavevmode\nobreak\ \textrm{\AA}$ is the lattice constant of
graphene.

At small twist angles, TBG undergoes significant atomic relaxations [34, 81,
35, 82, 83, 84, 85]. We calculate these relaxations using a classical force
field implemented in the LAMMPS software package [86]. The interlayer
interactions are modelled using the AIREBO-Morse potential [87], while
intralyer interactions are described with the Kolmogorov-Crespi potential
[88].

To investigate the electronic structure of TBG, we use the Hamiltonian

$\mathcal{\hat{H}}=\sum_{i\sigma}\varepsilon_{i\sigma}\hat{c}^{\dagger}_{i\sigma}\hat{c}_{i\sigma}+\sum_{ij\sigma}[t(\textbf{r}_{i}-\textbf{r}_{j})\hat{c}^{\dagger}_{j\sigma}\hat{c}_{i\sigma}+\text{H.c.}],$
(1)

where $\varepsilon_{i\sigma}$ and $\hat{c}^{\dagger}_{i\sigma}$
($\hat{c}_{i\sigma}$) denote the on-site energy of atom $i$ with spin $\sigma$
and the electron creation (annihilation) operator associated with atom $i$ and
spin $\sigma$, respectively. The hopping parameters between atoms $i$ and $j$,
$t(\textbf{r}_{i}-\textbf{r}_{j})$, are calculated using the Slater-Koster
rules [89]

$t(\textbf{r})=V_{pp\sigma}(\textbf{r})\bigg{(}\dfrac{\textbf{r}\cdot\textbf{e}_{z}}{|\textbf{r}|}\bigg{)}^{2}+V_{pp\pi}(\textbf{r})\bigg{(}1-\dfrac{\textbf{r}\cdot\textbf{e}_{z}}{|\textbf{r}|}\bigg{)}^{2},$
(2)

where
$V_{pp\sigma}(\textbf{r})=V_{pp\sigma}^{0}\exp\\{q_{\sigma}(1-|\textbf{r}|/d_{AB})\\}\Theta(R_{c}-|\mathbf{r}|)$
and
$V_{pp\pi}(\textbf{r})=V_{pp\pi}^{0}\exp\\{q_{\pi}(1-|\textbf{r}|/a)\\}\Theta(R_{c}-|\mathbf{r}|)$.
We take the pre-factor for the $pp\sigma$-hopping and $pp\pi$-hopping to be
$V_{pp\sigma}^{0}=0.48$ eV and $V_{pp\pi}^{0}=-2.7$ eV, respectively. The
carbon-carbon bond length is given by $a=a_{0}/\sqrt{3}$ and the interlayer
separation is taken to be $d_{\textrm{AB}}=3.35\leavevmode\nobreak\
\textrm{\AA}$. We take the decay parameters of the Slater-Koster rules to be
$q_{\sigma}=d_{\textrm{AB}}/0.184a_{0}$ and $q_{\pi}=1/0.184\sqrt{3}$ [36,
37]. Hoppings between carbon atoms separated by more than
$R_{c}=10\leavevmode\nobreak\ \textrm{\AA}$ are neglected [90].

To include the effects of short-range Hubbard interactions, the on-site energy
is determined by the mean-field Hubbard interaction

$\varepsilon_{i\sigma}=U[n_{i\sigma^{\prime}}-\textrm{const.}],$ (3)

where $U$ is the Hubbard parameter of the carbon pz orbital, and
$n_{i\sigma^{\prime}}$ is the mean-field electron density on atom $i$ with the
spin $\sigma^{\prime}$ being the opposite to $\sigma$.

The electron density can be determined from the Bloch eigenstates
$\psi_{n\mathbf{k}\sigma}(\mathbf{r})$ (with subscripts $n$ and $\mathbf{k}$
denoting a band index and the crystal momentum, respectively) according to

$\begin{split}n_{\sigma}(\textbf{r})&=\sum_{n\textbf{k}}f_{n\textbf{k}\sigma}|\psi_{n\textbf{k}\sigma}(\textbf{r})|^{2}\\\
&=\sum_{j}n_{j\sigma}\chi_{j}(\textbf{r}),\end{split}$ (4)

where
$f_{n\textbf{k}\sigma}=\Theta(\varepsilon_{F}-\varepsilon_{n\mathbf{k}\sigma})$
is the occupancy of state $\psi_{n\mathbf{k}\sigma}$ with eigenvalue
$\varepsilon_{n\mathbf{k}\sigma}$ (where $\varepsilon_{F}$ is the Fermi
energy),
$\chi_{j}(\textbf{r})=\sum_{\textbf{R}}\phi_{z}^{2}(\textbf{r}-\textbf{t}_{j}-\textbf{R})$
(with $\mathbf{R}$ denoting the moiré lattice vectors) and $n_{j\sigma}$ is
the total number of electrons in the $j$-th orbital with spin $\sigma$. To
characterizes the magnetic ordering, we calculate the spin polarisation

$\zeta=\dfrac{n_{\uparrow}-n_{\downarrow}}{n_{\uparrow}+n_{\downarrow}}.$ (5)

To start the mean-field calculations, instead of guessing random magnetic
configurations, we perform RPA spin-susceptibility calculations, following the
methods outlined in Refs. 63, 64. The eigenvalues of these calculations
provide the critical interaction strength of an instability ($U_{c}$) and the
eignevector is the form of the magnetic order. By using these eigenvectors as
an initial on-site interaction, we can induce spin polarisation, which can
then be used to perform self-consistent calculations.

To obtain a self-consistent solution of the atomistic Hubbard equation, we use
a simple mixing scheme with a mixing parameter of 0.1 typically (0.1 of the
new electron spin density is mixed into the same spin density). When
determining the Fermi energy, the total electron number is again forced to be
$N+\nu$, but this does not restrict the spin densities to be the same. We mix
the up and down spin density by the same amount, instead of choosing to work
with the total electron density and magnetic order parameter, as we find it is
sometimes more stable.

Only the leading instability was able to be stabilised with this method [91].
Therefore, we also perform constrained calculations. The constrained
calculations require an analytical form for the magnetic order to be specified
(which is outlined in Section III.1), say $\zeta^{(j)}$, and we self-
consistently determine the magnitude of this magnetic order through projecting
onto the magnetic order parameter, $\zeta$. The on-site energy term is then
given by

$\varepsilon_{i\sigma}=\pm\dfrac{U}{2}\zeta_{i}^{(j)},$ (6)

where the sign depends on the spin. A self-consistent solution is obtained
through a linear mixing of the spin densities.

###  II.2 Continuum Model Calculations

The mini-Brillouin Zone (mBZ) of the continuum model is spanned by the two
reciprocal lattice vectors given by,
$\mathbf{G}_{1}=\frac{2\pi}{L_{m}}\left(\frac{1}{\sqrt{3}},1\right)$ and
$\mathbf{G}_{2}=\frac{4\pi}{L_{m}}\left(-\frac{1}{\sqrt{3}},0\right)$, where
$L_{m}=\frac{a_{0}}{2\sin{(\theta/2)}}$ is the moiré period and $a_{0}$ is the
lattice constant of graphene. These vectors form the basis to define any
reciprocal lattice vector, $\mathbf{G}_{i}=n\mathbf{G}_{1}+m\mathbf{G}_{2}$
with $n,m\in\mathbb{Z}$, $i\in\mathbb{N}$. The first star of reciprocal
lattice vectors, i.e. the six first $\mathbf{G}_{i}$, are defined by
$n,m\in[-1,1]$.

The non-interacting continuum model is 4-fold degenerate, since it accounts
for valley and spin quantum numbers, with the Hamiltonian at crystal momentum
$\mathbf{k}$ being written as

$\mathcal{\hat{H}}_{TBG,\xi}\left(\mathbf{k}\right)=\matrixquantity(\hat{H}_{1,\xi}\left(\mathbf{k}\right)&\hat{T}\\\
\hat{T}^{\dagger}&\hat{H}_{2,\xi}\left(\mathbf{k}\right)),$ (7)

where $\hat{H}_{l,\xi}$ is the continuum single layer graphene Hamiltonian of
valley $\xi$, given by

$\hat{H}_{l,\xi}\left(\mathbf{k}\right)=\xi\hbar
v_{F}\left(\mathbf{k}-\xi\mathbf{K}_{l}\right)\tau_{\theta,l}.$ (8)

with $v_{F}=(\sqrt{3}V_{pp\pi}^{0}a)/(2\hbar)$ denoting the Fermi velocity,
$\mathbf{K}_{l}$ is the position of the Dirac point of layer $l$,
$\tau_{\theta,l}=e^{\mathrm{i}\xi\tau_{z}\theta/2}\matrixquantity(\tau_{x},\xi\tau_{y})e^{-\mathrm{i}\xi\tau_{z}\theta/2}$,
with $\tau_{i}$ being the Pauli matrices acting on the sublattice degree of
freedom. The matrix $\hat{T}$ is a periodic function in the moiré unit cell
that hybridises layers. For small angles, the main contribution comes from the
first three reciprocal lattice vectors, $\mathbf{G}=(0,0)$,
$\mathbf{G}=\mathbf{G}_{1}$ and $\mathbf{G}=\mathbf{G}_{1}+\mathbf{G}_{2}$
[38]

$\begin{split}\hat{T}=\sum_{\mathbf{G}}\hat{T}(\mathbf{G})&=\matrixquantity(u_{1}\hfil&u_{2}\\\
u_{2}&u_{1})+\matrixquantity(u_{1}&u_{2}e^{-2\mathrm{i}\theta\pi/3}\\\
u_{2}e^{2\mathrm{i}\theta\pi/3}&u_{1})\\\
&+\matrixquantity(u_{1}\hfil&u_{2}e^{2\mathrm{i}\theta\pi/3}\\\
u_{2}e^{-2\mathrm{i}\theta\pi/3}&u_{1}),\end{split}$ (9)

where $u_{1}=0.0797$ eV and $u_{2}=0.0975$ eV [92] are, respectively, the
hopping amplitudes between AB/BA and AA stacking, which takes into account the
atomic relaxation in the continuum model.

To account for electron-electron interactions, we include the mean-field
Hartree-Fock terms to the Hamiltonian. The Hartree contribution to the
Hamiltonian is given by

$\mathcal{\hat{H}}_{H}=\sum_{i,\xi,\sigma}\int_{\Omega}\differential^{2}{\mathbf{r}}\psi_{\xi,\sigma}^{i,\dagger}\left(\mathbf{r}\right)\psi_{\xi,\sigma}^{i}\left(\mathbf{r}\right)V_{H}\left(\mathbf{r}\right),$
(10)

where $i\in[1,4]$ labels the combined sublattice and layer degree of freedom,
$\sigma$ accounts for the spin, $\Omega$ is the area of the mBZ, and the local
Hartree potential is given by

$V_{H}\left(\mathbf{r}\right)=\int_{\Omega}\differential^{2}{\mathbf{r^{\prime}}}v_{C}\left(\mathbf{r}-\mathbf{r^{\prime}}\right)\expectationvalue{\delta\rho\left(\mathbf{r^{\prime}}\right)}.$
(11)

Here
$\delta\rho\left(\mathbf{r}\right)\equiv\rho\left(\mathbf{r}\right)-\rho_{CN}\left(\mathbf{r}\right)$
denotes the fluctuation in charge density, with
$\rho\left(\mathbf{r}\right)=\sum_{\xi,\sigma}{\psi_{\xi,\sigma}^{\dagger}\left(\mathbf{r}\right)\psi_{\xi,\sigma}\left(\mathbf{r}\right)}$
corresponding to the charge density and $\rho_{CN}\left(\mathbf{r}\right)$ is
the average density corresponding to the non-interacting TBG at charge
neutrality point. We assume that the Coulomb interaction is screened by a
double-metallic gate [45]

$v_{C}\left(\mathbf{q}\right)=\frac{2\pi
e^{2}}{\epsilon}\frac{\tanh{\left(d|\mathbf{q}|\right)}}{|\mathbf{q}|},$ (12)

where $d=40$ nm is the distance to the metallic gates and $\varepsilon=10$ is
the dielectric constant [49, 58, 67].

The Fock contribution to the Hamiltonian is given by

$\mathcal{\hat{H}}_{F}=\sum_{i,j,\xi,\sigma}\int_{\Omega}\differential^{2}{\mathbf{r}}\differential^{2}{\mathbf{r^{\prime}}}\psi_{\xi,\sigma}^{i,\dagger}\left(\mathbf{r}\right)V_{F}^{ij}\left(\mathbf{r},\mathbf{r^{\prime}}\right)\psi_{\xi,\sigma}^{j}\left(\mathbf{r^{\prime}}\right),$
(13)

where $i,j$ run over the sublattice and layer indexes. The non-local Fock
potential is described by,

$V_{F}^{ij}=-\expectationvalue{\psi_{\xi,\sigma}^{j,\dagger}\left(\mathbf{r^{\prime}}\right)\psi_{\xi,\sigma}^{i}\left(\mathbf{r}\right)}v_{C}\left(\mathbf{r}-\mathbf{r^{\prime}}\right),$
(14)

As we want to express the matrix elements of $\mathcal{\hat{H}}_{F}$, defined
in Eq. (13), in the reciprocal space for which we must transform the non-local
Fock potential into the Fourier space. By this procedure we compute the Fock
matrix elements as

$\begin{split}&\matrixelement{\mathbf{k}+\mathbf{G},\xi,\sigma,i}{\mathcal{H}_{F}}{\mathbf{k^{\prime}}+\mathbf{G^{\prime}},\xi^{\prime},\sigma^{\prime},i^{\prime}}=\\\
&-\sum_{i,j,\xi,\sigma,n}\sum_{\mathbf{k^{\prime\prime}},\mathbf{G^{\prime\prime}}}\psi_{n,\xi,\sigma}^{i}\left(\mathbf{k}+\mathbf{G^{\prime}}+\mathbf{G^{\prime\prime}}\right)\psi_{n,\xi,\sigma}^{j,*}\left(\mathbf{k^{\prime}}+\mathbf{G}+\mathbf{G^{\prime\prime}}\right)\\\
&\qquad\times
v_{C}\left(\mathbf{k}-\mathbf{k^{\prime\prime}}+\mathbf{G^{\prime\prime}}\right),\end{split}$
(15)

where the index $n$ runs over the occupied bands at a given Fermi energy. For
the Hartree-Fock calculations we work with a continuum model of TBG expanded
up to the third star. We use a density of points between $2-6\times
10^{5}\leavevmode\nobreak\ \textrm{\AA}^{2}$ in the mBZ, depending on the
twisting angle. The convergence of the Hartree-Fock potential is normally
reached after $5-6$ self-consistency steps and it has been proved for
increasing density of points.

To include the $\alpha$-magnetic potential ($\alpha=\text{M, N, H}$) in the
continuum model we use a scalar sublattice and spin dependent potential
expressed through its harmonic decomposition in the first star reciprocal
lattice vectors

$\hat{H}^{\alpha}(\mathbf{G}_{i},\mathbf{G}_{j})=\sum_{i,j=0}^{6}U_{\delta}^{\alpha}(\mathbf{G}_{i}-\mathbf{G}_{j}),$
(16)

The full form of $U_{\delta}^{\alpha}(\mathbf{G}_{i}-\mathbf{G}_{j})$ is
discussed in the Section III.3. The final Hamiltonian that combines both the
effective magnetic potential derived from atomistic calculations and the
Hartree-Fock potential at half-filling is given by,

$\mathcal{\hat{H}}\left(\mathbf{k}\right)=\mathcal{\hat{H}}_{TBG}\left(\mathbf{k}\right)+\mathcal{\hat{H}}_{F}\left(\mathbf{k}\right)+\mathcal{\hat{H}}^{\alpha}.$
(17)

Note that this final Hamiltonian is not treated in a self-consistent way since
the effective magnetic Hamiltonian is just added to the self-consistent
Hartree-Fock Hamiltonian, which would be equivalent to a first-order
approximation of the magnetic orderings in perturbation theory.

##  III Results

###  III.1 Short Range Atomistic Hubbard Interactions

From the RPA ($\textbf{q}=0$) spin-susceptibility calculations [63, 64], a
number of leading antiferromagnetic instabilities are found: modulated
antiferrogmagnetic order (MAFM), nodal antiferromagnetic order (NAFM) and
hexagonal antiferromagnetic order (HAFM). To find which instability is the
ground state at a given twist angle and doping level, (mean-field) atomistic
Hubbard calculations must be performed. As these atomistic calculations are
extremely computationally expensive at the magic angle, we focus on
1.54$\degree$ at charge neutrality. At this twist angle and doping level, TBGs
leading instability was found to be MAFM (with a critical Hubbard interaction
of $U_{c}\approx 5.1$ eV), with NAFM and HAFM having slightly larger critical
interaction strengths (of $U_{c}\approx 5.4$ eV). For these latter
instabilities, we perform unconstrained and constrained atomistic Hubbard
calculations, as outlined in Section II.1.

The MAFM instability, as seen in Fig. 1(b) and (c), is characterised by a sub-
lattice oscillation in the magnetic order parameter
$\zeta=(n_{\uparrow}-n_{\downarrow})/(n_{\uparrow}+n_{\downarrow})$ that is
modulated throughout the moiré supercell [63]. In Fig. 1(b), which plots the
magnetic structure along the diagonal of the moiré supperlattice, sublattice A
is shown in black and sublattice B is shown in grey for the top graphene
layer. To show its real-space structure more clearly, Fig. 1(c) plots the
magnetic order parameter on sublattice B of the top layer over the moiré
superlattice. The constant contribution of the MAFM order is larger than the
moiré-scale variation, which means the sub-lattice polarisation is the same
throughout the moiré unit cell. The moiré scale variation of the magnetic
order peaks in the AA regions, as might be expected from the LDOS of TBG
peaking peaked in the AA regions [36]. The MAFM order can be approximated with
the following analytical form

$\zeta^{M}(\mathbf{r})\approx\zeta_{s}^{\prime}+\dfrac{\zeta_{s}}{6}\sum_{i=1}^{6}\cos(\textbf{G}_{i}\cdot\textbf{r}),$
(18)

where $\mathbf{G}_{i}$ are reciprocal lattice vectors, and the sublattice
oscillation $\zeta_{s}^{\prime}$ is the constant sublattice polarisation and
$\zeta_{s}$ describes how this sublattice polarisation changes on the moré
scale. For sublattice $A_{l}$ ($B_{l}$), where $l=1,2$ is the layer index, the
sign of the polarisation is $-$ ($+$), or vice versa. Note this equation
assumes the AA region is located at the origin of the moiré unit cell.

Figure 1: (a,d,g) Self-consistent quasi-particle band structures for the
studied magnetic orders of TBG at 1.54$\degree$ and charge neutrality along
the high symmetry path. For MAFM, $U=5.4$ eV and unconstrained calculations
were performed; whereas, for NAFM and HAFM, $U=5.94$ eV and constrained
calculations were performed. (b,e,h) Corresponding self-consistent magnetic
order parameter plotted along the diagonal of the moiré supperlattice, where
$\textbf{R}_{1}$ and $\textbf{R}_{2}$ are the moiré lattice vectors.
Sublattice A is shown in black and sublattice B is shown in grey for the top
graphene layer (bottom layer not shown). (c,f,i) Corresponding plots in real
space for a single layer and sublattice, where only sublattice B of the top
layer is shown. Note, the $U$’s were chosen to be slightly above the critical
values. For MAFM, if $U=5.94$ eV is used, the gap at the Dirac point is
extremely large [91].

In Fig. 1(b) and (c) the unconstrained, self-consistent $\zeta$ for
1.54$\degree$ at charge neutrality and $U=5.4$ eV is shown. The peaks of
$\zeta\approx 0.1$ reside in the AA regions. If one electron is delocalized
across all atoms in the moiré unit cell, which contains 5548 atoms at
1.54$\degree$, then $\zeta\approx 10^{-4}$. Therefore, the magnetic structures
from these atomistic calculations are involving many more electrons than just
the 4 flat band electrons, which is in agreement with other works [70, 71,
69].

In Fig. 1(a), we show the corresponding self-consistent quasi-particle band
structure. The different valleys, K and K’, have been identified by applying
the valley operator to the states (see Refs. [93, 77, 76] for details of this
calculation), and shown in solid black and dotted grey, respectively. Since
the MAFM order breaks C2 symmetry, it causes a gap to open at the Dirac cones
at the K/K’ points of the moiré Brillouin zone of TBG. This instability was
not able to be stabilised at $U=5.1$ eV, but for $U$’s larger than $5.4$ eV,
the constant contribution dominates and it becomes graphene-like
($\zeta_{s}^{{}^{\prime}}$ $\gg$ $\zeta_{s}$, such that there is only a
constant sub-lattice polarisation in each graphene sheet), with the gap at the
K/K’ points becoming very large ($100$’s of meV) [91].

Similarly, NAFM also has a moiré-scale peak in the magnetic order in the AA
regions, but it does not possess a constant contribution to $\zeta$, as seen
in the self-consistent values plotted in Fig. 1(e). The corresponding real-
space structure is shown in Fig. 1(f), where nodes in the magnetic order
around the AA region separate the regions of opposite signs of $\zeta$ spin
polarisation. This magnetic order is referred to as nodal anti-ferromagnetic
order because $\zeta$ goes through $0$ between the AA and AB/BA regions,
causing the sign of $\zeta$ to change on each sub-lattice between these types
of stacking [63]. Therefore, it can be described by

$\zeta^{N}(\mathbf{r})\approx\dfrac{\zeta_{s}}{6}\sum_{i=1}^{6}\cos(\textbf{G}_{i}\cdot\textbf{r}),$
(19)

This instability was not the leading instability [63], and we found the
unconstrained calculations could never stabilise this order, as it would
always eventually revert to MAFM order. Therefore, we performed constrained
mean-field Hubbard calculations to find the “excited” magnetic order, as
explained in Section II.1. The quasi-particle band structure is shown in Fig.
1(d), and again it is a Mott insulator. A large gap is obtained in the Dirac
cone because of the Hubbard interaction parameter being taken as $U=5.94$ eV.
This instability could be stabilised with $U=5.67$ eV and also larger $U$’s in
the constrained method.

Finally, the HAFM instability is similar to NAFM, but where the magnetic order
parameter peaks on the AB/BA regions instead of the AA regions. The self-
consistent HAFM magnetic structure is shown in Fig. 1(h) and (i). The sign of
$\zeta$ changes between the AB and BA regions of the moiré unit cell. As the
peaks of $\zeta$ occur on the AB/BA regions, which form a hexagonal lattice on
the moiré scale, this ordering is referred to as hexagonal anti-ferromagnetic
order. This magnetic order is shown in Fig. 1(h) and (i), and can be described
by the analytical form

$\zeta^{H}(\mathbf{r})\approx\dfrac{\zeta_{s}}{6}\sum_{i=1}^{6}\sin(\textbf{G}_{i}\cdot\textbf{r}).$
(20)

This instability never appears to be a leading instability, but its critical
Hubbard interaction was $\sim 5.4$ eV, which is only slightly higher than the
leading instability. Therefore, we again found that constrained calculations
were required to obtain mean-field values of its order parameter, as explained
in Section II.1. In Fig. 1(g) we show the mean-field quasi-particle band
structure for $U=5.94$ eV (this was the smallest $U$ which could stabilise the
order with constrained calculations), where the different valleys have been
coloured solid black and dotted grey. We find that this magnetic order does
not create a gap at the K/K’ point, despite it having a sublattice
oscillation. This is because of the moiré-scale sine nature of its variation,
which means it does not break C2 on the moiré scale. It does cause the valleys
to split at the K/K’ points, however, resulting in one valley being pushed
higher in energy and the other to lower energies. This is analogous to the
effect of a perpendicular electric field has on the electronic structure.

These calculations give some insight into the magnetic order from Hubbard
interactions in TBG. However, we performed these well away from the magic
angle and for Hubbard interaction parameters that are large ($U\approx 5.5$ eV
for these calculations, but the value is thought to be $\sim 4$ eV [75]).
Therefore, to understand the role of these instabilities close to the magic
angle, we aim to include these magnetic states in the continuum model.

###  III.2 Long Range Interactions in the Continuum Model

In Fig. 2 we show the Hartree-Fock quasi-particle band structure (solid blue
line, for all subplots) for a number of twist angles at charge neutrality, in
addition to the non-interacting band structure (solid red line, for all
subplots). At 1.54$\degree$, Figs. 2(a), (d) and (g), we find that the
Hartree-Fock potential slightly modifies the non-interacting band structure,
indicating that this twist angle is too large for the formation of an
insulating state. At a twist angle of 1.25$\degree$, Figs. 2(b), (e) and (h),
the non-interacting bands are flat enough for the onset of a small gap due to
the Hartree-Fock potential, significantly smaller than the bandwidth, at the
Dirac cones K and K’ points in the moiré Brillouin zone. Right at the magic
angle of 1.05$\degree$, Figs. 2(c), (f) and (i), the Hartree-Fock interactions
induce a large gap at the K and K’ points, on a similar scale to the
bandwidth. At half filling (charge neutrality) these insulating states are
characterised by a broken sublattice symmetry and a preserved spin and valley
symmetry with respect to the non-interacting picture. This implies that the
mean-field ground state associated to the Hartree-Fock band structure is
actually a linear combination of 4 states with different spin and valley
indexes. These calculations are in good agreement with a large body of
literature which investigates these long ranged interactions …

Figure 2: Comparison of low energy band structure, along the path $\text{K}\rightarrow\Gamma\rightarrow\text{M}\rightarrow\text{K'}$ in the mBZ, of the non-interacting case in red, HF in blue and HF+MAFM (a, b, c), HF+NAFM (d, e, f), HF+HAFM (g, h, i) in orange, for $\theta=1.54^{\circ}$ (a, d, g) $\theta=1.25^{\circ}$ (b, e, f) and $\theta=1.05^{\circ}$ (c, f, i). Dashed lines correspond to valley flip. Order |  | $|\delta_{1}|$ | $|\delta_{2}|$
---|---|---|---
MAFM ($U=5.4$ eV) |  | $12.47$ | $12.54$
NAFM ($U=5.94$ eV) |  | $0$ | $237.09$
HAFM ($U=5.94$ eV) |  | $0$ | $211.89$
Table 1: Parameters $|\delta_{1}|$ and $|\delta_{2}|$ extracted form the
atomistic Hubbard calculations for the different orderings considered in this
work. All values are expressed in meV.

###  III.3 Short Range Interactions in the Continuum Model

In Section III.1, we described the leading anti-ferromagnetic instabilities
obtained from RPA calculations of the atomistic model, performed mean-field
Hubbard calculations of these, and found analytical forms which approximate
the spin polarisation well, with the magnetic potential they create simply
being $U\zeta/2$. Taking inspiration from how the Hartree potential is
calculated in the atomistic model relative to the continuum model, a general
expression for the potential induced by these anti-ferromagnetic instabilities
would be

$U_{\delta}(\mathbf{r})=\delta_{1}\tau_{z}+\dfrac{1}{6}\sum_{i=1}^{6}\delta_{2}\tau_{z}e^{\mathrm{i}\mathbf{G}_{i}\cdot\mathbf{r}},$
(21)

where $\delta_{1}$ corresponds to the constant sublattice polarization
strength, $\delta_{2}(\mathbf{G}_{i})$ is the moiré modulated part of the
potential, i.e. the weight associated to the expansion of the potential in the
$i$-th BZ reciprocal lattice vector, and $\tau_{z}$ is the Pauli matrix.

For the MAFM instability both $\delta_{1}$ and $\delta_{2}$ are real valued.
In the case of NAFM order there is not a constant contribution to $\zeta^{N}$
so $\delta_{1}=0$, and $\delta_{2}$ is a real number. Similarly for the HAFM
instability $\delta_{1}=0$ but $\delta_{2}$ is a purely imaginary number so
the modulation is a sine function. Both MAFM and NAFM orderings are degenerate
in the valley index but the HAFM is not and the modulated contribution to the
potential must be complex conjugated when exchanging valleys. The value of the
parameters $\delta_{1}$ and $\delta_{2}$ are obtained numerically within the
self-consistent atomistic Hubbard calculation at a twist angle of
1.54$\degree$ and at charge neutrality, as described in Section III.1. Their
values depending on the magnetic orders are summarized in Tab. 1.

These parameters should be linear in the interaction strength and the spin-
polarised electron density, meaning an overall non-linear dependence on $U$,
but we shall not seek self-consistent solutions to this potential in the
continuum model here. Instead, these potentials are included in the continuum
model through perturbation theory on the converged self-consistent Hartree-
Fock calculations, as explained in Section II.2. Note that the sublattice is
polarised in the Hartree-Fock calculations, which means the overall sign of
the potential can be chosen in two different ways, but we always choose the
sign such that the sublattice polarisation matches, as this should increase
any gaps, and further lower the energy. In the Appendix, we show the band
structures in the continuum model (without the Hartree-Fock contribution) at
1.54$\degree$ with these perturbed potentials, and find good agreement with
the atomistic calculations.

For twist angles smaller than 1.54$\degree$, the Hubbard interaction in the
continuum model should scale as
$|\delta_{1/2}\left(\theta\right)|=|\delta_{1/2}\left(1.54^{{}^{\circ}}\right)|\left(\frac{\theta}{1.54^{\circ}}\right)^{2}$,
[41]. However, extending the parameters to the magic angle in this way does
not yield very physical results, as the gaps in the flat bands are much larger
than is ever found in experiments, as shown in Fig. A2(c), (f) and (i). This
is because to obtain a mean-field solution of the magnetic order at
1.54$\degree$ a $U=5.4-6$ eV was required. However, a more physical value is
$U=4$ eV, or smaller [75]. Moreover, the large interaction strength gives rise
to large polarised spin densities, as seen in Section III.1, where the $\zeta$
values indicated that not just the flat band electrons were being polarised.
Thus, we rescale these parameters to obtain more suitable estimates for these
for the changes to the electronic structure. In Fig. 2 we show the results
with a scaling factor of $3$ in the main text, with other scaling factors
being shown in the Appendix. As these perturbed calculations are not self-
consistent, we cannot say what is the ground state. However, we can interpret
the changes to the electronic structure, which provides information about
where these interactions might be significant.

First we shall discuss MAFM order, as shown in Figs. 2(a), (b) and (c). At a
twist angle of 1.54$\degree$, Fig. 2(a), we find a small gap at the K/K’
points, which is more significant than the Hartree-Fock distortions. For a
smaller twist angle of 1.25$\degree$, Fig. 2(b), a similar situation is found,
the magnetic order opens a gap at the K/K’ points while the Hartree-Fock
potential is responsible for minor adjustments in the band structure. At the
magic angle of 1.05$\degree$, Fig. 2(c), the situation totally changes, now
the magnetic potential only slightly modifies the gap at the K/K’ points,
while the Hartree-Fock contribution dominates the deformations to the
electronic structure.

Next, we move on to describing the NAFM ordering, as seen in Fig. 2(d), (e)
and (f). At the largest twist angle of 1.54$\degree$, Fig. 2(d), this magnetic
order creates a significant gap at the K/K’ points. This could be attributed
to the non-self-consistent nature of the calculations, and form using large
values of the parameters. This large band deformations persists at the smaller
twist angles of 1.25$\degree$, Fig. 2(e), and 1.05$\degree$, Fig. 2(f). Even
if smaller values of the parameters are utilised, this NAFM order induces
large band deformations, lowering the energy of the occupied valence band.
Therefore, it appears that this magnetic order can compete with the Hartree-
Fock contribution.

Finally, we describe the effect of HAFM. As can be seen in Figs. 2(g), (h) and
(i), the HAFM order does not create a gap at the K/K’ points. Instead it
causes the Dirac cones at K and K’ to shift up and down, respectively, for the
single spin and valley channel. At 1.54$\degree$, Fig. 2(g), this effect is
almost imperceptible but stronger than the one induced by the Hartree-Fock
interactions. For the smaller twist angle of 1.25$\degree$, Fig. 2(h) the
situation is similar but the energy gap between K and K’ points increases.
While at the magic-angle, Fig. 2(i), it slightly contributes to reshaping the
band structure, which in contrast is heavily affected by the Hartree-Fock
contribution, so we can safely say that the HAFM is a secondary effect in this
case.

##  IV Discussion

Overall, it appears that these magnetic potentials are more significant away
from the magic angle, but at angles close enough that there could still be
broken symmetry phases [64, 26]. For example, 1.25$\degree$ seems to be the
most significantly affected by these Hubbard potentials relative to the
Hartree-Fock contribution. At the magic angle, the Hartree-Fock contribution
dominates, and at large twist angles the effects are small relative to the
bandwidth, suggesting that these magnetic orders are not significant at these
twist angles. This twist angle dependence could suggest why the predictions of
Klebl et al. [64], in terms of the twist angle and doping dependence of
magnetic states, agreed well with subsequent experiments [26]. This could be
because these Hubbard interactions are important close to the onset of broken
symmetry phases, but close to the magic angle these Hubbard interactions are
dominated by long-ranged Hartree-Fock interactions [41].

The NAFM order appears to effect the electronic structure most significantly,
significantly lowering the eigenvalues of the occupied valence band, and
therefore, it is a possible candidate for magnetic order in TBG. In contrast,
the MAFM order effects the electronic structure more weakly. Finally, the HAFM
appears to only effect the electronic structure slightly. This magnetic order
should, however, couple to perpendicular electric fields [75, 91], which could
make this ordering tendency more important. These perturbative calculations
are interesting to be performed because they are a very natural explanation
for the correlated insulating states in TBG [4]. From the Hartree-Fock
calculations, a spin-valley degenerate insulating state is obtained. The
atomistic Hubbard interaction, however, should break this symmetry and cause
the onset of magnetic order.

We have focused on TBG here, but many more moiré materials comprised of
graphene exist [94, 95]. Perhaps the most promising ones are where there is a
$\pm\theta$ twist between each adjacent graphene layer [96, 97, 98, 99, 100,
101, 102, 93]. These moiré graphene multilayers have been shown to host highly
tunable superconducting phases [103, 104, 105, 106, 107], and as the number of
layers increases, the superconducting phase occurs over wider and wider doping
ranges [108, 109]. Fischer et al. [100] has shown similar types of magnetic
order occur in these systems, which means it will be possible to use the
developed method. Another example of moiré structures is graphene twisted on a
graphene multilayer, such as twisted mono-bilayer graphene [110, 111]. The
magnetic structure of these systems was shown to be more complex by Goodwin et
al. [112], which suggests the approach described here could be difficult to
utilise. Finally, another class of moiré graphene multilayers is twisted
bilayers composed of graphene multilayers, such as twisted double bilayer
graphene [113, 114, 115, 116, 117, 118]. Further investigation of this system
would be of interest.

##  V Conclusion

In summary, starting from atomistic methods, we studied several leading
magnetic instabilities of charge neutral TBG at a large twist. These
calculations permitted analytical forms for the magnetic ordering. These
Hubbard potentials were investigated perturbatively from self-consistent
Hartree-Fock calculations in the continuum model, allowing a comparison of
long and short range exchange interactions. From these calculations, our take-
home conclusions are:

  1. 1.

These atomistic Hubbard interactions break the spin-valley degeneracy of the
insulating state at charge neutrality of TBG obtained from self-consistent
Hubbard calculations. Therefore, these insulating states are likely to have
some Mottness.

  2. 2.

These magnetic orders are most significant for intermediate twist angles
between the magic-angle and angles where non-interacting physics is
sufficient. At the magic-angle, the Hartree-Fock contribution dominates, and
at large angles the bandwidth dominates.

  3. 3.

Out of the studied magnetic orders, nodal anti-ferromagnetic order appears to
be the most significant for changes to the electronic structure.

It is hoped that these results further motivate inclusion of atomistic effects
in the continuum model. Moreover, performing self-consistent magnetic
calculations should also be possible, and investigating such ordering
tendencies in other moiré graphene multilayers is possible now.

##  VI Acknowledgments

A.J.P., P.A.P. and F.G. acknowledge support from the Severo Ochoa programme
for centres of excellence in R&D (Grant No. SEV-2016-0686, Ministerio de
Ciencia e Innovación, Spain); from the European Commission, within the
Graphene Flagship, Core 3, grant number 881603 and from grants NMAT2D
(Comunidad de Madrid, Spain) and SprQuMat (Ministerio de Ciencia e Innovación,
Spain). ZG was supported through a studentship in the Centre for Doctoral
Training on Theory and Simulation of Materials at Imperial College London
funded by the EPSRC (EP/L015579/1). We acknowledge funding from EPSRC grant
EP/S025324/1 and the Thomas Young Centre under grant number TYC-101. We
acknowledge the Imperial College London Research Computing Service
(DOI:10.14469/hpc/2232) for the computational resources used in carrying out
this work. This project has received funding from the European Union’s Horizon
2020 research and innovation programme under the Marie Skłodowska-Curie grant
agreement No. 101067977. The Deutsche Forschungsgemeinschaft (DFG, German
Research Foundation) is acknowledged for support through RTG 1995, within the
Priority Program SPP 2244 “2DMP” and under Germany’s Excellence Strategy-
Cluster of Excellence Matter and Light for Quantum Computing (ML4Q) EXC2004/1
- 390534769. We acknowledge support from the Max Planck-New York City Center
for Non-Equilibrium Quantum Phenomena. Spin susceptibility calculations were
performed with computing resources granted by RWTH Aachen University under
projects rwth0496 and rwth0589.

Figure A1: Low energy band structure for $\theta=1.54^{\circ}$, along the path
$\text{K}\rightarrow\Gamma\rightarrow\text{M}\rightarrow\text{K'}$ in the mBZ.
(a) for MAFM, (b) for NAFM, (c) for HAFM. Dashed lines correspond to valley
flip. Figure A2: Comparison of low energy band structure of the non-
interacting case in red, HF in blue and HF+MAFM (a, b, c), HF+NAFM (d, e, f),
HF+HAFM (g, h, i) in orange, for $\theta=1.54^{\circ}$ (a, d, g)
$\theta=1.25^{\circ}$ (b, e, f) and $\theta=1.05^{\circ}$ (c, f, i). Dashed
lines correspond to valley flip. This scenario correspond to a scaling factor
of 1. Figure A3: Comparison of low energy band structure of the non-
interacting case in red, HF in blue and HF+MAFM (a, b, c), HF+NAFM (d, e, f),
HF+HAFM (g, h, i) in orange, for $\theta=1.54^{\circ}$ (a, d, g)
$\theta=1.25^{\circ}$ (b, e, f) and $\theta=1.05^{\circ}$ (c, f, i). Dashed
lines correspond to valley flip. This scenario correspond to a scaling factor
of 2.

## References

  * Carr  _et al._ [2020a] S. Carr, S. Fang, and E. Kaxiras, Electronic-structure methods for twisted moiré layers, Nat. Rev. Mater 5, 748 (2020a).
  * Carr  _et al._ [2017] S. Carr, D. Massatt, S. Fang, P. Cazeaux, M. Luskin, and E. Kaxiras, Twistronics: Manipulating the electronic properties of two-dimensional layered structures through their twist angle, Phys. Rev. B 95, 075420 (2017).
  * Kennes  _et al._ [2021] D. M. Kennes, M. Claassen, L. Xian, A. Georges, A. J. Millis, J. Hone, C. R. Dean, D. N. Basov, A. Pasupathy, and A. Rubio, Moiré heterostructures: a condensed matter quantum simulator, Nat. Phys. 17, 155 (2021).
  * et al. [2018] Y. C. et al., Correlated insulator behaviour at half-filling in magic-angle graphene superlattices, Nature 556, 80 (2018).
  * Cao  _et al._ [2018] Y. Cao, V. Fatemi, S. Fang, K. Watanabe, T. Taniguchi, E. Kaxiras, and P. Jarillo-Herrero, Unconventional superconductivity in magic-angle graphene superlattices, Nature 556, 43 (2018).
  * Cao  _et al._ [2016] Y. Cao, J. Y. Luo, V. Fatemi, S. Fang, J. D. Sanchez-Yamagishi, K. Watanabe, T. Taniguchi, E. Kaxiras, and P. Jarillo-Herrero, Superlattice-induced insulating states and valley-protected orbits in twisted bilayer graphene, Phys. Rev. Lett. 117, 116804 (2016).
  * Balents  _et al._ [2020] L. Balents, C. R. Dean, D. K. Efetov, and A. F. Young, Superconductivity and strong correlations in moiré flat bands, Nat. Phys. 16, 725 (2020).
  * Andrei and MacDonald [2020] E. Y. Andrei and A. H. MacDonald, Graphene bilayers with a twist, Nat. Mater. 19, 1265 (2020).
  * Saito  _et al._ [2020] Y. Saito, J. Ge, K. Watanabe, T. Taniguchi, and A. F. Young, Independent superconductors and correlated insulators in twisted bilayer graphene, Nat. Phys. 16, 926 (2020).
  * Stepanov  _et al._ [2020] P. Stepanov, I. Das, X. Lu, A. Fahimniya, K. Watanabe, T. Taniguchi, F. H. L. Koppens, J. Lischner, L. Levitov, and D. K. Efetov, Untying the insulating and superconducting orders in magic-angle graphene, Nature 583, 375 (2020).
  * Oh  _et al._ [2021] M. Oh, K. P. Nuckolls, D. Wong, R. L. Lee, X. Liu, K. Watanabe, T. Taniguchi, and A. Yazdani, Evidence for unconventional superconductivity in twisted bilayer graphene, Nature 600, 240 (2021).
  * Cao  _et al._ [2019] Y. Cao, D. Chowdhury, D. Rodan-Legrain, O.Rubies-Bigordà, K. Watanabe, T. Taniguchi, T. Senthil, and P. Jarillo-Herrero, Strange metal in magic-angle graphene with near Planckian dissipation, Phys. Rev. Lett. 124, 076801 (2019).
  * Polshyn  _et al._ [2019] H. Polshyn, M. Yankowitz, S. Chen, Y. Zhang, K. Watanabe, T. Taniguchi, C. R. Dean, and A. F. Young, Large linear-in-temperature resistivity in twisted bilayer graphene, Nat. Phys. 15, 1011 (2019).
  * Cao  _et al._ [2021a] Y. Cao, D. Rodan-Legrain, J. M. Park, F. N. Yuan, K. Watanabe, T. Taniguchi, R. M. Fernandes, L. Fu, and P. Jarillo-Herrero, Nematicity and competing orders in superconducting magic-angle graphene, Science 372, 264 (2021a).
  * et al. [2019a] A. K. et al., Maximized electron interactions at the magic angle in twisted bilayer graphene, Nature 572, 95 (2019a).
  * Jiang  _et al._ [2019] Y. Jiang, X. Lai, K. Watanabe, T. Taniguchi, K. Haule, J. Mao, and E. Y. Andrei, Charge order and broken rotational symmetry in magic-angle twisted bilayer graphene, Nature 573, 91 (2019).
  * et al. [2019b] Y. C. et al., Electronic correlations in twisted bilayer graphene near the magic angle, Nat. Phys. 15, 1174 (2019b).
  * et al. [2020a] U. Z. et al., Cascade of phase transitions and Dirac revivals in magic-angle graphene, Nature 582, 203 (2020a).
  * Wong  _et al._ [2020] D. Wong, K. P. Nuckolls, M. Oh, B. Lian, Y. Xie, S. Jeon, K. Watanabe, T. Taniguchi, B. A. Bernevig, and A. Yazdani, Cascade of electronic transitions in magic-angle twisted bilayer graphene, Nature 582, 198 (2020).
  * et al. [2021a] A. R. et al., Entropic evidence for a Pomeranchuk effect in magic-angle graphene, Nature 592, 214 (2021a).
  * Saito  _et al._ [2021] Y. Saito, F. Yang, J. Ge, X. Liu, T. Taniguchi, K. Watanabe, J. I. A. Li, E. Berg, and A. F. Young, Isospin Pomeranchuk effect in twisted bilayer graphene, Nature 592, 220 (2021).
  * Wu  _et al._ [2020] S. Wu, Z. Zhang, K. Watanabe, T. Taniguchi, and E. Y. Andrei, Chern insulators, van Hove singularities and topological flat bands in magic-angle twisted bilayer graphene, Nat. Mater. 21, 488 (2020).
  * Nuckolls  _et al._ [2020] K. P. Nuckolls, M. Oh, D. Wong, B. Lian, K. Watanabe, T. Taniguchi, B. A. Bernevig, and A. Yazdani, Strongly correlated Chern insulators in magic-angle twisted bilayer graphene, Nature 588, 610 (2020).
  * Das  _et al._ [2021] I. Das, X. Lu, J. Herzog-Arbeitman, Z.-D. Song, K. Watanabe, T. Taniguchi, B. A. Bernevig, and D. K. Efetov, Symmetry broken Chern insulators and “magic series” of Rashbalike Landau level crossings in magic angle bilayer graphene, Nat. Phys. 17, 710 (2021).
  * et al. [2021b] Y. X. et al., Fractional chern insulators in magic-angle twisted bilayer graphene, Nature 600, 439 (2021b).
  * Choi  _et al._ [2021] Y. Choi, H. Kim, C. Lewandowski, Y. Peng, A. Thomson, R. Polski, Y. Zhang, K. Watanabe, T. Taniguchi, J. Alicea, and S. Nadj-Perge, Interaction-driven band flattening and correlated phases in twisted bilayer graphene, Nature Physics 17, 1375 (2021).
  * Yankowitz  _et al._ [2019] M. Yankowitz, S. Chen, H. Polshyn, Y. Zhang, K. Watanabe, T. Taniguchi, D. Graf, A. F. Young, and C. R. Dean, Tuning superconductivity in twisted bilayer graphene, Science 363, 1059 (2019).
  * Xie  _et al._ [2019] Y. Xie, B. Lian, B. Jäck, X. Liu, C.-L. Chiu, K. Watanabe, T. Taniguchi, B. A. Bernevig, and A. Yazdani, Spectroscopic signatures of many-body correlations in magic-angle twisted bilayer graphene, Nature 572, 101 (2019).
  * et al. [2019c] X. L. et al., Superconductors, orbital magnets, and correlated states in magic angle bilayer graphene, Nature 574, 653 (2019c).
  * Sharpe  _et al._ [2019] A. L. Sharpe, E. J. Fox, A. W. Barnard, J. Finney, K. Watanabe, T. Taniguchi, M. A. Kastner, and D. Goldhaber-Gordon, Emergent ferromagnetism near three-quarters filling in twisted bilayer graphene, Science 365, 605 (2019).
  * Serlin  _et al._ [2020] M. Serlin, C. L. Tschirhart, H. Polshyn, Y. Zhang, J. Zhu, K. Watanabe, T. Taniguchi, L. Balents, and A. F. Young, Intrinsic quantized anomalous Hall effect in a moiré heterostructure, Science 367, 900 (2020).
  * et al. [2019d] H. Y. et al., Atomic and electronic reconstruction at van der Waals interface in twisted bilayer graphene, Nat. Mater. 18, 448 (2019d).
  * et al. [2020b] A. U. et al., Mapping the twist-angle disorder and Landau levels in magic-angle graphene, Nature 581, 47 (2020b).
  * Uchida  _et al._ [2014] K. Uchida, S. Furuya, J.-I. Iwata, and A. Oshiyama, Atomic corrugation and electron localization due to moiré patterns in twisted bilayer graphenes, Phys. Rev. B 90, 155451 (2014).
  * Lucignano  _et al._ [2019] P. Lucignano, D. Alfé, V. Cataudella, D. Ninno, and G. Cantele, Crucial role of atomic corrugation on the flat bands and energy gaps of twisted bilayer graphene at the magic angle, Phys. Rev. B 99, 195419 (2019).
  * de Laissardière  _et al._ [2010] G. T. de Laissardière, D. Mayou, and L. Magaud, Localization of Dirac electrons in rotated graphene bilayers, Nano Lett. 10, 804 (2010).
  * de Laissardière  _et al._ [2012] G. T. de Laissardière, D. Mayou, and L. Magaud, Numerical studies of confined states in rotated bilayers of graphene, Phys. Rev. B 86, 125413 (2012).
  * Lopes dos Santos  _et al._ [2007] J. M. B. Lopes dos Santos, N. M. R. Peres, and A. H. Castro Neto, Graphene bilayer with a twist: Electronic structure, Phys. Rev. Lett. 99, 256802 (2007).
  * Lopes dos Santos  _et al._ [2012] J. M. B. Lopes dos Santos, N. M. R. Peres, and A. H. Castro Neto, Continuum model of the twisted graphene bilayer, Phys. Rev. B 86, 155449 (2012).
  * Bistritzer and MacDonald [2011] R. Bistritzer and A. H. MacDonald, Moiré bands in twisted double-layer graphene, Proceedings of the National Academy of Sciences 108, 12233 (2011).
  * Guinea and Walet [2018] F. Guinea and N. R. Walet, Electrostatic effects, band distortions, and superconductivity in twisted graphene bilayers, PNAS 115, 13174 (2018).
  * Xie and MacDonald [2020] M. Xie and A. H. MacDonald, Nature of the correlated insulator states in twisted bilayer graphene, Phys. Rev. Lett. 124, 097601 (2020).
  * Bultinck  _et al._ [2020] N. Bultinck, E. Khalaf, S. Liu, S. Chatterjee, A. Vishwanath, and M. P. Zaletel, Ground state and hidden symmetry of magic-angle graphene at even integer filling, Phys. Rev. X 10, 031034 (2020).
  * Liu and Dai [2021] J. Liu and X. Dai, Theories for the correlated insulating states and quantum anomalous hall effect phenomena in twisted bilayer graphene, Phys. Rev. B 103, 035427 (2021).
  * Cea  _et al._ [2019a] T. Cea, N. R. Walet, and F. Guinea, Electronic band structure and pinning of fermi energy to van hove singularities in twisted bilayer graphene: A self-consistent approach, Phys. Rev. B 100, 205113 (2019a).
  * Zhang  _et al._ [2020] Y. Zhang, K. Jiang, Z. Wang, and F. Zhang, Correlated insulating phases of twisted bilayer graphene at commensurate filling fractions: A hartree-fock study, Phys. Rev. B 102, 035136 (2020).
  * Cea  _et al._ [2022] T. Cea, P. A. Pantaleón, N. R. Walet, and F. Guinea, Electrostatic interactions in twisted bilayer graphene, Nano Materials Science 4, 27 (2022), special issue on Graphene and 2D Alternative Materials.
  * Cea and Guinea [2021] T. Cea and F. Guinea, Coulomb interaction, phonons, and superconductivity in twisted bilayer graphene, Proceedings of the National Academy of Sciences 118, e2107874118 (2021).
  * Cea and Guinea [2020] T. Cea and F. Guinea, Band structure and insulating states driven by coulomb interaction in twisted bilayer graphene, Phys. Rev. B 102, 045107 (2020).
  * Zou  _et al._ [2018] L. Zou, H. C. Po, A. Vishwanath, and T. Senthil, Band structure of twisted bilayer graphene: Emergent symmetries, commensurate approximants, and wannier obstructions, Phys. Rev. B 98, 085435 (2018).
  * Po  _et al._ [2018] H. C. Po, L. Zou, A. Vishwanath, and T. Senthil, Origin of mott insulating behavior and superconductivity in twisted bilayer graphene, Phys. Rev. X 8, 031089 (2018).
  * Kang and Vafek [2019] J. Kang and O. Vafek, Strong coupling phases of partially filled twisted bilayer graphene narrow bands, Phys. Rev. Lett. 122, 246401 (2019).
  * Seo  _et al._ [2019] K. Seo, V. N. Kotov, and B. Uchoa, Ferromagnetic mott state in twisted graphene bilayers at the magic angle, Phys. Rev. Lett. 122, 246402 (2019).
  * Calderón and Bascones [2020] M. Calderón and E. Bascones, Interactions in the 8-orbital model for twisted bilayer graphene, Phys. Rev. B 102, 155149 (2020).
  * Carr  _et al._ [2019] S. Carr, S. Fang, H. C. Po, A. Vishwanath, and E. Kaxiras, Derivation of Wannier orbitals and minimal-basis tight-binding Hamiltonians for twisted bilayer graphene: a first-principles approach, Phys. Rev. Research 1, 033072 (2019).
  * Goodwin  _et al._ [2019a] Z. A. H. Goodwin, F. Corsetti, A. A. Mostofi, and J. Lischner, Twist-angle sensitivity of electron correlations in moiré graphene bilayers, Phys. Rev. B 100, 121106(R) (2019a).
  * Goodwin  _et al._ [2019b] Z. A. H. Goodwin, F. Corsetti, A. A. Mostofi, and J. Lischner, Attractive electron-electron interactions from internal screening in magic angle twisted bilayer graphene, Phys. Rev. B 100, 235424 (2019b).
  * Goodwin  _et al._ [2020a] Z. A. H. Goodwin, V. Vitale, F. Corsetti, D. Efetov, A. A. Mostofi, and J. Lischner, Critical role of device geometry for the phase diagram of twisted bilayer graphene, Phys. Rev. B 101, 165110 (2020a).
  * Kennes  _et al._ [2018] D. M. Kennes, J. Lischner, and C. Karrasch, Strong correlations and $d+\mathit{id}$ superconductivity in twisted bilayer graphene, Phys. Rev. B 98, 241407 (2018).
  * Roy and Juričić [2019] B. Roy and V. Juričić, Unconventional superconductivity in nearly flat bands in twisted bilayer graphene, Phys. Rev. B 99, 121407 (2019).
  * Hsu  _et al._ [2020] Y.-T. Hsu, F. Wu, and S. Das Sarma, Topological superconductivity, ferromagnetism, and valley-polarized phases in moiré systems: Renormalization group analysis for twisted double bilayer graphene, Phys. Rev. B 102, 085103 (2020).
  * Laksono  _et al._ [2018] E. Laksono, J. N. Leaw, A. Reaves, M. Singh, X. Wang, S. Adam, and X. Gu, Singlet superconductivity enhanced by charge order in nested twisted bilayer graphene fermi surfaces, Solid State Communications 282, 38 (2018).
  * Klebl and Honerkamp [2019] L. Klebl and C. Honerkamp, Inherited and flatband-induced ordering in twisted graphene bilayers, Phys. Rev. B 100, 155145 (2019).
  * Klebl  _et al._ [2021] L. Klebl, Z. A. H. Goodwin, A. A. Mostofi, D. M. Kennes, and J. Lischner, Importance of long-ranged electron-electron interactions for the magnetic phase diagram of twisted bilayer graphene, Phys. Rev. B 103, 195127 (2021).
  * Fischer  _et al._ [2021a] A. Fischer, L. Klebl, C. Honerkamp, and D. M. Kennes, Spin-fluctuation-induced pairing in twisted bilayer graphene, Phys. Rev. B 103, L041103 (2021a).
  * Rademaker  _et al._ [2019] L. Rademaker, D. A. Abanin, and P. Mellado, Charge smoothening and band flattening due to Hartree corrections in twisted bilayer graphene, Phys. Rev. B 100, 205114 (2019).
  * Goodwin  _et al._ [2020b] Z. A. H. Goodwin, V. Vitale, X. Liang, A. A. Mostofi, and J. Lischner, Hartree theory calculations of quasiparticle properties in twisted bilayer graphene, Electron. Struct. 2, 034001 (2020b).
  * Cheung  _et al._ [2021] C. S. Cheung, Z. A. H. Goodwin, V. Vitale, J. Lischner, and A. A. Mostofi, Hartree theory and crystal field of twisted double bilayer graphene near the magic-angle, Electron. Struct. 4, 025001 (2021).
  * Vahedi  _et al._ [2021] J. Vahedi, R. Peters, A. Missaoui, A. Honecker, and G. T. de Laissardiére, Magnetism of magic-angle twisted bilayer graphene, SciPost Phys. 11, 083 (2021).
  * González and Stauber [2021] J. González and T. Stauber, Magnetic phases from competing Hubbard and extended Coulomb interactions in twisted bilayer graphene, Phys. Rev. B 104, 115110 (2021).
  * González and Stauber [2020] J. González and T. Stauber, Time-reversal symmetry breaking versus chiral symmetry breaking in twisted bilayer graphene, Phys. Rev. B 102, 081118 (2020).
  * Sboychakov  _et al._ [2019] A. O. Sboychakov, A. V. Rozhkov, A. L. Rakhmanov, and F. Nori, Many-body effects in twisted bilayer graphene at low twist angles, Phys. Rev. B 100, 045111 (2019).
  * Sboychakov  _et al._ [2020] A. O. Sboychakov, A. V. Rozhkov, A. L. Rakhmanov, and F. Nori, Spin density wave and electron nematicity in magic-angle twisted bilayer graphene, Phys. Rev. B 102, 155142 (2020).
  * Chichinadze  _et al._ [2020] D. V. Chichinadze, L. Classen, and A. V. Chubukov, Nematic superconductivity in twisted bilayer graphene, Phys. Rev. B 101, 224513 (2020).
  * Gonzalez-Arraga  _et al._ [2017] L. A. Gonzalez-Arraga, J. L. Lado, F. Guinea, and P. San-Jose, Electrically controllable magnetism in twisted bilayer graphene, Phys. Rev. Lett. 119, 107201 (2017).
  * Ramires and Lado [2019] A. Ramires and J. L. Lado, Impurity-induced triple point fermions in twisted bilayer graphene, Phys. Rev. B 99, 245118 (2019).
  * Wolf  _et al._ [2019] T. M. R. Wolf, J. L. Lado, G. Blatter, and O. Zilberberg, Electrically tunable flat bands and magnetism in twisted bilayer graphene, Phys. Rev. Lett. 123, 096802 (2019).
  * Carr  _et al._ [2018] S. Carr, S. Fang, P. Jarillo-Herrero, and E. Kaxiras, Pressure dependence of the magic twist angle in graphene superlattices, Phys. Rev. B 98, 085144 (2018).
  * Lopez-Bezanilla [2019] A. Lopez-Bezanilla, Emergence of flat-band magnetism and half-metallicity in twisted bilayer graphene, Phys. Rev. Materials 3, 054003 (2019).
  * Chen  _et al._ [2022] X. Chen, S. Liu, J. N. Fry, and H.-P. Cheng, First-principles calculation of gate-tunable ferromagnetism in magic-angle twisted bilayer graphene under pressure, J. Phys.: Condens. Matter. 34, 385501 (2022).
  * Nam and Koshino [2017] N. N. T. Nam and M. Koshino, Lattice relaxation and energy band modulation in twisted bilayer graphene, Phys. Rev. B 96, 075311 (2017).
  * Jain  _et al._ [2017] S. K. Jain, V. Juric̆ić, and G. T. Barkema, Structure of twisted and buckled bilayer graphene, 2D Mater. 4, 015018 (2017).
  * Gargiulo and Yazyev [2018] F. Gargiulo and O. V. Yazyev, Structural and electronic transformation in low-angle twisted bilayer graphene, 2D Mater. 5, 015019 (2018).
  * Guinea and Walet [2019] F. Guinea and N. R. Walet, Continuum models for twisted bilayer graphene: Effect of lattice deformation and hopping parameters, Physical Review B 99, 10.1103/physrevb.99.205134 (2019).
  * Liang  _et al._ [2020] X. Liang, Z. A. H. Goodwin, V. Vitale, F. Corsetti, A. Mostofi, and J. Lischner, Effect of bilayer stacking on the atomic and electronic structure of twisted double bilayer graphene, Phys. Rev. B 102, 155146 (2020).
  * Plimpton [1995] S. Plimpton, Fast parallel algorithms for short-range molecular dynamics, J. Comp. Phys. 117, 1 (1995).
  * O’Connor  _et al._ [2015] T. C. O’Connor, J. Andzelm, and M. O. Robbins, Airebo-m: A reactive model for hydrocarbons at extreme pressures, J. Chem. Phys. 142, 024903 (2015).
  * Kolmogorov and Crespi [2005] A. N. Kolmogorov and V. H. Crespi, Registry-dependent interlayer potential for graphitic systems, Phys. Rev. B 71, 235415 (2005).
  * Slater and Koster [1954] J. C. Slater and G. F. Koster, Simplified LCAO method for the periodic potential problem, Phys. Rev. 94, 1498 (1954).
  * Angeli  _et al._ [2018] M. Angeli, D. Mandelli, A. Valli, A. Amaricci, M. Capone, E. Tosatti, and M. Fabrizio, Emergent D6 symmetry in fully-relaxed magic-angle twisted bilayer graphene, Phys. Rev. B 98, 235137 (2018).
  * Goodwin [2022] Z. A. H. Goodwin, Theory and simulation of moiré graphene multilayers (2022).
  * Koshino  _et al._ [2018] M. Koshino, N. F. Q. Yuan, T. Koretsune, M. Ochi, K. Kuroki, and L. Fu, Maximally localized wannier orbitals and the extended hubbard model for twisted bilayer graphene, Phys. Rev. X 8, 031087 (2018).
  * Lopez-Bezanilla and Lado [2020] A. Lopez-Bezanilla and J. L. Lado, Electrical band flattening, valley flux, and superconductivity in twisted trilayer graphene, Phys. Rev. Research 2, 033357 (2020).
  * Cea  _et al._ [2019b] T. Cea, N. R. Walet, and F. Guinea, Twists and the electronic structure of graphitic materials, Nano Lett. 19, 8683 (2019b).
  * Pantaleón  _et al._ [2021] P. A. Pantaleón, T. Cea, R. Brown, N. R. Walet, and F. Guinea, Narrow bands and electrostatic interactions in graphene stacks, 2D Mater. 8, 044006 (2021).
  * Ledwith  _et al._ [2021] P. J. Ledwith, E. Khalaf, Z. Zhu, S. Carr, E. Kaxiras, and A. Vishwanath, Tb or not tb? contrasting properties of twisted bilayer graphene and the alternating twist n-layer structures (n=3,4,5…), arXiv:2111.11060 (2021).
  * Khalaf  _et al._ [2019] E. Khalaf, A. J.Kruchkov, G. Tarnopolsky, and A. Vishwanath, Magic angle hierarchy in twisted graphene multilayers, Phys. Rev. B 100, 085109 (2019).
  * Zhu  _et al._ [2020] Z. Zhu, S. Carr, D. Massatt, M. Luskin, and E. Kaxiras, Twisted trilayer graphene: A precisely tunable platform for correlated electrons, Phys. Rev. Lett. 125, 116404 (2020).
  * Carr  _et al._ [2020b] S. Carr, C. Li, Z. Zhu, E. Kaxiras, S. Sachdev, and A. Kruchkov, Ultraheavy and Ultrarelativistic dirac quasiparticles in sandwiched graphenes, Nano Lett. 20, 3030 (2020b).
  * Fischer  _et al._ [2021b] A. Fischer, Z. H. Goodwin, A. A. Mostofi, J. Lischner, D. M. Kennes, and L. Klebl, Unconventional superconductivity in magic-angle twisted trilayer graphene, npj Quantum Mater. 7, 5 (2021b).
  * o Phong  _et al._ [2021] V. o Phong, P. A. Pantaleón, T. Cea, and F. Guinea, Band structure and superconductivity in twisted trilayer graphene, Phys. Rev. B 104, L121116 (2021).
  * Shin  _et al._ [2021] J. Shin, B. L. Chittari, and J. Jung, Stacking and gate tunable topological at bands, gaps and anisotropic strip patterns in twisted trilayer graphene, Phys. Rev. B 104, 045413 (2021).
  * Park  _et al._ [2021] J. M. Park, Y. Cao, K. Watanabe, T. Taniguchi, and P. Jarillo-Herrero, Tunable strongly coupled superconductivity in magic-angle twisted trilayer graphene, Nature 590, 249 (2021).
  * Cao  _et al._ [2021b] Y. Cao, J. M. Park, K. Watanabe, T. Taniguchi, and P. Jarillo-Herrero, Large Pauli limit violation and reentrant superconductivity in magic-angle twisted trilayer graphene, Nature 595, 526 (2021b).
  * Hao  _et al._ [2021] Z. Hao, A. M. Zimmerman, P. Ledwith, E. Khalaf, D. H. Najafabadi, K. Watanabe, T. Taniguchi, A. Vishwanath, and P. Kim, Electric field tunable superconductivity in alternating twist magic-angle trilayer graphene, Science 371, 1133 (2021).
  * Turkel  _et al._ [2022] S. Turkel, J. Swann, Z. Zhu, M. Christos, K. Watanabe, T. Taniguchi, S. Sachdev, M. S. Scheurer, E. Kaxiras, C. R. Dean, and A. N. Pasupathy, Orderly disorder in magic-angle twisted trilayer graphene, Science 376, 193–199 (2022).
  * Kim  _et al._ [2022] H. Kim, Y. Choi, C. Lewandowski, A. Thomson, Y. Zhang, R. Polski, K. Watanabe, T. Taniguchi, J. Alicea, and S. Nadj-Perge, Evidence for unconventional superconductivity in twisted trilayer graphene, Nature 606, 494–500 (2022).
  * Zhang  _et al._ [2022] Y. Zhang, R. Polski, C. Lewandowski, A. Thomson, Y. Peng, Y. Choi, H. Kim, K. Watanabe, T. Taniguchi, J. Alicea, F. von Oppen, G. Refael, and S. Nadj-Perge, Promotion of superconductivity in magic-anglegraphene multilayers, Science 377, 1538–1543 (2022).
  * Burg  _et al._ [2022] G. W. Burg, E. Khalaf, Y. Wang, K. Watanabe, T. Taniguchi, and E. Tutuc, Emergence of correlations in alternating twist quadrilayer graphene, Nat. Mater 21, 884–889 (2022).
  * et al. [2021c] Y. S. et al., Tunable van Hove singularities and correlated states in twisted trilayer graphene, Nat. Phys. 17, 619 (2021c).
  * et al. [2021d] S. C. et al., Electrically tunable correlated and topological states in twisted monolayer–bilayer graphene, Nat. Phys. 17, 374 (2021d).
  * Goodwin  _et al._ [2021] Z. A. H. Goodwin, L. Klebl, V. Vitale, X. Liang, V. Gogtay, X. van Gorp, D. M. Kennes, A. A. Mostofi, and J. Lischner, Flat bands, electron interactions, and magnetic order in magic-angle mono-trilayer graphene, Phys. Rev. Materials 5, 084008 (2021).
  * et al. [2020c] C. S. et al., Correlated states in twisted double bilayer graphene, Nat. Phys. 16, 520 (2020c).
  * Liu  _et al._ [2020] X. Liu, Z. Hao, E. Khalaf, J. Y. Lee, K. Watanabe, T. Taniguchi, A. Vishwanath, and P. Kim, Tunable spin-polarized correlated states in twisted double bilayer graphene, Nature 583, 221 (2020).
  * G. W  _et al._ [2019] B. G. W, J. Zhu, T. Taniguchi, K. Watanabe, A. H. MacDonald, and E. Tutuc, Correlated insulating states in twisted double bilayer graphene, Phys. Rev. Lett. 123, 197702 (2019).
  * et al. [2021e] C. Z. et al., Visualizing delocalized correlated electronic states in twisted double bilayer graphene, Nat. Commun. 12, 2516 (2021e).
  * Rubio-Verdú  _et al._ [2021] C. Rubio-Verdú, S. Turkel, Y. Song, L. Klebl, R. Samajdar, M. S. Scheurer, J. W. F. Venderbos, K. Watanabe, T. Taniguchi, H. Ochoa, L. Xian, D. M. Kennes, R. M. Fernandes, Á. Rubio, and A. N. Pasupathy, Moiré nematic phase in twisted double bilayer graphene, Nature Physics 18, 196 (2021).
  * Cao  _et al._ [2020] Y. Cao, D. Rodan-Legrain, O. Rubies-Bigorda, J. M. Park, K. Watanabe, T. Taniguchi, and P. Jarillo-Herrero, Tunable correlated states and spin-polarized phases in twisted bilayer–bilayer graphene, Nature 583, 215 (2020).


# FrenchMedMCQA: A French Multiple-Choice Question Answering Dataset for
Medical domain

Yanis Labrak1,4 Adrien Bazoge2 Richard Dufour2 Béatrice Daille2 Pierre-Antoine
Gourraud3 Emmanuel Morin2 Mickael Rouvier1 LIA - Avignon University1
<EMAIL_ADDRESS>LS2N - Nantes University2
<EMAIL_ADDRESS>CHU de Nantes - La clinique des données - Nantes
University3 Zenidoc 4

###### Abstract

This paper introduces FrenchMedMCQA, the first publicly available Multiple-
Choice Question Answering (MCQA) dataset in French for medical domain. It is
composed of 3,105 questions taken from real exams of the French medical
specialization diploma in pharmacy, mixing single and multiple answers. Each
instance of the dataset contains an identifier, a question, five possible
answers and their manual correction(s). We also propose first baseline models
to automatically process this MCQA task in order to report on the current
performances and to highlight the difficulty of the task. A detailed analysis
of the results showed that it is necessary to have representations adapted to
the medical domain or to the MCQA task: in our case, English specialized
models yielded better results than generic French ones, even though
FrenchMedMCQA is in French. Corpus, models and tools are available online.

FrenchMedMCQA: A French Multiple-Choice Question Answering Dataset for Medical
domain

Yanis Labrak1,4 Adrien Bazoge2 Richard Dufour2 Béatrice Daille2

Pierre-Antoine Gourraud3 Emmanuel Morin2 Mickael Rouvier1

LIA - Avignon University1<EMAIL_ADDRESS>LS2N - Nantes
University2<EMAIL_ADDRESS>

CHU de Nantes - La clinique des données - Nantes University3 Zenidoc 4

##  1 Introduction

Multiple-Choice Question Answering (MCQA) is a natural language processing
(NLP) task that consists in correctly answering a set of questions by
selecting one (or more) of the given $N$ candidates answers (also called
options) while minimizing the number of errors. MCQA is one of the most
difficult NLP tasks because it requires more advanced reading comprehension
skills and external sources of knowledge to reach decent performance.

In MCQA, we can distinguish two types of answers: (1) single and (2) multiple
ones. Most datasets focus on single answer questions, such as MCTest
Richardson et al. (2013), ARC-challenge Clark et al. (2018), OpenBookQA
Mihaylov et al. (2018), QASC Khot et al. (2019), Social-IQA Sap et al. (2019),
or RACE Lai et al. (2017). To our knowledge, few studies have been done to
construct medical MCQA dataset. We can cite the MedMCQA Pal et al. (2022) and
HEAD-QA Vilares and Gómez-Rodríguez (2019) corpora which contain single answer
questions in Spanish and English respectively. For the multiple answer
questions, MLEC-QA Li et al. (2021) provides 136k questions in Chinese
covering various biomedical sub-fields, such as clinic, public health and
traditional Chinese medicine.

The French community has recently greatly increased its efforts to collect and
distribute medical corpora. Even if no open language model is currently
available, we can cite the named entity recognition Névéol et al. (2014) and
information extraction Grabar et al. (2018) tasks. However, they remain
relatively classic, current approaches already reaching a high level of
performance.

In this article, we introduce FrenchMedMCQA, the first publicly available MCQA
corpus in French related to the medical field, and more particularly in the
pharmacological domain. This dataset contains questions taken from real exams
of the French diploma in pharmacy. Among the difficulties related to the task,
the questions asked may require a single answer for some and multiple ones for
others. We also propose to evaluate state-of-the-art MCQA approaches,
including an original evaluation of several word representations across
languages.

Main contributions of the paper concern (1) the distribution of an original
MCQA dataset in French related to the medical field, (2) a state-of-the-art
approach on this task and a first analysis of the results, and (3) an open
corpus, including tools and models, all available online.

##  2 The FrenchMedMCQA Dataset

In this section, we detail the FrenchMedMCQA dataset and discuss data
collection and distribution.

###  2.1 Dataset collection

The questions and their associated candidate answer(s) were collected from
real French pharmacy exams on the
remede111http://www.remede.org/internat/pharmacie/qcm-internat.html website.
This site was built around a community linked to the medical field (medicine,
pharmacy, odontology…), offering multiple information (news, job offers,
forums…) both for students and also professionals in these sectors of
activity. Questions and answers were manually created by medical experts and
used during examinations. The dataset is composed of 2,025 questions with
multiple answers and 1,080 with a single one, for a total of 3,105 questions.
Each instance of the dataset contains an identifier, a question, five options
(labeled from A to E) and correct answer(s). The average question length is
14.17 tokens and the average answer length is 6.44 tokens. The vocabulary size
is of 13k words, of which 3.8k are estimated medical domain-specific words
(i.e. related to the medical field). We find an average of 2.5 medical domain-
specific words in each question (17% of words in average of a question) and
2.0 in each answer (36% of words in average of an answer). On average, a
targeted medical domain-specific word is present in 2 questions and in 8
answers.

###  2.2 Dataset distribution

Table 1 presents the proposed FrenchMedMCQA dataset distribution for the
train, development (dev) and test sets detailed per number of answers (i.e.
number of correct responses per question). Globally, 70% of the questions are
kept for the train, 10% for validation and last 20% for testing.

# Answers | Training | Validation | Test | Total
---|---|---|---|---
1 | 595 | 164 | 321 | 1,080
2 | 528 | 45 | 97 | 670
3 | 718 | 71 | 141 | 930
4 | 296 | 30 | 56 | 382
5 | 34 | 2 | 7 | 43
Total | 2171 | 312 | 622 | 3,105

Table 1: FrenchMedMCQA dataset distribution.

##  3 Methods

The use alone of the question to automatically find the right answer(s) is not
sufficient in the context of a MCQA task. State-of-the-art approaches then
require external knowledge to improve system performances Izacard and Grave
(2020); Khashabi et al. (2020). In our case, we decide to build a two-step
retriever-reader architecture comparable to UnifiedQA Khashabi et al. (2020),
where the retriever job is to extract knowledge from an external corpus and
using it by the reader to predict the correct answers for each question.
Figure 1 presents the two-step general pipeline, first step being the
retriever module, that extracts external context from the question (see
Section 3.1), and second step being the reader, called here question-answering
module (see Section 3.2), that automatically selects answer(s) to the targeted
question.

Figure 1: Steps of the pipeline.

###  3.1 Retriever module

An external medical-related corpus fully composed of French has first been
collected from two online sources: Wikipedia life science and HAL, the latter
being an open archive run by the French National Centre for Scientific
Research (CNRS) where authors can deposit scholarly documents from all
academic fields. In our case, we focus on extracting papers and thesis from
various specialization, such as Human health and pathology, Cancerology,
Public health and epidemiology, Immunology, Pharmaceutical sciences,
Psychiatric disorders and Drugs. This results in 1 million of passages (i.e. a
portion of text that contains at least 100 characters) in HAL and 286k
passages in Wikipedia.

This corpus is then used as a context extension for a question. We therefore
used a retriever pipeline to automatically assign questions to the most likely
passage in the external source. Two retrieval approaches are compared in this
article:

  * •

BM25 Okapi Trotman et al. (2014) for the implementation of the base BM25
algorithm Robertson and Sparck Jones (1988).

  * •

SentenceTransformers framework Reimers and Gurevych (2019) is used to perform
semantic search using state-of-the-art language representations taken from
Huggingface’s Transformers library Wolf et al. (2019).

For both approaches, the goal is to embed each passage of the external corpus
into a vector space using one of the two representations. On its side, the
question is concatenated with the five options (i.e. answers associated to the
question) to form a new query embedded in the same vector space. Embeddings
from question and passages are finally compared to return the closest passages
of a query (here, the cosine similarity is the distance metric). For the
SentenceTransformers approach, we used a fast and non domain specific model
called MiniLMv2 Wang et al. (2020). Note that the 1-best passage is only used
in these experiments.

###  3.2 Question-answer module

A goal of our experiments was to compare baseline approaches regarding two
different paradigms. The first one is referred to a discriminative approach
and consists in assigning one of $N$ classes to the input based on their
projection in a multidimensional space. We also referred to it as a multi-
class task. At the opposite, the second method is a generative one which
consists of generating a sequence of tokens, also called free text, based on a
sequence of input tokens identical to the one used for the discriminative
approach. The difference with the discriminative approach lies in the fact
that we are not outputting a single class, like ABE for the question
6234176387997480960, but a sequence of tokens following the rules of the
natural language and referring to a combination of classes like A + B + E in
the case of our studied generative model (see Section 3.2.2).

####  3.2.1 Discriminative representations

Four discriminative representations are studied in this paper. We firstly
propose to use CamemBERT Martin et al. (2020), a generic French pre-trained
language model based on RoBERTa Liu et al. (2019). Since no language
representation adapted to the medical domain are publicly available for
French, we propose to evaluate the two pre-trained representations BioBERT Lee
et al. (2019) and PubMedBERT Gu et al. (2022), both trained on English medical
data and reaching SOTA results on biomedical NLP tasks, including QA Pal et
al. (2022). Finally, we consider a multilingual generic pre-trained model,
XLM-RoBERTa Conneau et al. (2020) based on RoBERTa, to evaluate the gap in
terms of performance with CamemBERT.

####  3.2.2 Generative representation

Recently, generative models have demonstrated their interest on several NLP
tasks, in particular for text generation and comprehension tasks. Among these
approaches, BART Lewis et al. (2019) is a denoising autoencoder built with a
sequence-to-sequence model. Due to its bidirectional encoder and left-to-right
decoder, it can be considered as generalizing BERT and GPT Radford et al.
(2019), respectively. BART training has two stages: (1) a noising function
used to corrupt the input text, and (2) a sequence-to-sequence model learned
to reconstruct the original input text. We then propose to evaluate this
representation in this paper.

##  4 Experimental protocol

Each studied discriminative and generative model is fine-tuned on the MCQA
task with FrenchMedMCQA training data using an input sequence composed of a
question, its associated options (i.e. possible answers) and its additional
context, all separated with a "[SEP]" token, e.g. [CLS] <question> [SEP] (A)
<answer.a> [SEP] (B) <answer.b> [SEP] (C) <answer.c> [SEP] (D) <answer.d>
[SEP] (E) <answer.e> [SEP] <context> [EOS].

For each question, the context is the text passage with highest confidence
rate and can either be obtained using the BM25 algorithm or semantic search as
described in Section 3.1.

Concerning the outputs of the systems, we have for the BART generative model a
plain text containing the letter of the answers from $A$ to $E$ separated with
plus signs in case of the questions with multiple answers, e.g. A + D + E. For
the other architectures (i.e. discriminative approaches), we simplify the
multi-label problem into a multi-class one by classifying the inputs into one
of the 31 existing combinations in the corpus. Here, a class may be a
combination of multiple labels, e.g. if the correct answers are the $A$ and
$B$ ones, then we consider the correct class being $AB$, which explains the
number of 31 classes.

| Without Context | Wiki w/ BM25 | HAL w/ BM25 | Wiki w/ MiniLMv2 | HAL w/ MiniLMv2
---|---|---|---|---|---
Architecture | Hamming | EMR | Hamming | EMR | Hamming | EMR | Hamming | EMR | Hamming | EMR
BioBERT V1.1 | 36.19 | 15.43 | 38.72 | 16.72 | 33.33 | 14.14 | 35.13 | 16.23 | 34.27 | 13.98
PubMedBERT | 33.98 | 14.14 | 34.00 | 13.98 | 35.66 | 15.59 | 33.87 | 14.79 | 35.44 | 14.79
CamemBERT-base | 36.24 | 16.55 | 34.19 | 14.46 | 34.78 | 15.43 | 34.66 | 14.79 | 34.61 | 14.95
XLM-RoBERTa-base | 37.92 | 17.20 | 31.26 | 11.89 | 35.84 | 16.07 | 32.47 | 14.63 | 33.00 | 14.95
BART-base | 31.93 | 15.91 | 34.98 | 18.64 | 33.80 | 17.68 | 29.65 | 12.86 | 34.65 | 18.32

Table 2: Performance (in %) on the test set using the Hamming score and EMR
metrics.

###  4.1 Evaluation metrics

The majority of tasks concentrate either on multi-class or binary
classification since they have a single class at a time. However,
occasionally, we will have a task where each observation has many labels. In
this case, we would have different metrics to evaluate the system itself
because multi-label prediction has an additional notion of being partially
correct. Here, we focused on two metrics called the Hamming score (commonly
also multi-label accuracy) and Exact Match Ratio (EMR).

####  4.1.1 Hamming score

The accuracy for each instance is defined as the proportion of the predicted
correct labels to the total number (predicted and actual) of labels for that
instance. Overall accuracy is the average across all instances. It is less
ambiguously referred to as the Hamming score rather than Multi-label Accuracy.

####  4.1.2 Exact Match Ratio (EMR)

The Exact Match Ratio (EMR) is the percentage of predictions matching exactly
the ground truth answers. To be computed, we sum the number of fully correct
questions divided by the total number of questions available in the set. A
question is considered fully correct when the predictions are exactly equal to
the ground truth answers for the question (e.g. all multiple answers should be
correct to count as a correct question).

##  5 Results

Table 2 compiled the performance (in terms of Hamming score and EMR) of all
the studied architectures and retrievers pipelines. For sake of comparison,
the column Without Context has been added, considering that no retriever is
used (i.e. no external passage is present in the QA system).

As we can see, the best performing model is different according to the used
metric. BioBERT V1.1 reaches best performance using the Hamming score and
BART-base in the case of the EMR. These first observations are quite
surprising since both models are trained on English data. While we could
expect higher performance with French models (CamemBERT for example), the fact
that these models are trained on specialized data for one (BioBERT) and on a
model designed for the targeted task (SOTA on question-answering for BART)
finally shows that language models trained on generic data are inefficient for
the MCQA task on medical domain.

In all considered architectures, context seems to have a small impact on
systems performance, with a limited increase or drop depending on the
configurations. Clearly, the RoBERTa performance is much higher without
context (i.e. without the use of the retriever part), while models based on
BERT generally (8 times on 12) outperform their own baseline performances with
external context. The fact that we consider the 1-best passage only may
explain this impact.

Concerning XLM-RoBERTa-base (cross lingual representation), we obtain in the
case of the context extracted using BM25 from Wikipedia, the worst Hamming
score and EMR out of all the discriminative approaches. This confirms our
first observation that a non-specialized model does not allow to achieve the
best performance on this task.

Using BM25 promotes better context than semantic search using MiniLMv2 on both
Wikipedia and HAL for most of the runs. Finally, the source depends of the
retriever and model used. A majority of the experiments demonstrate that HAL
outperforms Wikipedia on BM25 despite the fact that the best model was
obtained using Wikipedia.

The scripts to replicate the
experiments222https://github.com/qanastek/FrenchMedMCQA as well as the pre-
trained models333https://huggingface.co/qanastek/FrenchMedMCQA-
BioBERT-V1.1-Wikipedia-BM25/tree/main are available online.

##  6 Conclusion

We proposed in this paper FrenchMedMCQA, an original, open and publicly
available Multiple-Choice Question Answering (MCQA) dataset in the medical
field. This is the first French corpus in this domain, including single and
multiple answers to questions. Several state-the art systems have been
evaluated to show current performance on the dataset. The analysis of these
first results notably highlighted the fact that language models specialized to
the medical domain allow us to reach better performance than generic models,
even if these have been trained in a different language (here, English
biomedical models applied to French).

In future works, we will focus on improving the existing methods for the task
of MCQA, considering other strategies for the retriever module (multiple
passages, combining contexts…). Likewise, we will also consider the
construction of data representation models for French specialized for medical
domain.

##  7 Acknowledgments

This work was financially supported by Zenidoc, the DIETS project financed by
the Agence Nationale de la Recherche (ANR) under contract ANR-20-CE23-0005 and
the ANR AIBy4 (ANR-20-THIA-0011). This work was performed using HPC resources
from GENCI-IDRIS (Grant 2022-AD011013061R1 and 2022-AD011013715).

## References

  * Clark et al. (2018) Peter Clark, Isaac Cowhey, Oren Etzioni, Tushar Khot, Ashish Sabharwal, Carissa Schoenick, and Oyvind Tafjord. 2018.  Think you have solved question answering? try arc, the ai2 reasoning challenge.  _ArXiv_ , abs/1803.05457.
  * Conneau et al. (2020) Alexis Conneau, Kartikay Khandelwal, Naman Goyal, Vishrav Chaudhary, Guillaume Wenzek, Francisco Guzmán, Edouard Grave, Myle Ott, Luke Zettlemoyer, and Veselin Stoyanov. 2020.  Unsupervised cross-lingual representation learning at scale.  In _Proceedings of the 58th Annual Meeting of the Association for Computational Linguistics_ , pages 8440–8451, Online. Association for Computational Linguistics.
  * Grabar et al. (2018) Natalia Grabar, Vincent Claveau, and Clément Dalloux. 2018.  CAS: French corpus with clinical cases.  In _Proceedings of the Ninth International Workshop on Health Text Mining and Information Analysis_ , pages 122–128, Brussels, Belgium. Association for Computational Linguistics.
  * Gu et al. (2022) Yu Gu, Robert Tinn, Hao Cheng, Michael Lucas, Naoto Usuyama, Xiaodong Liu, Tristan Naumann, Jianfeng Gao, and Hoifung Poon. 2022.  Domain-specific language model pretraining for biomedical natural language processing.  _ACM Transactions on Computing for Healthcare_ , 3(1):1–23.
  * Izacard and Grave (2020) Gautier Izacard and Edouard Grave. 2020.  Leveraging passage retrieval with generative models for open domain question answering.  _arXiv preprint arXiv:2007.01282_.
  * Khashabi et al. (2020) Daniel Khashabi, Sewon Min, Tushar Khot, Ashish Sabharwal, Oyvind Tafjord, Peter Clark, and Hannaneh Hajishirzi. 2020.  Unifiedqa: Crossing format boundaries with a single qa system.
  * Khot et al. (2019) Tushar Khot, Peter Clark, Michal Guerquin, Peter Jansen, and Ashish Sabharwal. 2019\.  Qasc: A dataset for question answering via sentence composition.
  * Lai et al. (2017) Guokun Lai, Qizhe Xie, Hanxiao Liu, Yiming Yang, and Eduard Hovy. 2017.  RACE: Large-scale ReAding comprehension dataset from examinations.  In _Proceedings of the 2017 Conference on Empirical Methods in Natural Language Processing_ , pages 785–794, Copenhagen, Denmark. Association for Computational Linguistics.
  * Lee et al. (2019) Jinhyuk Lee, Wonjin Yoon, Sungdong Kim, Donghyeon Kim, Sunkyu Kim, Chan Ho So, and Jaewoo Kang. 2019.  BioBERT: a pre-trained biomedical language representation model for biomedical text mining.  _Bioinformatics_.
  * Lewis et al. (2019) Mike Lewis, Yinhan Liu, Naman Goyal, Marjan Ghazvininejad, Abdelrahman Mohamed, Omer Levy, Ves Stoyanov, and Luke Zettlemoyer. 2019.  Bart: Denoising sequence-to-sequence pre-training for natural language generation, translation, and comprehension.
  * Li et al. (2021) Jing Li, Shangping Zhong, and Kaizhi Chen. 2021.  MLEC-QA: A Chinese Multi-Choice Biomedical Question Answering Dataset.  In _Proceedings of the 2021 Conference on Empirical Methods in Natural Language Processing_ , pages 8862–8874, Online and Punta Cana, Dominican Republic. Association for Computational Linguistics.
  * Liu et al. (2019) Yinhan Liu, Myle Ott, Naman Goyal, Jingfei Du, Mandar Joshi, Danqi Chen, Omer Levy, Mike Lewis, Luke Zettlemoyer, and Veselin Stoyanov. 2019.  Roberta: A robustly optimized bert pretraining approach.
  * Martin et al. (2020) Louis Martin, Benjamin Muller, Pedro Javier Ortiz Suárez, Yoann Dupont, Laurent Romary, Éric de la Clergerie, Djamé Seddah, and Benoît Sagot. 2020.  CamemBERT: a tasty French language model.  In _Proceedings of the 58th Annual Meeting of the Association for Computational Linguistics_ , pages 7203–7219, Online. Association for Computational Linguistics.
  * Mihaylov et al. (2018) Todor Mihaylov, Peter Clark, Tushar Khot, and Ashish Sabharwal. 2018.  Can a suit of armor conduct electricity? a new dataset for open book question answering.  In _EMNLP_.
  * Névéol et al. (2014) Aurélie Névéol, Cyril Grouin, Jeremy Leixa, Sophie Rosset, and Pierre Zweigenbaum. 2014.  The QUAERO French medical corpus: A ressource for medical entity recognition and normalization.  In _Proc of BioTextMining Work_ , pages 24–30.
  * Pal et al. (2022) Ankit Pal, Logesh Kumar Umapathi, and Malaikannan Sankarasubbu. 2022.  Medmcqa: A large-scale multi-subject multi-choice dataset for medical domain question answering.  In _Proceedings of the Conference on Health, Inference, and Learning_ , volume 174 of _Proceedings of Machine Learning Research_ , pages 248–260. PMLR.
  * Radford et al. (2019) Alec Radford, Jeff Wu, Rewon Child, David Luan, Dario Amodei, and Ilya Sutskever. 2019.  Language models are unsupervised multitask learners.
  * Reimers and Gurevych (2019) Nils Reimers and Iryna Gurevych. 2019.  Sentence-bert: Sentence embeddings using siamese bert-networks.
  * Richardson et al. (2013) Matthew Richardson, Christopher J.C. Burges, and Erin Renshaw. 2013.  MCTest: A challenge dataset for the open-domain machine comprehension of text.  In _Proceedings of the 2013 Conference on Empirical Methods in Natural Language Processing_ , pages 193–203, Seattle, Washington, USA. Association for Computational Linguistics.
  * Robertson and Sparck Jones (1988) Stephen E. Robertson and Karen Sparck Jones. 1988.  _Relevance Weighting of Search Terms_ , page 143–160. Taylor Graham Publishing, GBR.
  * Sap et al. (2019) Maarten Sap, Hannah Rashkin, Derek Chen, Ronan LeBras, and Yejin Choi. 2019.  Socialiqa: Commonsense reasoning about social interactions.
  * Trotman et al. (2014) Andrew Trotman, Antti Puurula, and Blake Burgess. 2014.  Improvements to bm25 and language models examined.  In _Proceedings of the 2014 Australasian Document Computing Symposium_ , ADCS ’14, page 58–65, New York, NY, USA. Association for Computing Machinery.
  * Vilares and Gómez-Rodríguez (2019) David Vilares and Carlos Gómez-Rodríguez. 2019.  HEAD-QA: A healthcare dataset for complex reasoning.  In _Proceedings of the 57th Annual Meeting of the Association for Computational Linguistics_ , pages 960–966, Florence, Italy. Association for Computational Linguistics.
  * Wang et al. (2020) Wenhui Wang, Hangbo Bao, Shaohan Huang, Li Dong, and Furu Wei. 2020.  Minilmv2: Multi-head self-attention relation distillation for compressing pretrained transformers.
  * Wolf et al. (2019) Thomas Wolf, Lysandre Debut, Victor Sanh, Julien Chaumond, Clement Delangue, Anthony Moi, Pierric Cistac, Tim Rault, Rémi Louf, Morgan Funtowicz, Joe Davison, Sam Shleifer, Patrick von Platen, Clara Ma, Yacine Jernite, Julien Plu, Canwen Xu, Teven Le Scao, Sylvain Gugger, Mariama Drame, Quentin Lhoest, and Alexander M. Rush. 2019.  Huggingface’s transformers: State-of-the-art natural language processing.