Ferromagnetic Mott State in Twisted Graphene Bilayers at the Magic Angle

TL;DR

This paper investigates the Mott insulating state in twisted bilayer graphene at the magic angle, revealing a ferromagnetic insulator with distinct experimental signatures based on an effective spin model.

Contribution

It derives an effective spin model for the ferromagnetic Mott state in twisted bilayer graphene, highlighting its exotic magnetic properties and experimental signatures.

Findings

Identifies a ferromagnetic Mott insulator in twisted bilayer graphene.

Derives an effective spin Hamiltonian from the tight-binding model.

Predicts experimental signatures of the ferromagnetic Mott state.

Abstract

We address the effective tight-binding Hamiltonian that describes the insulating Mott state of twisted graphene bilayers at a magic angle. In that configuration, twisted bilayers form a honeycomb superlattice of localized states, characterized by the appearance of flat bands with four-fold degeneracy. After calculating the maximally localized superlattice Wannier wavefunctions, we derive the effective spin model that describes the Mott state. We suggest that the system is an exotic ferromagnetic Mott insulator, with well defined experimental signatures.

Tables

Table 1. Table 1: Electronic hopping amplitude | t α ​ α | subscript 𝑡 𝛼 𝛼 |t_{\alpha\alpha}| , direct Coulomb interaction V α ​ β subscript 𝑉 𝛼 𝛽 V_{\alpha\beta} and exchange interaction J α ​ β subscript 𝐽 𝛼 𝛽 J_{\alpha\beta} for various nearest neighbor sites: on-site ( 0 0 ), and n 𝑛 n -th nearest neighbors ( n 𝑛 n ), with n = 1 − 6 𝑛 1 6 n=1-6 . Energies in meV calculated for κ = 5 𝜅 5 \kappa=5 . n = 1 ,  3 𝑛 1 3 n=1,\>3 and 4 4 4 correspond to sites in opposite sublattices. J α ​ α ≈ ± J α , − α subscript 𝐽 𝛼 𝛼 plus-or-minus subscript 𝐽 𝛼 𝛼 J_{\alpha\alpha}\approx\pm J_{\alpha,-\alpha} , with + + ( − - ) for sites in the same (opposite) sublattice.

(meV)		0	1	2	3	4	5	6
|$t_{\alpha \alpha }|$		0	0.384	0.005	0.447	0.162	0.084	0.007
$V_{\alpha \beta }$		21.2	16.9	16.7	15.6	12.6	11.58	9.68
$J_{\alpha \alpha }$		0	5.09	1.11	0.52	0.25	0.16	0.09
$J_{\alpha ,-\alpha }$		0	$-$4.93	1.02	$-$0.51	$-$0.18	0.12	0.08

Equations

$$\mathcal{H}=-\sum _{ij}{J}_{ij}\left(\frac{1}{2}+2{\tau }_i^z{\tau }_j^z-2{\eta }_{ij}{\mathbf{\tau }}_i^{\perp }\cdot {\mathbf{\tau }}_j^{\perp }\right)\left(\frac{1}{2}+2{\mathbf{S}}_i\cdot {\mathbf{S}}_j\right),$$

$$\mathcal{H}=\left(\begin{array}[]{cc}\mathcal{H}_{g}^{(1)}&\mathcal{H}_{\perp}\\ \mathcal{H}_{\perp}^{\dagger}&\mathcal{H}_{g}^{(2)}\end{array}\right),$$

$${\Phi }_{\mathbf{k},\sigma }\equiv (|{\phi }_{\mathbf{k},A,\sigma }^{(1)}⟩,|{\phi }_{\mathbf{k},B,\sigma }^{(1)}⟩,|{\phi }_{\mathbf{k},A,\sigma }^{(2)}⟩,|{\phi }_{\mathbf{k},B,\sigma }^{(2)}⟩)$$

