Gonzalez-Arraga L.A., Lado J.L., Guinea F., San-Jose P., Gonzalez-Arraga Luis A. 1 Lado J.?L. 2 Guinea Francisco 1,3 San-Jose Pablo 4 1 IMDEA Nanociencia , Calle de Faraday, 9, Cantoblanco, 28049 Madrid, Spain QuantaLab, 2 International Iberian Nanotechnology Laboratory (INL) , Avenida Mestre Jose Veiga, 4715-330 Braga, Portugal School of Physics and Astronomy, 3 University of Manchester , Oxford Road, Manchester M13 9PL, United Kingdom 4 Instituto de Ciencia de Materiales de Madrid (ICMM-CSIC) , Cantoblanco, 28049 Madrid, Spain 5 September 2017 8 September 2017 119 10 107201 27 March 2017 6 July 2017 © 2017 American Physical Society 2017 American Physical Society Twisted graphene bilayers develop highly localized states around A A -stacked regions for small twist angles. We show that interaction effects may induce either an antiferromagnetic or a ferromagnetic (FM) polarization of said regions, depending on the electrical bias between layers. Remarkably, FM-polarized A A regions under bias develop spiral magnetic ordering, with a relative 120° misalignment between neighboring regions due to a frustrated antiferromagnetic exchange. This remarkable spiral magnetism emerges naturally without the need of spin-orbit coupling, and competes with the more conventional lattice-antiferromagnetic instability, which interestingly develops at smaller bias under weaker interactions than in monolayer graphene, due to Fermi velocity suppression. This rich and electrically controllable magnetism could turn twisted bilayer graphene into an ideal system to study frustrated magnetism in two dimensions. Marie-Curie-ITN 607904-SPINOGRAPH Ministerio de Economía y Competitividad http://dx.doi.org/10.13039/501100003329 Ministry of Economy and Competitiveness MINECO http://sws.geonames.org/2510769/ FIS2015-65706-P Ramon y Cajal RYC-2013-14645 Magnetism in 2D electronic systems is known to present a very different phenomenology from its three-dimensional counterpart due to the reduced dimensionality and the increased importance of fluctuations. Striking examples are the impossibility of establishing long range magnetic order in a 2D system without magnetic anisotropy [1] or the emergence of unique finite-temperature phase transitions that are controlled by the proliferation of topological magnetic defects [2] . In the presence of magnetic frustration, in, e.g., Kagome [3,4] or triangular lattices [5–8] , 2D magnetism may also lead to the formation of remarkable quantum spin-liquid phases [3,9,10] . The properties of these states remain under active investigation, and have recently been shown to develop exotic properties, such as fractionalized excitations [11] , long-range quantum entanglement of their ground state [12,13] , topologically protected transport channels [14] , or even high- T C superconductivity upon doping [4,15,16] . The importance of 2D magnetism extends also beyond fundamental physics into applied fields. One notable example is data storage technologies. Recent advances in this field are putting great pressure on the magnetic memory industry to develop solutions that may remain competitive in speed and data densities against new emerging platforms. Magnetic 2D materials are thus in demand as a possible way forward [17] . Of particular interest for applications in general are 2D crystals and van der Waals heterostructures. These materials have already demonstrated a great potential for a wide variety of applications, most notably nanoelectronics and optoelectronics [18–20] . Some of them have been shown to exhibit considerable tunability through doping, gating, stacking, and strain. Unfortunately, very few 2D crystals have been found to exhibit intrinsic magnetism [21,22] , let alone magnetic frustration and potential spin-liquid phases. In this work we predict that twisted graphene bilayers could be a notable exception, realizing a peculiar magnetism on an effective triangular superlattice, and with exchange interactions that may be tuned by an external electric bias. We show that, at a mean-field level, spontaneous magnetization of two different types may develop for small enough twist angles ? ? 2 ° as a consequence of the moiré pattern in the system. This effect is a consequence of the high local density of states generated close to neutrality at moiré regions with A A stacking, triggering a Stoner instability when electrons interact. The local order is localized at A A regions but may be either antiferromagnetic (AFM) or ferromagnetic (FM). The two magnetic orders can be switched electrically by applying a voltage bias between layers. Interestingly, the relative ordering between different A A regions in the FM ground state is predicted to be spiral, despite the system possessing negligible spin-orbit coupling. This type of magnetism combines a set of unique features: electric tunability, magnetic frustration, the interplay of two switchable magnetic phases with zero net magnetization, spatial localization of magnetic moments, and an adjustable period of the magnetic superlattice. Finally, we show that our mean-field treatment allows us to cast the system into an effective spin Hamiltonian that could be tackled beyond the mean-field level to evaluate the effects of spin fluctuations. The type of frustrated spin Hamiltonian obtained suggests that twisted graphene bilayers should be a prime playground for studies of spin-liquid phases. We discuss some of these possibilities in our concluding remarks. Description of the system.— Twisted graphene bilayers are characterized by a relative rotation angle ? between the two layers [23] . The rotation produces a modulation of the relative stacking at each point, following a moiré pattern of period L M ? a 0 / ? at small ? , where a 0 = 0.24 ? ? nm is graphene’s lattice constant [24] . The stacking smoothly interpolates between three basic types, A A (perfect local alignment of the two lattices), and A B or B A (Bernal stackings related by point inversion) [25] . The stacking modulation leads to a spatially varying coupling between layers. This results in a remarkable electronic reconstruction [26,27] , particularly at small angles ? ? 1 ° – 2 ° [28,29] , for which the interlayer coupling ? 1 ? 0.3 ? ? eV exceeds the moiré energy scale ? M = ? v F ? K [here, ? K = 4 ? / ( 3 L M ) is the rotation-induced wave vector shift between the Dirac points in the two layers, and v F ? 10 6 ? ? m / s is the monolayer Fermi velocity]. It was shown [24,29–33] that in such a regime the Fermi velocity of the bilayer becomes strongly suppressed, and the local density of states close to neutrality becomes dominated by quasilocalized states in the A A regions [28] . The confinement of these states is further enhanced by an interlayer bias V b , which effectively depletes the A B and B A regions due to the opening of a local gap [34,35] . At sufficiently small angles this was also shown to result in the formation of a network of helical valley currents flowing along the boundaries of depleted A B and B A regions [36] . The quasilocalized A A states form a weakly coupled triangular superlattice of period L M , analogous to a network of quantum dots. Each A A 'dot' has space for eight degenerate electrons, due to the sublattice, layer, and spin degrees of freedom. A plot of their spatial distribution under zero and large bias V b = 300 ? ? meV is shown in Figs. 1(a) and 1(b) , respectively. These A A states form an almost flat band at zero energy [37] , see Figs. 1(c) and 1(d) , which gives rise to a zero-energy peak in the DOS. The small but finite width of this zero-energy A A resonance represents the residual coupling between adjacent A A dots due to their finite overlap. A comparison of Figs. 1(a) and 1(b) shows that a finite interlayer bias leads to a suppression of said overlap and a depletion of the intervening A B and B A regions, as described above. The electronic structure presented here was computed using the tight-binding approach described in the Supplemental Material [38] , which includes a scaling approximation that allows the accurate and efficient computation of the low-energy band structure in low-angle twisted bilayers [compare the solid and dashed curves in Figs 1(c) and 1(d) ]. Our scaling approach makes the problem much more tractable computationally, which is a considerable advantage when dealing with the interaction effects, discussed below. 1 10.1103/PhysRevLett.119.107201.f1 FIG. 1. Zero-energy local density of states in real space (a),(b), band structure (c),(d), and density of states (e),(f) for a ? = 1.5 ° twisted graphene bilayer. The left column has no interlayer bias, and the right column has a bias V b = 300 ? ? meV . This enhances the localization of the A A quasibound states, red in (a) and (b). The said states arise from almost flat subbands at zero energy, which show up as large DOS peaks in (e) and (f). The solid (dashed) lines in (c) and (d) correspond to a scaled (unscaled) tight-binding model, see the main text. Moiré-induced magnetism.— It is known that in the presence of sufficiently strong electronic interactions, a honeycomb tight-binding lattice may develop a variety of ground states with spontaneously broken symmetry [42–46] . The simplest one is the lattice antiferromagnetic phase in the honeycomb Hubbard model. The Hubbard model is a simple description relevant to monolayer graphene with strongly screened interactions (the screening may arise intrinsically at high doping or, e.g., due to a metallic environment). Above a critical value of the Hubbard coupling, U > U c ( 0 ) ? 5.7 ? ? eV (value within the mean-field), the system favors a ground state in which the two sublattices are spin polarized antiferromagnetically. This is known as lattice-AFM (or Néel) order. In the absence of adsorbates [47] , edges [48] , vacancies [49] , or magnetic flux [50] isolated graphene monolayers, with their vanishing density of states at low energies, are known experimentally not to suffer any interaction-induced magnetic instability. In contrast, Bernal ( ? = 0 ) bilayer graphene and A B C trilayer graphene have been suggested [51–54] to develop magnetic order, due to their finite low-energy density of states, although some controversy remains [55–60] . Twisted graphene bilayers at small angles exhibit an even stronger enhancement of the low-energy density of states associated with A A confinement and the formation of quasiflat bands. It is thus natural to expect some form of interaction-induced instability in this system with realistic interactions, despite the lack of magnetism in the monolayer [61] . By analyzing the Hubbard model in twisted bilayers we now explore this possibility, and describe the different magnetic orders that emerge in the U , V b parameter space. We consider the Hubbard model in a low angle ? ? 1.5 ° twisted bilayer for a moderate [62] value of U = 3.7 , quite below the monolayer lattice-AFM critical interaction U c ( 0 ) . We use a self-consistent mean-field approximation to compute the system’s ground state, and use the same parameters of Fig. 1 . Self-consistency involves the iterative computation of charge and spin density on the moiré supercell, integrated over Bloch momenta, see the Supplemental Material [38] for details. Since U is repulsive we neglect superconducting symmetry breaking, and concentrate on arbitrary normal solutions instead [63] . In Fig. 2 we show the resulting real-space distribution of the ground-state spin polarization M ( r ? ) of the converged solution. The top and bottom rows correspond, respectively, to the lattice-FM and lattice-AFM components M A + M B and M A - M B , where the polarization density is defined as M ? = ? ? ? n ? ? ? ( r ? ) - n ? ? ? ( r ? ) ? . Here, ? = A , B are the two sublattices and ? = ± are the two layers. 2 10.1103/PhysRevLett.119.107201.f2 FIG. 2. Spatial distribution of the magnetic moment in the ground state of an interacting twisted bilayer with Hubbard U = 3.7 ? ? eV . In the first row (a),(b) we show the ferromagnetic component of the two sublattices, M A + M B , in units of electrons per (monolayer) unit cell, both for zero interlayer bias V b = 0 (a) and V b = 200 ? ? meV (b). Analogous plots of the lattice-AFM component M A - M B are shown in (c) and (d). The scale in all color bars is expressed in units of one electron spins per supercell. Panels (e) and (f) show the variation of the total electronic energy per supercell as a function of the angle ? M between polarizations of adjacent A A regions, indicating parallel alignment of the lattice-AFM order (e), and a spiral misalignment of 120° for the lattice-FM case (f). We obtain two distinct solutions for the magnetization, depending on the interlayer bias V b . At small interlayer bias and for the chosen U = 3.7 ? ? eV we see that the ferromagnetic polarization [Fig. 2(a) ] is small and collinear, and spatially integrates to zero. Thus, the unbiased bilayer remains nonferromagnetic in the small V b case. However, the lattice-AFM component of the polarization, Fig. 2(c) , is large and integrates to a nonzero value of around 0.5 electron spins per unit cell. This is the analogue of the monolayer lattice-AFM phase, with two important differences. On the one hand, we find that the lattice-AFM density is strongly concentrated at the A A regions instead of being spatially uniform like in the monolayer. On the other hand the lattice-AFM ground state is found to arise already for U ? 2 ? ? eV , i.e., for much weaker interactions than in the monolayer. The reason for the reduction of U c can be traced to the suppression of the Fermi velocity v F at small twist angles [29,32] , which controls the critical U for the lattice-AFM instability. The dependence of U c and v F as a function of angle ? is shown in Fig. 3(a) . This result already points to strong magnetic instabilities of twisted graphene bilayers as the angle falls below the 1°–2° threshold. 3 10.1103/PhysRevLett.119.107201.f3 FIG. 3. (a) Critical value U c of the Hubbard U beyond which the twisted bilayer develops lattice-AFM order at the mean-field level. The red dots show U c as a function of the twist angle ? , and the dashed line shows the corresponding Fermi velocity at the Dirac point, normalized to the monolayer value v F 0 . At high twist angles both U c and v F converge to the monolayer values, while they become strongly suppressed at smaller angles. (b) Phase diagram for the ground state magnetic order in a ? = 1.5 ° twisted bilayer as a function of Hubbard U and interlayer bias V b . The blue and red regions denote the spatial integral of the lattice-AFM and spiral-FM polarizations, respectively, while the yellow region is nonmagnetic. Under a large electric bias between layers, the ground state magnetization for the same U is dramatically different, see Figs. 2(b) and 2(d) . In this case, the lattice-AFM polarization, Fig. 2(d) , is strongly suppressed and integrates to zero spatially, while the lattice-FM component, Fig. 2(b) , becomes large around the A A regions, and integrates to a finite value of approximately four electron spins per moiré supercell. The A A regions are thus found to become ferromagnetic under sufficient interlayer bias. This type of magnetic order is the result of the increased confinement of A A states at high V b , and can be interpreted as an instance of flat-band ferromagnetism driven by the Stoner mechanism. The lattice-AFM and lattice-FM states are also different when comparing the relative orientations of neighboring A A regions. By computing the total energy per supercell in each case as a function of the polarization angle ? M between adjacent regions [Figs. 2(e) and 2(f) ], we find that the energy is minimized for ? M = 0 ° in the lattice-AFM case (parallel alignment), but for ? M = 12 0 ° in the lattice-FM case (spiralling polarization). The equilibrium polarization is depicted by white arrows in Figs. 2(c) and 2(b) . The depth of the energy minimum, ranging from ? 2 – 100 ? ? K in our simulations, represents the effective exchange coupling of neighboring A A regions, which is ferromagnetic for lattice-AFM states and antiferromagnetic for lattice-FM states (see the Supplemental Material [38] for the next-nearest neighbor exchange). In the lattice-FM phase, which from now on we denote the spiral-FM phase, the spiral order arises as a result of the triangular symmetry of A A regions that frustrates a globally antiferromagnetic A A alignment. The same spiral order has been described in studies of the Hubbard model in the triangular lattice. It is a rather remarkable magnetic state, as the polarization at different points becomes noncollinear [7,64,65] despite the complete absence of spin-orbit coupling in the system. To better understand the onset of the spiral magnetism, we have computed the integrated FM and AFM polarization across the U , V b plane. We find first-order phase transitions separating the two types of ground states. The result is shown in Fig. 3 . The regions in red and blue denote, respectively, a finite spatial integral of the ferro M A + M B and lattice-AFM M A - M B polarizations. It is apparent that an electric interlayer bias of around 120 meV is able to switch between the lattice-AFM and spiral-FM orders for values of U between 2 and 3 eV. The precise thresholds for such electric switching of magnetic order depend on the specific twist angle and on other details not considered in this work (longer-range interactions, spontaneous deformations, or interlayer screening), although our simulations suggest they are likely within reach of current experiments for sufficiently small ? . Our mean-field analysis neglects thermal and quantum spin fluctuations around the mean-field solution. Thermal spin excitations in the magnetically isotropic case under study (from gapless Goldstone modes) are expected to destroy long-range spiral order, which then survives only locally, in keeping with the Mermin-Wagner theorem [1] . Breaking the magnetic isotropy (by allowing for a hard magnetic axis due to, e.g., spin-orbit coupling or coupling to a suitable magnetic substrate) gaps the Goldstone modes and stabilizes the mean-field solution. Otherwise, even at zero temperature, quantum spin fluctuations are known to produce spin-liquid-like ground states [5–8] . An efficient way to explore such nontrivial effects in this moiré system is to cast our mean-field results into an effective spin Hamiltonian on the triangular A A moiré pattern, which could be tackled using more sophisticated approaches (e.g., matrix-product states). The procedure is described in the Supplemental Material [38] . Conclusion.— For a long time unmodified graphene was thought to be relatively uninteresting from the point of view of magnetism. Twisted graphene bilayers, however, could prove to be a surprisingly rich playground for nontrivial magnetic phases. We have shown that two different types of mean-field magnetic solutions arise in twisted graphene bilayers at small angles. The two types of magnetic order, lattice antiferromagnetism and spiral ferromagnetism, are both concentrated at A A -stacked regions. The spiral-FM phase is favored over the lattice-AFM phase when applying a sufficient electric bias between layers. This phase constitutes a form of electrically controllable, noncollinear, and spatially nonuniform magnetism in a material with a negligible spin-orbit coupling. This possibility is of fundamental interest, as it realizes electrically tunable 2D magnetism on a triangular superlattice, a suitable platform to explore spin-liquid phases. Indeed, it is known that next-nearest neighbor interactions in a magnetic triangular lattice should transform spiral order into a spin-liquid phase [5–8] , as long as the system remains magnetically isotropic. Moreover, in the spin-liquid state, electronic doping can give rise to high T C superconductivity [66,67] . The possibility of modifying the electronic filling of our emergent frustrated triangular lattice by means of an electric gate offers a unique platform to realize this possibility, avoiding the detrimental effects of chemical doping in conventional compounds [68] . While the above is highly speculative at this point and would require a careful nonperturbative analysis of our effective spin Hamiltonian, it highlights the interesting fundamental possibilities afforded by the rich magnetic phase diagram of twisted graphene bilayers. We acknowledge financial support from the Marie-Curie-Initial Training Networks (ITN) program through Grant No. 607904-SPINOGRAPH, and the Spanish Ministry of Economy and Competitiveness through Grants No. FIS2015-65706-P (MINECO/FEDER) and No. RYC-2013-14645 (Ramon y Cajal program). L.?A.?G.-A. is grateful for the hospitality of the Applied Physics Department in the University of Alicante and to N. Garcia for useful discussions. We specially thank J. Fernandez Rossier for his help settling the environment and the initial idea for this work. [1] 1 N.?D. Mermin and H. Wagner , Phys. Rev. Lett. 17 , 1133 ( 1966 ). PRLTAO 0031-9007 10.1103/PhysRevLett.17.1133 [2] 2 D.?R. Nelson and J.?M. Kosterlitz , Phys. Rev. Lett. 39 , 1201 ( 1977 ). PRLTAO 0031-9007 10.1103/PhysRevLett.39.1201 [3] 3 M. Fu , T. Imai , T.-H. Han , and Y.?S. Lee , Science 350 , 655 ( 2015 ). SCIEAS 0036-8075 10.1126/science.aab2120 [4] 4 S.-H. Lee , H. Kikuchi , Y. Qiu , B. Lake , Q. Huang , K. Habicht , and K. Kiefer , Nat. Mater. 6 , 853 ( 2007 ). NMAACR 1476-1122 10.1038/nmat1986 [5] 5 T. Isono , H. Kamo , A. Ueda , K. Takahashi , M. Kimata , H. Tajima , S. Tsuchiya , T. Terashima , S. Uji , and H. Mori , Phys. Rev. Lett. 112 , 177201 ( 2014 ). PRLTAO 0031-9007 10.1103/PhysRevLett.112.177201 [6] 6 L. Seabra , T. Momoi , P. Sindzingre , and N. Shannon , Phys. Rev. B 84 , 214418 ( 2011 ). PRBMDO 1098-0121 10.1103/PhysRevB.84.214418 [7] 7 W.-J. Hu , S.-S. Gong , W. Zhu , and D.?N. Sheng , Phys. Rev. B 92 , 140403 ( 2015 ). PRBMDO 1098-0121 10.1103/PhysRevB.92.140403 [8] 8 Z. Zhu and S.?R. White , Phys. Rev. B 92 , 041105 ( 2015 ). PRBMDO 1098-0121 10.1103/PhysRevB.92.041105 [9] 9 L. Savary and L. Balents , Rep. Prog. Phys. 80 , 016502 ( 2017 ). RPPHAG 0034-4885 10.1088/0034-4885/80/1/016502 [10] 10 Y. Xu , J. Zhang , Y.?S. Li , Y.?J. Yu , X.?C. Hong , Q.?M. Zhang , and S.?Y. Li , Phys. Rev. Lett. 117 , 267202 ( 2016 ). PRLTAO 0031-9007 10.1103/PhysRevLett.117.267202 [11] 11 T.-H. Han , J.?S. Helton , S. Chu , D.?G. Nocera , J.?A. Rodriguez-Rivera , C. Broholm , and Y.?S. Lee , Nature (London) 492 , 406 ( 2012 ). NATUAS 0028-0836 10.1038/nature11659 [12] 12 T. Grover , Y. Zhang , and A. Vishwanath , New J. Phys. 15 , 025002 ( 2013 ). NJOPFM 1367-2630 10.1088/1367-2630/15/2/025002 [13] 13 M. Pretko and T. Senthil , Phys. Rev. B 94 , 125112 ( 2016 ). PRBMDO 2469-9950 10.1103/PhysRevB.94.125112 [14] 14 N.?Y. Yao , C.?R. Laumann , A.?V. Gorshkov , H. Weimer , L. Jiang , J.?I. Cirac , P. Zoller , and M.?D. Lukin , Nat. Commun. 4 , 1585 ( 2013 ). NCAOBW 2041-1723 10.1038/ncomms2531 [15] 15 Z.?A. Kelly , M.?J. Gallagher , and T.?M. McQueen , Phys. Rev. X 6 , 041007 ( 2016 ). PRXHAE 2160-3308 10.1103/PhysRevX.6.041007 [16] 16 P.?W. Anderson , Mater. Res. Bull. 8 , 153 ( 1973 ). MRBUAC 0025-5408 10.1016/0025-5408(73)90167-0 [17] 17 X. Wang , K. Du , Y.?Y.?F. Liu , P. Hu , J. Zhang , Q. Zhang , M.?H.?S. Owen , X. Lu , C.?K. Gan , P. Sengupta , C. Kloc , and Q. Xiong , 2D Mater. 3 , 031009 ( 2016 ). DMATB7 2053-1583 10.1088/2053-1583/3/3/031009 [18] 18 K.?S. Novoselov , A. Mishchenko , A. Carvalho , and A.?H. Castro Neto , Science 353 , aac 9439 ( 2016 ). SCIEAS 0036-8075 10.1126/science.aac9439 [19] 19 A. Castellanos-Gomez , Nat. Photonics 10 , 202 ( 2016 ). NPAHBY 1749-4885 10.1038/nphoton.2016.53 [20] 20 J. Quereda , P. San-Jose , V. Parente , L. Vaquero-Garzon , A.?J. Molina-Mendoza , N. Agraït , G. Rubio-Bollinger , F. Guinea , R. Roldán , and A. Castellanos-Gomez , Nano Lett. 16 , 2931 ( 2016 ). NALEFD 1530-6984 10.1021/acs.nanolett.5b04670 [21] 21 B. Huang , G. Clark , E. Navarro-Moratalla , D.?R. Klein , R. Cheng , K.?L. Seyler , D. Zhong , E. Schmidgall , M.?A. McGuire , D.?H. Cobden , W. Yao , D. Xiao , P. Jarillo-Herrero , and X. Xu , Nature (London) 546 , 270 ( 2017 ). NATUAS 0028-0836 10.1038/nature22391 [22] 22 C. Gong , L. Li , Z. Li , H. Ji , A. Stern , Y. Xia , T. Cao , W. Bao , C. Wang , Y. Wang , Z.?Q. Qiu , R.?J. Cava , S.?G. Louie , J. Xia , and X. Zhang , Nature (London) 546 , 265 ( 2017 ). NATUAS 0028-0836 10.1038/nature22060 [23] 23 J.?M.?B. Lopes dos Santos , N.?M.?R. Peres , and A.?H. Castro Neto , Phys. Rev. Lett. 99 , 256802 ( 2007 ). PRLTAO 0031-9007 10.1103/PhysRevLett.99.256802 [24] 24 J.?M.?B. Lopes dos Santos , N.?M.?R. Peres , and A.?H. Castro Neto , Phys. Rev. B 86 , 155449 ( 2012 ). PRBMDO 1098-0121 10.1103/PhysRevB.86.155449 [25] 25 J.?S. Alden , A.?W. Tsen , P.?Y. Huang , R. Hovden , L. Brown , J. Park , D.?A. Muller , and P.?L. McEuen , Proc. Natl. Acad. Sci. U.S.A. 110 , 11256 ( 2013 ). PNASA6 0027-8424 10.1073/pnas.1309394110 [26] 26 G. Li , A. Luican , J.?M.?B. Lopes dos Santos , A.?H. Castro Neto , A. Reina , J. Kong , and E.?Y. Andrei , Nat. Phys. 6 , 109 ( 2010 ). NPAHAX 1745-2473 10.1038/nphys1463 [27] 27 A. Rozhkov , A. Sboychakov , A. Rakhmanov , and F. Nori , Phys. Rep. 648 , 1 ( 2016 ) PRPLCM 0370-1573 10.1016/j.physrep.2016.07.003 [28] 28 G.?T. de Laissardière , D. Mayou , and L. Magaud , Nano Lett. 10 , 804 ( 2010 ). NALEFD 1530-6984 10.1021/nl902948m [29] 29 R. Bistritzer and A.?H. MacDonald , Proc. Natl. Acad. Sci. U.S.A. 108 , 12233 ( 2011 ). PNASA6 0027-8424 10.1073/pnas.1108174108 [30] 30 A. Luican , G. Li , A. Reina , J. Kong , R.?R. Nair , K.?S. Novoselov , A.?K. Geim , and E.?Y. Andrei , Phys. Rev. Lett. 106 , 126802 ( 2011 ). PRLTAO 0031-9007 10.1103/PhysRevLett.106.126802 [31] 31 P. San-Jose , J. González , and F. Guinea , Phys. Rev. Lett. 108 , 216802 ( 2012 ). PRLTAO 0031-9007 10.1103/PhysRevLett.108.216802 [32] 32 G. Trambly de Laissardiere , D. Mayou , and L. Magaud , Phys. Rev. B 86 , 125413 ( 2012 ). PRBMDO 1098-0121 10.1103/PhysRevB.86.125413 [33] 33 A.?O. Sboychakov , A.?L. Rakhmanov , A.?V. Rozhkov , and F. Nori , Phys. Rev. B 92 , 075402 ( 2015 ). PRBMDO 1098-0121 10.1103/PhysRevB.92.075402 [34] 34 K.?F. Mak , C.?H. Lui , J. Shan , and T.?F. Heinz , Phys. Rev. Lett. 102 , 256405 ( 2009 ). PRLTAO 0031-9007 10.1103/PhysRevLett.102.256405 [35] 35 Y. Zhang , T.-T. Tang , C. Girit , Z. Hao , M.?C. Martin , A. Zettl , M.?F. Crommie , Y.?R. Shen , and F. Wang , Nature (London) 459 , 820 ( 2009 ). NATUAS 0028-0836 10.1038/nature08105 [36] 36 P. San-Jose and E. Prada , Phys. Rev. B 88 , 121408 ( 2013 ). PRBMDO 1098-0121 10.1103/PhysRevB.88.121408 [37] 37 E. Suárez Morell , J.?D. Correa , P. Vargas , M. Pacheco , and Z. Barticevic , Phys. Rev. B 82 , 121407 ( 2010 ). PRBMDO 1098-0121 10.1103/PhysRevB.82.121407 [38] 38 See Supplemental Material at http://link.aps.org/supplemental/10.1103/PhysRevLett.119.107201 for details on modeling and effective spin Hamiltonians, which includes Refs. [23,24,28,33,37,39–41]. [39] 39 P. Moon and M. Koshino , Phys. Rev. B 87 , 205404 ( 2013 ). PRBMDO 1098-0121 10.1103/PhysRevB.87.205404 [40] 40 P. Moon and M. Koshino , Phys. Rev. B 90 , 155406 ( 2014 ). PRBMDO 1098-0121 10.1103/PhysRevB.90.155406 [41] 41 R. Shankar , Rev. Mod. Phys. 66 , 129 ( 1994 ). RMPHAT 0034-6861 10.1103/RevModPhys.66.129 [42] 42 S. Sorella and E. Tosatti , Europhys. Lett. 19 , 699 ( 1992 ). EULEEJ 0295-5075 10.1209/0295-5075/19/8/007 [43] 43 I.?F. Herbut , Phys. Rev. Lett. 97 , 146401 ( 2006 ). PRLTAO 0031-9007 10.1103/PhysRevLett.97.146401 [44] 44 S. Sorella , Y. Otsuka , and S. Yunoki , Sci. Rep. 2 , 992 ( 2012 ). SRCEC3 2045-2322 10.1038/srep00992 . [45] 45 F.?F. Assaad and I.?F. Herbut , Phys. Rev. X 3 , 031010 ( 2013 ). PRXHAE 2160-3308 10.1103/PhysRevX.3.031010 [46] 46 N.?A. García-Martínez , A.?G. Grushin , T. Neupert , B. Valenzuela , and E.?V. Castro , Phys. Rev. B 88 , 245123 ( 2013 ). PRBMDO 1098-0121 10.1103/PhysRevB.88.245123 [47] 47 H. González-Herrero , J. Gómez-Rodríguez , P. Mallet , M. Moaied , J.?J. Palacios , C. Salgado , M.?M. Ugeda , J.-Y. Veuillen , F. Yndurain , and I. Brihuega , Science 352 , 437 ( 2016 ). SCIEAS 0036-8075 10.1126/science.aad8038 [48] 48 G.?Z. Magda , X. Jin , I. Hagymasi , P. Vancso , Z. Osvath , P. Nemes-Incze , C. Hwang , L.?P. Biro , and L. Tapaszto , Nature (London) 514 , 608 ( 2014 ). NATUAS 0028-0836 10.1038/nature13831 [49] 49 J.?J. Palacios , J. Fernández-Rossier , and L. Brey , Phys. Rev. B 77 , 195428 ( 2008 ). PRBMDO 1098-0121 10.1103/PhysRevB.77.195428 [50] 50 A.?F. Young , C.?R. Dean , L. Wang , H. Ren , P. Cadden-Zimansky , K. Watanabe , T. Taniguchi , J. Hone , K.?L. Shepard , and P. Kim , Nat. Phys. 8 , 550 ( 2012 ). NPAHAX 1745-2473 10.1038/nphys2307 [51] 51 W. Bao , L. Jing , J. Velasco, Jr , Y. Lee , G. Liu , D. Tran , B. Standley , M. Aykol , S. Cronin , D. Smirnov , Nat. Phys. 7 , 948 ( 2011 ). NPAHAX 1745-2473 10.1038/nphys2103 [52] 52 Y. Lee , D. Tran , K. Myhro , J. Velasco , N. Gillgren , C.?N. Lau , Y. Barlas , J.?M. Poumirol , D. Smirnov , and F. Guinea , Nat. Commun. 5 , 5656 ( 2014 ). NCAOBW 2041-1723 10.1038/ncomms6656 [53] 53 J. Velasco , L. Jing , W. Bao , Y. Lee , P. Kratz , V. Aji , M. Bockrath , C.?N. Lau , C. Varma , R. Stillwell , D. Smirnov , F. Zhang , J. Jung , and A.?H. MacDonald , Nat. Nanotechnol. 7 , 156 ( 2012 ). NNAABX 1748-3387 10.1038/nnano.2011.251 [54] 54 M. Kharitonov , Phys. Rev. B 86 , 195435 ( 2012 ). PRBMDO 1098-0121 10.1103/PhysRevB.86.195435 [55] 55 E.?V. Castro , N.?M.?R. Peres , T. Stauber , and N.?A.?P. Silva , Phys. Rev. Lett. 100 , 186803 ( 2008 ). PRLTAO 0031-9007 10.1103/PhysRevLett.100.186803 [56] 56 R. Nandkishore and L. Levitov , Phys. Rev. Lett. 104 , 156803 ( 2010 ). PRLTAO 0031-9007 10.1103/PhysRevLett.104.156803 [57] 57 O. Vafek and K. Yang , Phys. Rev. B 81 , 041401 ( 2010 ). PRBMDO 1098-0121 10.1103/PhysRevB.81.041401 [58] 58 A.?S. Mayorov , D.?C. Elias , M. Mucha-Kruczynski , R.?V. Gorbachev , T. Tudorovskiy , A. Zhukov , S.?V. Morozov , M.?I. Katsnelson , V.?I. Fal’ko , A.?K. Geim , and K.?S. Novoselov , Science 333 , 860 ( 2011 ). SCIEAS 0036-8075 10.1126/science.1208683 [59] 59 Y. Lemonik , I. Aleiner , and V.?I. Fal’ko , Phys. Rev. B 85 , 245451 ( 2012 ). PRBMDO 1098-0121 10.1103/PhysRevB.85.245451 [60] 60 R.?E. Throckmorton and S. Das Sarma , Phys. Rev. B 90 , 205407 ( 2014 ). PRBMDO 1098-0121 10.1103/PhysRevB.90.205407 [61] 61 A.?L. Rakhmanov , A.?V. Rozhkov , A.?O. Sboychakov , and F. Nori , Phys. Rev. Lett. 109 , 206801 ( 2012 ). PRLTAO 0031-9007 10.1103/PhysRevLett.109.206801 [62] 62 M. Schüler , M. Rösner , T.?O. Wehling , A.?I. Lichtenstein , and M.?I. Katsnelson , Phys. Rev. Lett. 111 , 036601 ( 2013 ). PRLTAO 0031-9007 10.1103/PhysRevLett.111.036601 [63] 63 A. Auerbach , Interacting Electrons and Quantum Magnetism ( Springer Science & Business Media , New York, 2012 ). [64] 64 B. Bernu , C. Lhuillier , and L. Pierre , Phys. Rev. Lett. 69 , 2590 ( 1992 ). PRLTAO 0031-9007 10.1103/PhysRevLett.69.2590 [65] 65 L. Capriotti , A.?E. Trumper , and S. Sorella , Phys. Rev. Lett. 82 , 3899 ( 1999 ). PRLTAO 0031-9007 10.1103/PhysRevLett.82.3899 [66] 66 P.?W. Anderson , Mater. Res. Bull. 8 , 153 ( 1973 ). MRBUAC 0025-5408 10.1016/0025-5408(73)90167-0 [67] 67 Z. Meng , T. Lang , S. Wessel , F. Assaad , and A. Muramatsu , Nature (London) 464 , 847 ( 2010 ). NATUAS 0028-0836 10.1038/nature08942 [68] 68 Z.?A. Kelly , M.?J. Gallagher , and T.?M. McQueen , Phys. Rev. X 6 , 041007 ( 2016 ). PRXHAE 2160-3308 10.1103/PhysRevX.6.041007, and Marie-Curie-ITN 607904-SPINOGRAPH Ministerio de Economía y Competitividad http://dx.doi.org/10.13039/501100003329 Ministry of Economy and Competitiveness MINECO http://sws.geonames.org/2510769/ FIS2015-65706-P Ramon y Cajal RYC-2013-14645. We acknowledge financial support from the Marie-Curie-Initial Training Networks (ITN) program through Grant No. 607904-SPINOGRAPH, and the Spanish Ministry of Economy and Competitiveness through Grants No. FIS2015-65706-P (MINECO/FEDER) and No. RYC-2013-14645 (Ramon y Cajal program). L.?A.?G.-A. is grateful for the hospitality of the Applied Physics Department in the University of Alicante and to N. Garcia for useful discussions. We specially thank J. Fernandez Rossier for his help settling the environment and the initial idea for this work.