Kinetic magnetism in triangular moiré materials

Magnetic properties of materials ranging from conventional ferromagnetic metals to strongly correlated materials such as cuprates originate from Coulomb exchange interactions. The existence of alternate mechanisms for magnetism that could naturally facilitate electrical control has been discussed theoretically1–7, but an experimental demonstration8 in an extended system has been missing. Here we investigate MoSe2/WS2 van der Waals heterostructures in the vicinity of Mott insulator states of electrons forming a frustrated triangular lattice and observe direct evidence of magnetic correlations originating from a kinetic mechanism. By directly measuring electronic magnetization through the strength of the polarization-selective attractive polaron resonance9,10, we find that when the Mott state is electron-doped, the system exhibits ferromagnetic correlations in agreement with the Nagaoka mechanism.

Magnetic properties of materials ranging from conventional ferromagnetic metals to strongly correlated materials such as cuprates originate from Coulomb exchange interactions.The existence of alternate mechanisms for magnetism that could naturally facilitate electrical control have been discussed theoretically [1][2][3][4][5][6] but an experimental demonstration [7] in an extended system has been missing.Here, we investigate MoSe 2 /WS 2 van der Waals heterostructures in the vicinity of Mott insulator states of electrons forming a frustrated triangular lattice and observe direct evidence for magnetic correlations originating from a kinetic mechanism.By directly measuring electronic magnetization through the strength of the polarization-selective attractive polaron resonance [8,9], we find that when the Mott state is electron doped the system exhibits ferromagnetic correlations in agreement with the Nagaoka mechanism.
Moiré heterostructures of two-dimensional materials provide a new paradigm for the investigation of the physics of strongly correlated electrons.In contrast to well-studied quantum materials, they provide a very high degree of tunability of the parameters relevant for controlling correlations, such as the carrier density and the ratio of interaction energy to hopping strength.Moreover, unlike cold-atom quantum simulators, physics and functionality of moiré materials can be varied using readily accessible external electric and magnetic fields, creating a platform where different many-body phases compete.During the five years since the first realization of a moiré material, a wealth of correlation physics ranging from correlated Mott-Wigner states through quantum anomalous Hall effect to superconductivity has been observed both in magic angle twisted bilayer graphene (MATBG) and in bilayers of transition metal dichalcogenides (TMDs) [10][11][12][13][14][15][16][17][18][19][20][21][22][23].With the exception of orbital magnetism in MATBG [14,24,25] as well as early spin susceptibility and scanning probe measurements in TMD bilayers [26][27][28][29] however, quantum magnetism in moiré materials has until recently remained experimentally unexplored.Theoretical works on the other hand have investigated the magnetic properties of the correlated Mott state in a moiré lattice with one electron per lattice site [30,31] and focused on the possibility to realize quantum spin liquids [31][32][33][34].
Here we investigate the magnetic properties of electrons in MoSe 2 /WS 2 heterobilayers using lowtemperature confocal microscopy.We focus on the magnetization as a function of temperature T and out-ofplane magnetic field B z at dopings around one electron per moiré lattice site (ν = 1).For ν > 1, our experiments show that the itinerant electrons exhibit a positive Curie-Weiss constant θ CW .The abrupt jump in spin susceptibility as soon as the electron density is tuned beyond ν = 1 at T ≈ 170 mK is consistent with kinetic ferromagnetic correlations linked to the Nagaoka mechanism [1].
We study two R-type MoSe 2 /WS 2 heterostructures encapsulated in hBN.The lattice mismatch and twist angle between the TMD monolayers creates a moiré superlattice with a lattice constant of about 7.5 nm.The minima of the resulting electronic potential for the conduction band are located at the high symmetry points where the metal atoms in the two layers are aligned (MM sites).Injected electrons occupy the triangular lattice of MM sites as illustrated in Fig. 1a.In Sample I, the charge density and the electric field in the heterostructure can be tuned independently using top and bottom graphene gates, while Sample II is only single-gated.
The reflection spectrum as a function of electron density in Fig. 1b shows multiple resonances close to the energy of the optical transitions in monolayer MoSe 2 .Intensity maxima and cusps in the resonance energies appear at equally spaced gate voltages (see Extended Data Fig. 1 for extended range).These voltages correspond to commensurate filling of the moiré superlattice with one or two electrons per site (ν = 1 and 2, respectively), where incompressible states are formed.We focus here on the resonance at 1.58 eV, which can be identified as an attractive polaron (AP) resonance associated with collective excitation of bound electron-exciton pairs (tri- ons) located at the moiré potential minima [8,9,35].As shown in Fig. 1c, the area, or equivalently the oscillator strength, of the AP resonance increases linearly as function of electron density up to ν = 1 and subsequently decreases again linearly between ν = 1 and 2. This behaviour suggests the presence of an isolated Hubbard band where all electrons occupy the same lattice sites, forming double occupancies (doublons) for ν > 1.Since the AP resonance is associated with the bound trion of an exciton and a resident electron, it can only be optically excited on moiré lattice sites already occupied by a single electron.Consequently, the densities ν = ε and ν = 2 − ε provide the same number of sites for AP formation and hence lead to an identical oscillator strength of the AP resonance.
The optical selection rules of monolayer MoSe 2 are retained in the heterostructure, giving rise to circularly polarized resonances for B z = 0, corresponding to transitions in the K and K' points of the MoSe 2 Brillouin zone [36].The linear dependence of the AP peak area on the electron density, together with the optical valley selection rules and strong spin-orbit coupling leading to spin-valley locking, allows us to employ the polarization resolved AP resonance as a quantitative probe of the degree of spin polarization of the electrons.Since the AP is only formed by excitons in the K valley and spin-down electrons in the K' valley or vice versa [35], the AP oscillator strength in σ + -polarization (σ − -polarization) is proportional to the density n ↓ (n ↑ ) of spin-down (spin-up) electrons.The degree of spin polarization is then given by where A σ ± is the area of the AP resonance in σ ± polarization and ρ AP denotes the degree of circular polarization of the AP resonance.The polarization resolved spectrum in Fig. 1d, measured at B z = 1 T, T = 4.2 K, and ν = 1, highlights how the AP resonance becomes partially polarized in a moderate magnetic field.Note that the resonance at 1.635 eV is also sensitive to the spin polarization, mainly through a splitting with giant effective g-factor g eff = 31, as previously reported for other moiré heterostructures [26,27].
In order to gain insight into the interactions between spins of the electrons residing in the superlattice potential, we measure the AP degree of polarization ρ AP as a function of B z for electron densities satisfying 0.5 < ν < 1.8.An excitation power of 4.7 pW is used to avoid laserinduced spin depolarization (see Methods) and thereby ensure that we probe magnetic properties of the electronic ground state [37].We perform a linear fit to extract the slope at B z = 0, as shown in Fig. 2a, which is related to the magnetic susceptibility through where M (ν) is the magnetization, µ 0 the vacuum permeability, and M s (ν) = gµ B n e,ν (1 − |ν − 1|) /2 the saturation magnetization for each density, with µ B the Bohr Strikingly, there is a similarly sharp decrease again at ν = 1.5, even though this is not a characteristic filling of the triangular lattice and an ordered state underlying the abrupt change can only emerge with additional symmetry breaking.We focus here on the density dependence of magnetic interactions around ν = 1 and further investigate them by performing temperature-dependent measurements shown in Fig. 2c.Each curve measured at a temperature T is normalized by 1/T such that for paramagnetic behaviour the curves collapse onto one value, as observed for ν < 1.The deviation from this universal behaviour for ν > 1 is evidence for the presence of ferromagnetic interactions.
In general, exchange interactions are expected to play a key role in determining the magnetic order of the system.For a moiré structure with lattice constant of 7.5 nm in particular, first principle calculations predict reduced superexchange due to the large on-site repulsion, such that direct exchange is the relevant interaction to consider [31].The magnetic properties at ν = 1, where the electrons form an incompressible Mott insulator and are localized on moiré lattice sites, should be exclusively determined by exchange interactions.Surprisingly, we do not find a significant deviation from paramagnetic behaviour at ν = 1, as can be seen in Fig. 2d from the fit- ted Curie-Weiss constant θ CW as function of doping.This suggests that exchange interactions do not play a significant role in our system.Possible explanations for this are a cancellation of competing direct and superexchange interactions or that the electrons are more strongly localized than predicted by theory due to a deeper moiré potential.
In addition to exchange, effective magnetic interactions based on a kinetic mechanism appear when the electrons are mobile [38]: minimization of the kinetic energy gives rise to magnetic order in a Hubbard band close to half filling (ν = 1) even in the limit of vanishing spin interactions [1].We attribute the observed ferromagnetic correlations with a sharp onset at ν = 1 to kinetic magnetism.This is corroborated by the dip in susceptibility at ν = 4/3, marked by the dashed line in Fig. 2b: At fractional filling factors commensurate with the moiré superlattice, electrons have been shown to form incompressible Mott-Wigner states [18,19,21,22,28], reducing their mobility and suppressing kinetic magnetism.
In a simple single-band Fermi-Hubbard model in the strongly interacting regime, the kinetic mechanism leads to a transition from antiferromagnetic to ferromagnetic interactions close to ν = 1 on a triangular lattice [4][5][6].To qualitatively understand the experimentally observed absence of antiferromagnetic correlations for ν < 1, we consider an extended model (see Methods) taking into account long-range Coulomb interactions V , while setting exchange interactions J to zero, motivated by the paramagnetic response at ν = 1.The Coulomb interaction term modifies the hopping of electrons onto sites that are already occupied, which renormalizes the doublon hopping, while leaving the hole hopping t unchanged, as illustrated in Fig. 3a.The effective doublon hopping is given by t+A, with the assisted hopping term A = − w i , w i | V |w i , w j , where w i and w j denote states localized on neighbouring sites [30,31].Due to the asymmetry in hopping between holes and doublons, the kinetic magnetism is enhanced for ν > 1.Therefore, in an intermediate temperature range t k B T t + A, we expect a sizeable modification of the susceptibility for ν > 1, but only negligible deviations from a paramagnetic response for ν ≤ 1.This asymmetry of the susceptibility around ν = 1 is captured by our finite temperature tensor network simulations shown in Fig. 3b.Details on the theoretical model employed and the simulations can be found in the Methods section.
Furthermore, we find that the overall hopping strength is renormalized by the presence of long-range interactions and/or disorder, reducing the strength of kinetic magnetism, particularly for ν < 1.In Fig. 3c we show simulated magnetization curves at ν = 0.89 and T = 0 comparing the cases with and without long-range interactions or disorder.Introducing disorder or interactions leads to an increased slope at low fields, corresponding to an enhanced susceptibility or suppressed antiferromagnetic correlations.This limits the temperatures required to observe kinetic magnetism in the moiré structure to smaller values than what would be expected from theoretically predicted hopping strengths on the order 1 meV for such moiré structures.
Conclusion-We have demonstrated accurate determination of the spin polarization through polarized AP oscillator strength measurements at picowatt power levels.We have used this method to study the densitydependent spin susceptibility and found a sudden appearance of ferromagnetic correlations for ν > 1. Supported by tensor network simulations, our observations can be attributed to kinetic magnetism in an extended single-band Fermi-Hubbard model on a triangular lattice.Even though prior studies found good agreement with ab initio calculations, our observation of a paramgentic response at ν = 1 suggests that direct exchange and superexchange interactions are weak and the spin state is dominated by effective kinetic interactions.The strong asymmetry between ν < 1 and ν > 1 suggests the presence of a large Coulomb-assisted hopping of doublons.Moreover, Coulomb interaction and disorder renormalize the effective hopping for holes and doublons, which in turn reduces the scale of magnetic correlations to low temperatures below 100 mK.
Note added -During the preparation of this manuscript, we became aware of several parallel works exploring different aspects of magnetism in bilayer MoTe 2 moiré structures [39][40][41].
Extended Data Figure 1.Normalized reflectance spectrum as function of doping, Sample I, extended range.