$${\mathcal{H}}_{\mathbf{k}}(\mathbf{r},{\mathbf{r}}^{\prime })=\sum _i\mathcal{H}(\mathbf{r},{\mathbf{r}}^{\prime }+{\mathbf{T}}_i){\text{e}}^{i\mathbf{k}\cdot {\mathbf{T}}_i}$$

$$[{\mathcal{H}}_{\mathbf{k}}{]}_{ab}=t_{\mathbf{k}}^{ab}(\mathbf{r},{\mathbf{r}}^{\prime })=\sum _j{\text{e}}^{i\mathbf{k}\cdot {\mathbf{T}}_j}{t}^{ab}(\mathbf{r},{\mathbf{r}}^{\prime }+{\mathbf{T}}_j),$$

$$|\mathbf{R}\nu ⟩=\frac{1}{N_s}\sum _{\mathbf{k}}{\text{e}}^{-i\mathbf{k}\cdot \mathbf{R}}{U}_{n\nu }(\mathbf{k})|{\hat{\psi }}_{n,\mathbf{k}}⟩,$$

$$\mathcal{H}=\sum _{i,j}{t}_{\alpha \beta }({\mathbf{R}}_{ij})d_{\alpha ,\sigma }^{†}({\mathbf{R}}_i)d_{\beta ,\sigma }({\mathbf{R}}_j),$$

$$t_{\alpha \beta }({\mathbf{R}}_i)=⟨\mathbf{R}\alpha |\mathcal{H}|\mathbf{R}+{\mathbf{R}}_i,\beta ⟩.$$

$${\mathcal{H}}_d=\sum _{i,j}{V}_{\alpha \beta }({\mathbf{R}}_{ij})n_{\alpha }({\mathbf{R}}_i)n_{\beta }({\mathbf{R}}_j),$$

$${\mathcal{H}}_e=\sum _{i\ne j}{J}_{\alpha {\alpha }^{\prime },\beta {\beta }^{\prime }}({\mathbf{R}}_{ij})d_{\alpha ,\sigma }^{†}({\mathbf{R}}_i)d_{\beta ,{\sigma }^{\prime }}^{†}({\mathbf{R}}_j)d_{{\beta }^{\prime },{\sigma }^{\prime }}({\mathbf{R}}_i)d_{{\alpha }^{\prime },\sigma }({\mathbf{R}}_j),$$

$$J_{\alpha {\alpha }^{\prime },\beta {\beta }^{\prime }}({\mathbf{R}}_{ij})$$

$$\times {\hat{W}}_{\beta }^{†}({\mathbf{r}}_j^{\prime })\cdot {\hat{W}}_{{\beta }^{\prime }}({\mathbf{r}}_i^{\prime })$$

Full text

*Introduction.$-$ *Mott insulators describe materials that exhibit insulating behavior as a result of strong local interactions imada . In those systems, strong on site repulsion penalizes the kinetic energy for electrons to hop between sites, rendering the electronic orbitals localized. The strong degree of localization of the electronic wavefunction favors antiferromagnetic alignment of the spins due to Pauli principle Auerbach . Recent experiments Cao1 ; Cao2 indicate that twisted graphene bilayers have a Mott state with an activation gap of $\Delta\approx 0.3$ meV that undergoes a metal-insulator transition in the vicinity of a superconducting phase Cao2 ; Yankowitz . This system is purely made of carbon atoms, with additional degrees of freedom inherited from graphene Castro Neto . That has motivated the question of whether the observed state could be described by a novel Mott insulator Po2 or other exotic correlated states Padhi ; Xu ; Irkhin ; Thomson ; dodaro . Unveiling the nature of the insulating state may be key to explain some of the the remarkable properties in the metallic phase.

By twisting two graphene sheets at a small angle of the order of $\theta\sim 1.1^{\circ}$, what was dubbed a “magic” angle, interference due to hopping between the layers leads to a Moire pattern and to a significant reconstruction of the mini bands in the Moire Brillouin zone, which become flat santos ; Bistritzer . Those flat bands have four-fold degeneracy, which is reminiscent of the valley and spin quantum numbers of the graphene sheets. In general, the confinement of interacting Dirac fermions in flat bands is expected to create an emergent SU(4) symmetry, as previously predicted in graphene heterostructures Uchoa1 ; Dou ; Xu3 and in graphene Landau levels Goerbig ; Young0 ; Young ; Nomura ; Alicea ; Abanin ; Sodemann . Here, the Moire pattern forms a superlattice of quasi-localized states with the size of the unit cell set by the twist angle, as shown in Fig. 1.

In this Letter, we show that the low energy Hamiltonian of the flat bands at quarter filling maps into the ferromagnetic spin exchange Hamiltonian on a honeycomb superlattice,

[TABLE]

where $\mathbf{S}_{i}$ is the localized spin on a superlattice site $i$, $\bm{\tau}_{i}=(\tau^{x},\tau^{y},\tau^{z})\equiv(\bm{\tau}^{\perp},\tau^{z})$ is an orbital pseudospin operator that is reminiscent of the valley quantum numbers, and $J_{ij}>0$ is the exchange coupling. The parameter $\eta_{ij}=-1$ when $i,j$ belong to the same sublattice, in which case the exchange interaction has SU(4) symmetry, and $\eta_{ij}=1$ otherwise, including nearest neighbor (NN) sites. This Hamiltonian acts in the Hilbert space which is spanned by four degenerate states per site, $|\alpha,\sigma\rangle$, with $\alpha=\pm$ and $\sigma=\uparrow,\downarrow$ for the two orbital pseudospins and spin quantum numbers respectively.

The existence of direct exchange ferromagnetism in an insulating state is uncommon Uchoa1 and reflects the very unusual shape of the Wannier orbitals in this system. Ferromagnetism has been recently observed in insulating van der Waals heterostructures of magnetic chromium trihalide materials, CrX3 (X$=$I, Br, Cl) Huang ; Yelon ; Samuelson , which have crystalline field anisotropies that produce an ordered Ising state. To the best of our knowledge, we are not aware of any examples of ferromagnetic Mott states which do not involve orbital ordering via a superexchange mechanism Erickson ; Kugel2 .

After performing calculations of the maximally localized Wannier orbitals of the Moire superstructure, we establish the parameters of a minimal interacting tight-binding model that captures the Mott physics near the magic angle. We show that even though the orbitals are well localized in the Mott regime at quarter filling, surprisingly the direct exchange interaction between different sites is dominant and favors ferromagnetic spin order at zero temperature. While charging effects Pizarro ; Guinea , which were not taken into account, may change our conclusions, the scenario of zero temperature ferromagnetism in twisted graphene bilayers seems in line with the reduced degeneracy of the Landau levels measured with Shubnikov de Haas experiments near quarter filling Cao1 . We discuss the experimental signatures of this state.

*Bloch Hamiltonian.$-$ *The free Hamiltonian for twisted graphene bilayers can be constructed at the lattice level using a parametrization for the hopping amplitudes between sites in the two different sheets,

[TABLE]

where *$\mathcal{H}_{g}$ *is the graphene Hamiltonian and *$\mathcal{H}_{\perp}$ *is the interlayer hopping between the two sheets in real space. The Moire pattern can be used to construct Bloch states that are periodic in the superlattice vectors $\mathbf{T}_{i}$. For commensurate structures, the Moire lattice vectors are parametrized by two integers $m$ and $r$, and correspond to the twist angle $\cos\theta=1-r^{2}/2(3m^{2}+3mr+r^{2})$, or equivalently $\theta\approx r/\sqrt{3}m$ for small angles.

In a basis for Bloch states

[TABLE]

defined in the two sublattices $A$ and $B$ of each of the two layers $(1,2)$, the Bloch Hamiltonian of the twisted system

[TABLE]