Sample Fabrication
Graphene, hBN, MoSe 2 , and WS 2 were exfoliated from bulk crystals onto Si/SiO 2 (285 nm) substrates and assembled in heterostructures using a standard dry-transfer technique with a poly(bisphenol A carbonate) film on a polydimethylsiloxane (PDMS) stamp.Both samples were encapsulated between two approximately 30 nm thick hBN flakes.Optical lithography and electron beam metal deposition were used to fabricate electrodes for the electrical contacts.

Experiment Setup
Experiments were done in a liquid helium bath cryostat and a dilution refrigerator with free-space optical access, with windows at the still, 4 K, and room temperature stages.The samples were mounted on three piezoelectric nanopositioners.A fiber-coupled confocal microscope was used for optical measurements, with a either a single aspheric lens or an objective (NA = 0.7 for both) focusing the light to a diffraction-limited spot.A schematic of the optical setup is shown in Extended Data Fig. 2. Reflection spectra were measured using a supercontinuum laser with a variable filter as light source and a spectrometer with a liquid-nitrogen-or Peltier-cooled CCD camera as detector.For resonant single-frequency measurements of the magnetic circular dichroism (MCD), we used a tunable CW titanium sapphire laser and a Geiger-mode avalanche photodiode for detection of low-power signals.
Extended Data Figure 2. Schematic of the optical setup.
The laser polarization was switched between σ + and σ − at kilohertz rates using an electro-optico modulator in order to reduce the sensitivity of the MCD measurement to slow drifts.All measurements were power-stabilized with feedback from a photodiode to an acousto-optic modulator or a fiber-coupled variable optical attenuator.
Extended Data Figure 3. High temperature data, Sample II. a Inverse slope dρAP/dB measured at ν = 0.75 as function of temperature for different excitation powers.For each power, the data are precisely reproduced by the Curie-Weiss formula with the Curie-Weiss constant getting sizably lower for larger powers due to light-induced electron spin depolarization.b Fitted Curie-Weiss constant as function of electron filling factor for three different excitation powers.For ν < 1, the light-induced spin depolarization leads to apparent negative Curie-Weiss constants.c Ratio of the slope dρAP/dBz measured at two different powers, plotted as function of temperature and filling factor.The power dependence changes with both temperature and filling factor and is most prominent at small ν and low temperatures.The grey dashed lines mark contours of fixed susceptibility ratios (as indicated).Error bars correspond to the standard error of the fit.

Background Subtraction
Differential reflectance presented in the plots is defined as ∆R/R 0 = (R − R 0 )/R 0 , where R is the measured reflection spectrum of the heterostructure and R 0 is the background reflection spectrum on the hBN flakes away from the TMD flakes.For the MCD measurements, the background reflectance at the laser frequency is measured in both polarizations at charge neutrality or high electron density (ν > 2) where there is no AP resonance.The degree of circular polarization is then given by

Power Dependence of Spin Polarization
To access the true magnetic ground state properties of the system, it is essential to ensure that the intensity of the probe light is sufficiently low to not perturb the system.Similarly to monolayer MoSe 2 [42], light illumination leads to depolarization of the spin population in the moiré heterostructure.Since the strength of the depolarizing effect depends on both temperature and charge density, it can give rise to misleading artefacts in the measured electronic magnetism.This is directly revealed by our filling-factor-dependent measurements of the Curie-Weiss constant carried out at high temperatures T > 4 K in the bath cryostat on Sample II.In these experiments, the sample was illuminated with a white light of tunable power.The magnetic susceptibility of the electron system was extracted based on the degree of circular polarization of the AP resonance that was in turn determined by fitting its spectral profile with a dispersive Lorentzian lineshape [43].On this basis, we were able to analyze the temperature dependence of the inverse magnetic susceptibility for various filling factors and excitation powers.As seen in Extended Data Fig. 3a (for ν = 0.75), even though the employed powers remain in the submicrowatt range, they still markedly affect the magnetic response.More specifically, the Curie-Weiss constant is clearly lower for larger excitation powers.This effect is most prominent for low filling factors and becomes indiscernible at ν 1 (see Extended Data Fig. 3b).This power dependence originates primarily from the changes in spin-valley relaxation dynamics of the electron system.As it has been demonstrated in prior studies of TMD monolayers [44], the spin relaxation time becomes shorter for larger electron densities and higher temperatures.As a result, if a certain number of electrons undergo a spin flip due to the interaction with optically injected excitons, it takes longer for them to relax back to their ground state when ν and T are low.For this reason, the magnetic susceptibility determined upon exciton injection into the system is lower compared to its unperturbed value.Moreover, the deviation between these two quantities becomes larger for higher excitation powers and lower ν and T (see Extended Data Fig. 3c), which explains the striking power dependence of the Curie-Weiss constant at Extended Data Figure 4. Magnetization curves measured with two excitation powers at T ≈ 170 mK in Sample I.The curves measured at powers differing by a factor of 4 overlap, indicating that the exciation power is sufficiently low to not perturb the system.ν < 1 in Extended Data Fig. 3b.In particular, the data in this figure directly reveal that in order for the excitons to constitute a nondestructive probe of the electron spin system at T > 4 K and 0.5 ν 1.5, the excitation power needs to be around a few nW.
Owing to the aforementioned temperature dependence of the spin-valley relaxation time, accessing the true magnetic ground state properties of the electron system at mK temperatures requires us to further reduce the excitation power.We find that the requisite power is on the order of a few pW, as seen in Extended Data Fig. 4, where increasing the light intensity by a factor of four did not affect the outcome of the measurement and the two magnetization curves overlap.Taking this into account, we used a resonant laser with 4.7 pW incident power on the sample in our mK measurements.Note that this level of power is about six orders of magnitude below the level at which the laser measurably heats the cold finger in the cryostat.Because the line shape and energy of the AP are almost constant with respect to gate voltage and magnetic field, measuring the reflectance at a single frequency is equivalent to measuring the area of the peak.