satisfies $\mathcal{H}_{\mathbf{k}}(\mathbf{r},\mathbf{r}^{\prime}+\mathbf{T}_{i})=\mathcal{H}_{\mathbf{k}}(\mathbf{r},\mathbf{r}^{\prime})e^{-i\mathbf{k}\cdot\mathbf{T}_{i}}$. In that basis,

[TABLE]

are the matrix elements of (2), with $a,b$ indexes running over the four components of basis (3). The hopping amplitudes $t^{ab}(\mathbf{r},\mathbf{r}^{\prime})=\cos^{2}\theta_{z}V_{\sigma}(\mathbf{r}-\mathbf{r}^{\prime})+\sin^{2}\theta V_{\pi}(\mathbf{r}-\mathbf{r}^{\prime})$, where $\cos\theta_{z}=d/\sqrt{d^{2}+(\mathbf{r}-\mathbf{r}^{\prime})^{2}}$ with $d$ the distance between the planes. $V_{\sigma}(\mathbf{r})$ and $V_{\pi}(\mathbf{r})$ are Slater-Koster functions V , which decay exponentially and were parameterized following previous ab initio works Lin ; Tang .

Diagonalization of the Bloch Hamiltonian results in a set of four-component Bloch eigenspinors $\hat{\psi}_{n,\mathbf{k}}(\mathbf{r})\equiv\langle\mathbf{r}|\hat{\psi}_{n,\mathbf{k}}\rangle$ that satisfy $\hat{\psi}_{n,\mathbf{k}}(\mathbf{r}+\mathbf{T})=\hat{\psi}_{n,\mathbf{k}}(\mathbf{r})\text{e}^{i\mathbf{k}\cdot\mathbf{T}}$ and correspond to the energy spectrum $\varepsilon_{n}(\mathbf{k})$. We calculate the bands for a small twist angle of $\theta=1.0845^{\circ}$ ($m=30$, $r=1$) near the experimental magic angle $\theta_{0}\sim 1.1^{\circ}$. At that angle, the Bloch Hamiltonian is a $N_{s}\times N_{s}$ matrix with $N_{s}=11164$ sites inside the Moire unit cell. The low energy bands ($n=1,\ldots,4$), shown in Fig. 2b, are four-fold degenerate at the $K$ points (excluding the spin). They have a two-fold degeneracy at the other two high symmetry points of the Brillouin zone, $\Gamma^{\prime}$ and $M^{\prime}$, where they open up a gap between particle and hole branches. At the $\Gamma^{\prime}$ point, the Bloch states have $C_{3}$ and $C_{2y}^{\prime}$ symmetry, which involves at $\pi$ rotation around the $y$ axis placed half-way between the two layers (shown in Fig. 1b). We also find numerically that all Bloch eigenspinors satisfy the time reversal symmetry (TRS) relation $\mathcal{T}\hat{\psi}_{n,\mathbf{k}}(\mathbf{r})=\hat{\psi}_{n,-\mathbf{k}}^{*}(\mathbf{r})$, with $\mathbf{k}$ measured from the center of the Moire Brillouin zone at $\Gamma^{\prime}$. The $K$ and $K^{\prime}$ points are hence related by TRS, and must have opposite $\pi$ Berry phases. This fact indicates that the Bloch states of the twisted structure do not suffer from Wannier obstructions Po , and hence could be reconstructed through a proper basis of Wannier states.

Wannier orbitals.$-$ From the Bloch states of the four low energy bands, one can extract the Wannier wave functions in the Moire unit cell,

[TABLE]