Temperature Calibration
To calibrate the lowest temperatures used in the dilution refrigerator, we first carried out Curie-Weiss fits for T > 300 mK relying on the built-in temperature readout based on a resistance measurement.From these fits, we find Curie-Weiss constants close to 0 for ν ≤ 1, which allows us to assume paramagnetic behaviour dρ AP (T )/dB z = gµ B /(2k B T ), where µ B is the Bohr magneton, k B the Boltzmann constant, and g the electronic gfactor.This is further confirmed by measured magnetization curves that follow ρ AP (B z ) = tanh(gµ B B z /(2k B T )).The value g = 4.5 of the g-factor can be fixed from this relation using the measured magnetization slope and temperature.We then use the same relation to extract the temperature at T < 300 mK using the measured slope at ν = 1.
The slope of dρ AP (T )/dB z at ν = 1 was also utilized to determine the sample temperature in high-temperature measurements that were performed in the bath cryostat.In this case, the obtained temperature values were further verified by analysis of the temperature-induced redshift of the exciton resonance in a MoSe 2 monolayer region of Sample II.As shown in Extended Data Fig. 5, the measured energy E X (T ) of this resonance decreases quadratically with temperature, following aptly the Varshni formula E X (T ) = E 0 − γT 2 [45].The corresponding γ = 1.6 µeV/K 2 agrees very well with the values reported in previous studies of MoSe 2 monolayers carried out in wider temperature ranges [46].This finding provides a strong confirmation of the validity of our temperature calibration procedure.

Theoretical Model
To explain the experimental results we consider a single-band extended t-J model, where t is the hopping strength, J is the spin-spin interaction, A is the assisted hopping, V is the strength of Coulomb interaction projected into the lowest Wannier orbital, h is the external magnetic field in units of gµ B , ∆ i is the on-site potential energy, and P is a projector that projects out doublons (holes) in the hole (electron) doped regime.We consider a null spin-spin interaction J = 0 motivated by the experimental results pointing to a paramagnetic response at ν = 1.Moreover, to implement the long-range coupling proportional to V we cut the range of the interaction at third neighbours.The on-site potential energy ∆ i takes into account the spatial variations of the moiré potential.We consider a uniformly distributed disorder ∆ i ∈ [−∆/2, ∆/2) of width ∆ with a corresponding RMS parameter ∆/ √ 12.

Model Parameters
To estimate the relevant parameters of the Hamiltonian used in the tensor network simulations, we start from the finite discrete Fourier expansion of the moiré potential, where V n = −V 0 exp i(−1) n−1 ϕ and we introduce the reciprocal lattice vectors cos(πn/3) sin(πn/3) .
The parameters V 0 = 6.3 meV and ϕ = 0 are obtained from first principles calculations.The single-electron problem is described by the lowenergy Hamiltonian, where we introduce the effective mass m * = 0.7m e of the MoSe 2 conduction band electrons.Since the moiré potential has a periodic structure we can employ Bloch's theorem to write the wavefunctions as where n is the band index, k is restricted to the first moiré Brillouin zone, and u (n) k are the Bloch functions.Since the Bloch functions have the same periodicity as the moiré potential, u where G is the set of all reciprocal lattice vectors.Therefore, the Hamiltonian can be written in the basis of reciprocal lattice vectors as which can be diagonalized for each quasi-momentum k by using a large set of reciprocal lattice vectors.The ground state solution corresponds to the lowest band n = 0.The associated Wannier wavefunction w i ( r) localized at site R i is obtained by performing the change of basis where we drop the band index and introduce the normalization factor N .The interaction potential between charges within the TMDs is given by the Rytova-Keldysh potential [47,48] V RK (r) = e 2 8 0 r 0 H 0 r r r 0 − Y 0 r r r 0 , Tensor Network Simulations Our finite temperature tensor network simulations are based on a purification scheme performed in the canonical ensemble.We implement the U (1) symmetry associated with the conservation of the total number of electrons, but we do not fix the net magnetization of the system.The finite temperature density matrix is represented as a matrix product state (MPS) in a doubled Hilbert space.The MPS maximum bond dimension is set to χ = 768.The cooling process is performed similarly to Ref. [4].We progressively apply an infinitesimal (δβ = 0.1) Boltzmann factor e −δβ/2 by employing the W II technique [50].The finite temperature calculations are performed in a triangular cylinder of size L = L x × L y = 15 × 3.
To obtain the ground state of the system we employ the density matrix renormalization group (DMRG) algorithm.We perform simulations in a triangular cylinder of size L = L x × L y = 15 × 6 and we fix the maximum bond dimension of our MPS to χ = 1024.To capture the effects of a disordered on-site potential we have performed calculations in three different disorder realizations and average over them.