where $\mathbf{R}$ is the center of the Wannier orbitals and $U_{n\nu}(\mathbf{k})$ some $4\times 4$ unitary transformation. The four component Wannier spinors $\hat{W}_{\nu}(\mathbf{r}-\mathbf{R})\equiv\langle\mathbf{r}|\mathbf{R}\nu\rangle$ are not unique since adding a phase to the Bloch state $\text{e}^{-i\mathbf{k}\cdot\mathbf{r}}\hat{\psi}_{n\mathbf{k}}(\mathbf{r})$ corresponds to a new set of Wannier orbitals. We choose the set of maximally localized Wannier orbitals in finding the unitary transformation that minimizes their spread, $\Omega=\sum_{\nu}[\langle r^{2}\rangle_{\nu}-\langle r\rangle_{\nu}^{2}]$, with $\langle X\rangle_{\nu}\equiv\langle\mathbf{R}\nu|X|\mathbf{R}\nu\rangle$. The minimization was carried with the Wannier90 package Mostofi . The momentum space $\mathbf{k}$ mesh points are generated by the reciprocal supercell lattice vectors with 300 $\times$ 300 grid points using periodic boundary conditions, including all high symmetry points.

Following the symmetry arguments outlined in Ref. Kang , we perform the minimization of the spread enforcing the $C_{3}$ and $C_{2y}^{\prime}$ symmetry for the Bloch states around the $\Gamma^{\prime}$ points. Those two symmetries describe a $D_{3}$ point symmetry group, which is a local symmetry of the lattice at $AA$ site regions when the two graphene layers are rotated around a site Yuan , as depicted in Fig. 1b. In agreement with earlier results Kang ; Koshino , the Wannier functions that satisfy those symmetries have three sharp peaks centered around either the $AB$ or $BA$ sites, forming a honeycomb superlattice with two-fold degenerate orbitals per site, as shown in Fig. 3.

On a given Moire unit cell, we label the Wannier orbitals by the four-component spinors $\hat{W}_{\nu}=(w_{\nu,1},w_{\nu,2},w_{\nu,3},w_{\nu,4})^{T}$. Among the four orbitals, $\hat{W}_{\nu}(r-\mathbf{R}_{j})$, two are centered at $R_{j}\in AB$ sites and are eigenstates of the $C_{3}$ rotation operator, with eigenvalues $\epsilon=\text{e}^{2\pi i/3}$ and $\epsilon^{*}$. The other two are centered at $R_{j}\in BA$ sites and also have the same eigenvalues $\epsilon$ and $\epsilon^{*}$. From now on, we will label the Wannier orbital spinors based on their $C_{3}$ rotation eigenvalues, $C_{3}\hat{W}_{\alpha}(\mathbf{r}-\mathbf{R}_{j})=\text{e}^{\alpha 2\pi i/3}\hat{W}_{\alpha}(\mathbf{r}-\mathbf{R}_{j})$, with $\alpha=\pm$ and $R_{j}\in AB$ or $BA$. The two degenerate orbitals centered at a given superlattice site $\mathbf{R}_{j}$ are related by TRS, $\mathcal{T}\hat{W}_{\alpha}(\mathbf{r}-\mathbf{R}_{j})=\hat{W}_{-\alpha}(\mathbf{r}-\mathbf{R}_{j})$. Orbitals in NN superlattice sites $\mathbf{R}_{i}$ and $\mathbf{R}_{j}$ are related by the $C_{2}^{\prime}$ rotation, $C_{2}^{\prime}\hat{W}_{\alpha}(\mathbf{r}-\mathbf{R}_{i})=\hat{W}_{-\alpha}(\mathbf{r}-\mathbf{R}_{j})$.

*Tight binding Hamiltonian.$-$ *The effective lattice model of this problem can be constructed by rewriting the Bloch Hamiltonian (4) into a kinetic energy term of the form

[TABLE]

where $\mathbf{R}_{i}$ indexes the sites of the honeycomb superlattice, $\mathbf{R}_{ij}\equiv\mathbf{R}_{i}-\mathbf{R}_{j}$ and the $d_{\alpha}(\mathbf{R})$ annihilates an electron with orbital of type $\alpha$ and spin $\sigma$ at a given superlattice site. The hopping matrix elements between superlattice sites can be extracted from the matrix elements of Hamiltonian (2) in a basis of maximally localized Wannier functions,