Figure 1 .
Figure 1. a Moiré potential in the conduction band.Electrons occupy the potential minima at the MM sites.b Normalized reflectance spectrum as function of gate voltage, tuning the doping of the heterostructure Sample II.At integer fillings, intensity maxima and cusps in the resonance energies appear.c Area of the AP resonance as function filling factor at Bz = 0.The linear increase and decrease confirm that electrons occupy a single minimum in the moiré unit cell, forming an isolated Hubbard band.d Polarization resolved reflection spectrum at ν = 1, Bz = 1 T, and T = 4.2 K, Sample II.The AP resonance at 1.58 eV is sensitive to the spin polarization of the electrons through its degree of circular polarization.

Figure 2 .
Figure 2. Low temperature data, Sample I. a Degree of polarization of the AP resonance as function of magentic field and linear fit around Bz = 0 yielding the susceptibility.b Slope of ρAP as function of electron filling factor.A sharp increase in spin susceptibility is observed at ν = 1.c Doping dependence of dρAP/dB at different temperatures, normalized by 1/T .At 1 < ν < 1.5, the susceptibility deviates from 1/T due to ferromagnetic interactions.d Fitted Curie-Weiss constant as function doping.Error bars correspond to the standard error of the fit.

Figure 3 .
Figure 3. a Illustration of hopping processes for holes and doublons.The presence of long-range interactions introduces an assisted hopping A that modifies the doublon hopping strength.b Simulated spin susceptibility as function of filling factor for assisted hopping A = 10t and temperature kBT = A, but absent long-range Coulomb interaction V = 0, disorder ∆ = 0, and exchange interactions J = 0.The assisted hopping leads to an asymmetry between doublons and holes.The normalization is chosen such that the value 1 corresponds to a paramagnetic response.c Simulated degree of spin polarization as function of magnetic field at ν = 0.89 and T = 0.Both disorder with distribution width ∆ and longrange Coulomb interactions V suppress the antiferromagnetic correlations.

where H 0
is the Struve function, Y 0 the Bessel function of the second kind, r 0 = 3.5 nm the screening length for MoSe 2 , and r = 4.5 the relative permittivity of hBN as the surrounding medium[49].The matrix elementst = − w i | Ĥ |w j = 0.75 meV U = w i , w i | V RK (| r 2 − r 1 |) |w i , w i = 157 meV V = w i , w j | V RK (| r 2 − r 1 |) |w i , w j = 44.6 meV (14) J = − w i , w j | V RK (| r 2 − r 1 |) |w j , w i = −0.61meV A = − w i , w i | V RK (| r 2 − r 1 |) |w i , w j = 6.1 meVare evaluated numerically, where |w i , w j denotes a state where two electrons occupy neighbouring Wannier orbitals.