[TABLE]

Due to the translational invariance of the superlattice, $t_{\alpha\beta}(\mathbf{R})=\langle 0,\alpha|\mathcal{H}|\mathbf{R},\beta\rangle$. For NN sites, we find that $|t_{\alpha\alpha}|(1)\approx 0.384$ meV whereas for $n$-th NN sites $t_{\alpha,-\alpha}(n)=0$. Hence, hopping between sites conserves the orbital pseudospin quantum number $\alpha=\pm$. $|t_{\alpha\alpha}|(n)$ has a non-trivial dependence with the distance between sites (see table I), in qualitative agreement with the findings of Ref. Kang for a significantly larger twist angle.

The Coulomb interactions between lattice sites can be written as $\mathcal{H}_{C}=\frac{1}{2}\int\text{d}\text{{r}}\text{d}\mathbf{r}^{\prime}\rho(\mathbf{r})\frac{e^{2}}{\kappa|\mathbf{r}-\mathbf{r}^{\prime}|}\rho(\mathbf{r}^{\prime})$, where $\rho(\mathbf{r})$ is the electron density and $\kappa\approx 5$ the dielectric constant of twisted bilayers encapsulated in boron nitride. We can rewrite this term in terms of $d_{\alpha,\sigma}$ operators by expressing the density $\rho(\mathbf{r})=\sum_{\sigma}\hat{\Psi}_{\sigma}^{\dagger}(\mathbf{r})\hat{\Psi}_{\alpha}(\mathbf{r})$ in terms of field operators $\hat{\Psi}_{\sigma}(\mathbf{r})=\sum_{\alpha,j}\hat{W}_{\alpha}(\mathbf{r}-\mathbf{R}_{j})d_{\alpha,\sigma}.$ The resulting Coulomb Hamiltonian has a direct term and also an exchange part, $\mathcal{H}_{C}=\mathcal{H}_{d}+\mathcal{H}_{e}$. The first term,

[TABLE]

with $n_{\alpha}(\mathbf{R})=\sum_{\sigma}d_{\alpha,\sigma}^{\dagger}(\mathbf{R})d_{\alpha,\sigma}(\mathbf{R})$ the density operator and repeated $\alpha,\,\beta$ indexes to be summed. The Coulomb coupling is cast as an overlap integral of Wannier orbital spinors, $V_{\alpha\beta}(\mathbf{R}_{ij})=\frac{1}{2}\int\text{d}\mathbf{r}\text{d}\mathbf{r}^{\prime}|\hat{W}_{\alpha}(\mathbf{r}_{i})|^{2}\frac{e^{2}}{\kappa|\mathbf{r}-\mathbf{r}^{\prime}|}|\hat{W}_{\beta}(\mathbf{r}_{j}^{\prime})|^{2},$ with $\mathbf{R}_{ij}=\mathbf{R}_{i}-\mathbf{R}_{j}$ and $\mathbf{r}_{j}\equiv\mathbf{r}-\mathbf{R}_{j}$. The exchange part is

[TABLE]

where

[TABLE]

is the exchange coupling between lattice sites. In general, we find that the combinations $J_{\alpha\beta,\beta\alpha}(\mathbf{R}_{ij})=J_{\alpha\beta,\alpha\beta}(\mathbf{R}_{ij})=0$ for $\alpha\neq\beta$, within the numerical precision. That includes the on site exchange (Hund’s coupling), which is zero due to the orthogonality between same site Wannier spinors Uchoa1 ; Koshino . From now on, we define the only non-zero combination $J_{\alpha\alpha,\beta\beta}\equiv J_{\alpha\beta}$.

The numerical values of the hopping energy, Coulomb interaction and the exchange interaction for $n$-th NNs, is shown in table I, which is the first main result of the paper. We find the on-site Hubbard $U_{\alpha\beta}\equiv V_{\alpha\beta}(0)=21.2$ meV, which is much larger than the first NN hopping $t(1)$, and hence the ratio $U/t(1)\sim 55$ falls comfortably in the realm of the Mott regime.

The exchange interaction for first NN sites $(n=1)$ is $|J_{\alpha\beta}(1)|\approx 5$ meV. In general, the diagonal terms $J_{\alpha\alpha}(n)>0$ are positive definite, whereas the off diagonal ones can be either positive or negative, $J_{\alpha,-\alpha}(n)\approx\pm J_{\alpha\alpha}(n)$, with $+$ ($-$) for $i,\,j$ sites in the same (opposite) sublattice, as shown in table I. For sites in the same sublattice, the fact that $J_{\alpha\beta}(n)\approx J_{\alpha\alpha}(n)>0$ is the same for all four combinations of $\alpha,\beta=\pm$ indexes hints at an emergent SU(4) symmetry between spin and orbital degrees of freedom at quarter filling. For sites in opposite sublattices, the exchange interaction has SU(2) symmetry in the spin. It has also both ferro $(J_{\alpha\alpha}>0$) and antiferromagnetic $(J_{\alpha,-\alpha}<0)$ correlations in the orbital sector, depending on the orientation of the pseudospins.

Since Hund’s coupling is zero, at quarter filling the lower flat bands are in the unitary limit Coqblin , with each Moire superlattice site $\mathbf{R}_{j}$ being singly occupied and having a well defined spin $\sigma$ and orbital quantum number $\alpha=\pm$. Mapping the exchange term in terms of spin $\mathbf{S}_{i}=\frac{1}{2}d_{\alpha,\sigma}^{\dagger}(\mathbf{R}_{i})\vec{\sigma}_{\sigma\sigma^{\prime}}d_{\alpha,\sigma^{\prime}}(\mathbf{R}_{i})$ and pseudospin $\bm{\tau}_{i}=\frac{1}{2}d_{\alpha,\sigma}^{\dagger}(\mathbf{R}_{i})\vec{\sigma}_{\alpha\beta}d_{\beta,\sigma}(\mathbf{R}_{i})$ operators, the result is the ferromagnetic exchange interaction announced in Eq. (1), with $J_{ij}\equiv J_{\alpha\alpha}(n)>0$ SM1 . This Hamiltonian favors ferromagnetic alignment of the spins at zero temperature ($T=0$). In the orbital sector different states are possible, including canted magnetism with ferromagnetic order in the pseudospin $\tau^{z}$ component, accompanied by staggered (antiferromagnetic) order in the transverse, $\bm{\tau}^{\perp}$ direction.

The superexchange interaction follows from second order perturbation theory in the hopping energy kk ; Vanderbos and has the same form as the exchange term in Eq. (1) for $\eta_{ij}=-1$ Uchoa1 . The superexchange term has SU(4) symmetry and favors antiferromagnetic alignment between nearest neighbor sites due to Pauli principle. It’s coupling $J\to-t^{2}/U\approx-0.01$ meV is very small compared to the exchange one, and can be safely igonored.

Ferromagnetic Mott state.$-$ Mott-Hubbard insulators have strongly localized states and are known to be overwhelmingly antiferromagnetic due to strong superexchange interactions ($t^{2}/U\gg J$) khomskii ; anderson ; anderson-2 . Ferromagnetism occurs mostly either in metallic systems or in metallic bands hybridized with localized moments via the Anderson impurity mechanism anderson ; anderson-2 ; moriya . Within the Hubbard model framework, the only credible mechanism for spin ferromagnetism exists for multi-orbital systems in the context of the Kugel-Khomskii model khomskii ; kk , where superexchange can become effectively ferromagnetic in the presence of staggered orbital ordering.

We conjecture that the flat bands in twisted graphene bilayers are in a way intermediate between ferromagnetic bad metals and antiferromagnetic Mott-Hubbard insulators. Due to the exotic shape of the Wannier orbitals, the hierarchy between hopping, direct exchange and the local Hubbard interaction, $t\ll J\ll U,$ leads to an anomalously small superexchange. In the charge sector the Mott gap is also anomalously small, $\Delta\sim 0.3\text{ meV}\ll W\ll U$, where $W\sim 5$ meV is the bandwidth, and the system undergoes an insulator-metal transition at $T\approx 4$K Cao1 .

In spite of the fact that $U/t$ is large, the strong overlap between the orbitals found in the non-interacting theory suggests that the system is potentially close to an insulator-metal transition imada due to a charge fluctuation mechanism which presently is not well-understood Pizarro ; Guinea . Nevertheless, the effective spin model we propose in this work should not depend on the details of this mechanism, as long as the system remains quarter filled and does not undergo a charge-ordering transition (potentially accompanied by dimerization) due to Coulomb interactions. In carbon lattices, which are notoriously stiff Lee , charge density wave instabilities are hindered by the high elastic energy cost for the system to deform the lattice and restore charge neutrality.

Experimental signatures.$-$ Since the honeycomb superlattice is not frustrated, it will exhibit ferromagnetic spin order at $T=0$ in the universality class of the ferromagnetic (spin $S$) Heisenberg model. It is well known that the magnetization $M$, correlation length $\xi$ and the spin susceptibility $\chi$ exhibit peculiar features in two dimensions, since for any $T\neq 0$ the system is disordered, with zero Curie temperature. The model has been extensively studied both in zero and finite external magnetic field $H$ on various lattices kopietz ; antsygina ; junger ; note2 . At finite field $H\neq 0$, $M(H)$ is finite and strongly temperature dependent. In the regime $T/J\ll 1$, which can take place for $T\approx 2$K (where $T/J\approx 1/25$), a weak magnetic field of $H\approx 0.2$T (i.e. $H/J\approx 1/250)$ already provides nearly maximum magnetization antsygina ; junger . The susceptibility $\chi(H)$ is zero for $T=0$ and $H\neq 0$ and exhibits a characteristic finite-temperature peak at $T=T_{\chi}$ which scales in a well-defined way with external field.

It has been established experimentally that doping away from the Mott insulating phase leads to metallic (and even superconducting) behavior Cao1 ; Cao2 . Therefore the structure of the ground state and excitation spectrum of this unconventional metallic state is of great experimental and theoretical interest. A profound new feature has emerged at finite magnetic field, which persists both in weak (Shubnikov-de Haas oscillations) and strong field limits (Quantum Hall effect), for hole doping Cao1 ; Cao2 . Those measurements suggest a small Fermi surface that develops from doping the correlated insulating phase, accompanied by a possible symmetry breaking of yet unknown origin. The resulting state has a fermionic degeneracy of 2, indicating a reduction of the original four-fold band degeneracy by a factor of 2.

This behavior is consistent with the system being in the proximity to a ferromagnetic Mott state, in which the spins align when nudged by an infinitesimally weak field. At the same time, any long-range order in the orbital sector is expected to be much more fragile and disappear quickly due to charge disorder and motion of holes in the metallic state. Therefore we conjecture that in the weak field limit, the ground state emerging from doping the ferromagnetic insulator would be a ferromagnetic, spin-polarized, strongly-correlated metal, with the orbital pseudospin symmetry preserved.

*Conclusions.$-$ *We have derived the effective spin model that describes the Mott phase of twisted graphene bilayers at the magic angle. After calculating the maximally localized Wannier wavefunctions from the lattice, we propose that the system forms a novel ferromagnetic Mott state at quarter filling, with clear experimental signatures.

Acknowledgements.$-$ BU Acknowledges P. Jarillo-Harillo, T. Senthil, and K. Beach for discussions. K.S. and B. U. acknowledge NSF CAREER grant No. DMR-1352604 for support.

Note added.$-$ After the submission of this work, we became aware of Ref. Oskar , which also found a ferromagnetic ground state using different arguments. We thank O. Vafek for pointing it out.