In situ controllable magnetic phases in doped twisted bilayer transition-metal dichalcogenides

We study the electronic structure of hole-doped transition metal dichalcogenides for small twist-angels, where the onsite repulsion is extremely strong. Using unbiased diagrammatic Monte Carlo simulations, we find evidence for magnetic correlations which are driven by delocalization and can be controlled in situ via the dielectric environment. For weak spin-orbit coupling, we find that the moderately doped system becomes anti-ferromagnetic, whilst the regime of strong spin-orbit coupling features ferromagnetic correlations. We show that this behavior is accurately predicted by an analytical theory based on moment expansion of the Hamiltonian, and analysis of corresponding particle trajectories.

Moiré materials-resulting from the incommensurate stacking of single atomic layers-have emerged as an important platform for exploring strongly correlated physics. This field was initiated by the discovery of superconductivity [1] and strongly correlated phases [2] in twisted bilayer graphene (TBG) at "magic angles" where virtually flat bands occur in the spectrum [3]. While TBG remains a very active research field, it presents certain problems due to its theoretical complexity and the fragile nature of the electronic structure. TBG defies the construction of symmetric Wannier functions, meaning that it cannot be accurately described by a conventional lattice model [4], and this complicates theoretical analysis. Experiments are faced with the problem that the interesting correlated phases of TBG only exist for certain twistangles, thus requiring extensive fine tuning.
The twisted bilayer transition metal dichalcogenides (tTMD) have been identified as an attractive alternative to TBG [5] due being robust, highly tunable via the twist-angle and the dielectric environment, and relatively simple to model theoretically [6]. To the first order approximation they can be described by a spin-orbit coupled Hubbard model on the triangular lattice. The bandwidth of the tTMD's can be tuned over several orders of magnitude by the twist angle, while spinorbit coupling and doping can be controlled in situ [7]. The electronic properties also evolve continuously over an extensive range of twist-angels, allowing for greater experimental control.
Experiments on the tTMD's have revealed a rich phenomenology, and the combination of a high degree of control and a simple theoretical description suggests that a more comprehensive understanding of strongly correlated physics may now be possible. In particular, several enigmatic phenomena observed in the cuprate superconductors have been reproduced recently. At a 4.2 • twist-angle, a strange metal regime has been identified at small doping, giving way to a fermi liquid at higher carrier concentrations [8]. At a 5.1 • twist, superconductivity has been observed with a critical temperature of T c ≈ 3K [9] . These twist angles correspond to a parameter regime where the onsite repulsion is comparable to the bandwidth [7,10], implying strong magnetic correlations. For physically realizable model parameters, super-exchange processes are predicted to favor anti-ferromagnetism [7,11], and so this scenario bears a striking similarity to the cuprate superconductors [12]. Recently, strange metal behavior was also predicted theoretically using dynamic mean-field theory [10].
Thus far, both theory and experiments have focused on the scenario of a large twist angle where the bandwidth is comparable to the onsite repulsion. At half-filling or small doping, this results in super-exchange processes with an energy scale that is comparable to the hopping integral. At the smallest twist-angles, this situation changes dramatically as the onsite repulsion becomes much larger than the bandwidth, effectively suppressing super-exchange.
In this paper we focus on magnetism in the regime of small twist angles. Using diagrammatic Monte Carlo [13], we establish that delocalization of the charge carriers can induce either ferromagnetism or anti-ferromagnetism, depending on the spin-orbit coupling, which is in turn controlled by the dielectric environment.
Model-The low-energy physics of the tTMD's can approximately be described by a generalized "Moiré Hubbard model" on the triangular lattice with spin-orbit coupling [6,7]: Here, σ refers to both spin and valley, as these are locked in transition metal dichalcogenides [14,15]. The phase φ is determined by the potential difference between the two layers, and is thus controllable in situ. DFT calculations indicate that realistically achievable parameters correspond to approximately 0 ≤ φ ≤ π/3 [9]. In this model, interchanging the spin corresponds to taking φ → −φ, while particle-hole transformation is equivalent to taking φ → φ + π. The hopping integral |t| decays exponentially with the moiré period, which is in turn determined by the twist angle θ. As a result, the bandwidth can be tuned over several orders of magnitude in the interval 1 • ≤ θ ≤ 5 • [7]. By contrast, the onsite repulsion changes by less that a factor 2 in the same interval, thus permitting the relative interaction U/t to change dramatically. At θ = 4 − 5 • , the hopping integral has been estimated to |t| ∼ 100K while U/t ∼ 8 [10], and so the regime θ ≈ 1 • corresponds to extremely large onsite interactions with a corresponding suppression of super-exchange processes. Higher order corrections to the model (1) consist arXiv:2204.03917v1 [cond-mat.str-el] 8 Apr 2022 of longer range hopping and also non-local repulsive interactions.
In this work, we focus on the physics corresponding to small twist-angles, where U is much greater than the bandwidth. For our calculations, we will specifically consider the strong coupling limit U → ∞, though our results remain valid for a range of finite couplings as well. In this regime, magnetic correlations emerge as the effective density of states of a charge carrier is renormalized by the mean-free path [16]. On bipartite lattices, the lowest energy for a single charge carrier is obtained on a polarized background, leading to the well known Nagaoka theorem [17], with generalizations to finite doping as well [18].
For the Moiré Hubbard model, no analytical results exist in Nagaoka's scenario, though certain qualitative assessments can be made based on the moments of the kinetic energy. First, we note that the moments of the dispersion k can be connected to an expansion of H for a state consisting of a single localized hole on a polarized background (denoted by ψ), Here, the expansion in H can be regarded as a quantum walk, with non-vanishing contributions resulting when the system returns to the initial state. For a polarized background, the mean-free path is infinite, and the dopant is correspondingly an ideal fermion. By contrast, other types of backgrounds lead to a reduced mean-free path that renormalize the effective density of states seen by the carrier [16]. There are two principal types of trajectories in the quantum walks (2): Self-retracing paths that inherently preserve the spin-background, and non-retracing paths that exchange elements of the background, see also Fig. 1. It is only the latter of the two that depends on the spin-background. Since all odd terms in the expansion (2) must enclose a non-vanishing area, it follows that they are non-retracing. Hence, all odd terms in the expansion (2) are significantly reduced if the mean-free path is finite, as occurs in the absence of ferromagnetic correlations. Furthermore, most of the effective bandwidth of a carrier on a Mott insulating background is contained in the self-retracing paths [16,19].
From the preceding considerations, we arrive at the following conjecture for the interplay of delocalization and magnetic correlations: When the system is weakly to moderately holedoped, the charge carriers will occupy the top of the effective band and provide a kinetic energy that is dE/N ≈ − max , where N is the number of carriers and max is the effective band top. If the odd terms in (2) increase the band top max , then delocalization will drive drive the system towards ferromagnetism, whilst if they decrease the band top, then antiferromagnetic correlations will be preferred, as this minimizes the mean-free path. Fig. 2 shows how the band top and bottom depend on the spin orbit coupling. At the edge of the parameter space (φ = ±π/3) the band top is situated at 6t, while the bottom corresponds to −3t, implying that the odd terms in (2) increase max . At φ = 0, this situation is reversed. Based on the preceding conjecture, this implies that the system becomes ferromagnetic at the edges of the parameter space, and antiferromagnetic at the origin, with a crossover occurring at φ c ≈ ±π/6.
Because of spin-orbit coupling, SU(2) symmetry is explicitly broken down to Z 2 -corresponding to time-reversal-when φ = nπ. This means that the Mermin-Wagner theorem [20] does not apply, and that a long-range ferromagnetic order may persist to finite temperatures, despite that the system is 2D. for |φ| ≤ π/3 and σ = 1. At the edge of the parameter region (φ = ±π/3), the band top is situated at k = 6t while the bottom corresponds to −3t. This asymmetry results from odd terms in (2), and allows a hole to delocalize with an energy of −6t on a ferromagnetic background, providing a strong case for ferromagnetism. At φ = 0, the situation is the precise opposite.
Numerical treatment-To the test the preceding conjecture on how delocalization interacts with magnetism, we employ diagrammatic Monte Carlo, which is a numerical protocol based on stochastic sampling of Feynman type graphs [13]. This method has the advantage that it can be employed directly in the macroscopic limit, and is asymptotically exact, as the only systematic source of error is truncation of the series. To be able to address systems with strong interactions we use a particular formulation known as strong-coupling diagrammatic Monte Carlo (SCDMC) [21][22][23][24], where the diagrammatic elements are connected vertices of propagating electrons that are non-perturbative in U . The computational protocol employed here is outlined in detail in [21].
The expansion parameter is the hopping integral t. The principal observable that we compute is the polarization operator of the hopping integral, which we denote Π t (ω, k). The dressed hopping integral is then obtained via the Bethe-Salpiter equation according tõ By conducting a skeleton expansion in the dressed hopping integralt, and iterating until convergence, we obtain a selfconsistent solution that implicitly takes into account certain classes of diagrams up to infinite order. The dressed Greens function can be derived from the polarization according to Because SU(2) symmetry is explicitly broken, we expect the onset of ferromagnetism to be accompanied by critical behavior, which implies a divergent susceptibility with respect to symmetry-breaking perturbations. Furthermore, even in the proximity of a phase transition, the correction from diagrams at higher order tend to be large. This presents a challenge for diagrammatic methods, which are based on the expansion of an analytical function. To overcome this difficulty, we impose a strong symmetry-breaking perturbation on the system, in the form of an external field where t is the bare hopping integral, and β is the inverse temperature. At half-filling or in the atomic limit, the response to the perturbation (5) is that of a paramagnet, so that which defines the relative density of the minority component, ρ min . We choose B 0 ≈ 2.946, so that ρ min = 0.05 for a paramagnet. Next, we consider hole-doping the system. We choose our chemical potential such that the carrier density is 20% when the system is fully polarized, though the resulting equation of state will depend on renormalization of the kinetic energy due to a finite mean-free path in the system. Now, we conduct a self-consistent expansion in t to obtain the Greens function. In Fig. 3 we display the relative density of the minority component for values of the spin-orbit coupling in the interval 0 ≤ φ ≤ π/3 and an inverse temperature βt = 3. In the atomic limit we expect ρ min = 0.05, and so a deviation from this figure is a result of delocalization. For small values of φ, we find that the relative density of the minority component increases, implying that a short mean-free path decreases the kinetic energy of charge carriers, thus driving the system to form anti-ferromagnetic correlations. At strong spin-orbit coupling, we instead see a reduction of the minority component density, implying a preference for a long mean-free path and thus ferromagnetism. The crossover occurs at φ c ≈ 0.5 which is very close to our prediction that φ c = π/6 ≈ 0.52.
Next, we focus on the point φ = π/3, where we expect, and also find, the largest reduction of the minority component density. In Fig. 3 (g) we observe that ρ min is oscillating with an amplitude that falls off rather slowly, suggesting that we are close to the convergence radius of the expansion. To reach lower temperatures we therefore employ a resummation scheme of the form Here, W i ({x}) denotes the weight of the topology i with internal variables {x}, while N i is the expansion order of this topology, and ξ is the resummation parameter. From the reweighted series we obtain observables which formally depend on the expansion order and the resummation parameter ξ, so that For a series with a finite convergence radius we expect that ρ min (N, ξ) converges as N → ∞ for a range ξ ≥ ξ min . Our aim is now to estimate ρ min (N → ∞, ξ ≥ ξ min ) and extrapolate this result to ξ = 0. Thus, we compute ρ min (N, ξ) for N = 8, 9, 10 and a range of ξ. Fig. 4 (a) shows ρ min as a function of the resummation parameter for inverse temperatures 3 ≤ β ≤ 3.75. For ξ ≥ ξ min ∼ 0.45 − 0.5 the lines corresponding to different expansion orders N coincide, indicating convergence. Thus, we fit a second order polynomial to the simulation data in the region ξ ≥ ξ min . Extrapolating this function to ξ = 0 gives an estimate for ρ min (N → ∞, ξ = 0). A summary of the estimate is given in Fig. 4 (b): As the temperature is reduced, ρ min decreases monotonically. It should be stressed, that this calculation is conducted for a magnetic field that decreases with temperature according to Eq. (5), so that the reduction of the minority density is entirely the result of delocalization. The error bars in Fig. 4 (b) were obtained by varying the cutoff ξ min and observing how the result changes. However, it should be stressed that there is also an implicit error associated with the choice of fitting function. We also tried other fitting functions, including a third order polynomial, though this resulted in greater sensitivity to noise in the simulation data. Nonetheless, the key qualitative feature of the seriesthat ρ min decreases as the temperature is reduced-is reproduced for all fitting functions we tried.
Discussion-We have identified a mechanism by which delocalization of charge carriers gives rise to magnetic correlations in the moiré Hubbard model. Based on analysis of particle paths, we conjecture a cross over from anti-ferromagnetic correlations at weak spin-orbit coupling to a ferromagnetic regime at strong coupling. Comparison of this theory to diagrammatic Monte Carlo simulations show excellent agree- . Relative density of the minority component (ρmin = n ↑ /(n ↑ + n ↓ )) obtained from diagrammatic Monte Carlo expansion up to an order N = 8 or N = 9 (a-g). The expansion is conducted for spin-orbit couplings in the interval 0 ≤ φ ≤ π/3 and in a strong external field according to (Eq. 5). The external field is chosen such that ρmin = 0.05 in the atomic limit (see Eq. 6 ), and the deviation from this figure is thus a result of delocalization of charge carriers. The system is hole-doped, and the chemical potential in (a-g) is chosen such that the carrier density is 20% if the system is fully polarized. The shaded region is an estimate for ρmin as N → ∞. The results in (a-g) are summarized and compared to the paramagnetic response in (h). At small spin-orbit coupling, we observe an increase in the minority component, implying that delocalization favors a short mean-free path, and thus anti-ferromagnetism. At large spin-orbit coupling, the preference is for ferromagnetism. The critical point separating these regions can be identified by the crossing of the Monte Carlo data (red dashed line) and the paramagnetic response (green dotted line), which occurs at φc ≈ 0.5. We may compare this result to the initial conjecture: We expect that the type of magnetism favored by delocalization is primarily determined by the magnitude difference between the band top ( max) and bottom ( min). Correspondingly, we expect the critical point to occur at φc = π/6 ≈ 0.52. The predicted phase diagram shown in (i) is thus in excellent agreement with our findings. ment. As the temperature is reduced, we confirm a rapid reduction of the minority component in the strong-coupling regime, indicating that the system eventually becomes fully polarized.
While the results reported here were obtained in the strongcoupling limit, they remain relevant for physical parameter regimes of the TTMD's. For example, at twist-angles of θ = 4 − 5 • , estimates of effective model parameters yield |t| ∼ 100K and U/t ∼ 8 [10]. If the hopping integral is shrunk by two orders of magnitude via the twist-angle, and the onsite repulsion changes by less than a factor 2 [7], then this gives a super-exchange J/t ∼ 4t/U ∼ 100 and t ∼ 1K. The energy scale of a magnetic bond is thus M ∼ t/400, implying that M β ∼ 10 −2 in our simulations. Thus, a conservative estimate in this scenario is that at the energy scale of magnetic bonds would be at least two orders of magnitude smaller than the simulated temperature, which corresponds to 0.25 − 0.33K. This clearly demonstrates the existence of parameter regimes where super-exchange is irrelevant, and the strong coupling limit is an accurate description of the system.  . Temperature dependence of the relative density of the minority component (ρmin = n ↑ /(n ↑ + n ↓ )), obtained by resummation of the series. By reweighting diagram topologies according to Eq. (7), we obtain a series for ρmin(N, ξ). We computed ρmin(N, ξ) for N = 8, 9, 10 and a range of resummation parameters ξ (a). For ξ ≥ ξmin ∼ 0.45 − 0.5, the solutions for different expansion orders N coincide, indicating that the series for ρmin(N, ξ) has converged to ρmin(N → ∞, ξ). Thus, we fit a second order polynomial to the simulation data in the region ξ ≥ ξmin where ρ has converged with respect to N , indicated by the solid orange line.
The polynomial provides an extrapolation to ξ = 0, indicated by the dashed orange line. We extracted these results for inverse temperatures 3 ≤ βt ≤ 3.75. The estimate for ρ(N → ∞, ξ = 0) is summarized in (b). As the temperature decreases, the relative density of the minority component drops. It should be noted that the applied field is taken to be B = B0t/β, and so the decrease of the minority component density is only driven by the increasing energy scale of delocalization as compared to temperature. The system considered is hole-doped at ∼ 20%.
In the polarized regime, the moiré Hubbard model (1) reduces to a system of ideal single-component fermions. However, if we also include non-contact interactions, then we obtain an effective "polarized moiré model" of the form Here, we have not included beyond nearest-neighbour hopping as these terms are extremely small, and most likely not important. The nearest neighbour interactions can be quite large however: An approximative assessment suggests that U/V N N ∼ 15 at a twist angle of 1 • and U/V N N ∼ 3 at 5 • [7]. In our example, where U/t ∼ 400, this would imply that V N N > 25t, so that the interaction is much larger than the bandwidth. Even beyond nearest neighbour interaction, the energy scale can be comparable to or larger than the bandwidth. We expect nonlocal interactions to have limited impact on kinetically driven magnetism, since these terms do not affect the interaction between charge carriers and the spinbackground. However, at much lower temperatures they may drive instabilities in the system.
Finally, we note that the magnetic phases that appear for small twist-angles can be controlled via the dielectric environment. In the ferromagnetic regime, we expect that the minority component becomes gapped at sufficiently low temperature, implying that the system only transmits currents of the majority component. This would allow TTMD's to be used to construct "spin-filters" that are controllable in situ.
This work was supported by the Swedish Research Council (VR) through grant 2018-03882. Computations were performed on resources provided by the Swedish National Infrastructure for Computing (SNIC) at the National Supercomputer Centre in Linköping, Sweden. The author would like to thank Ahmed Abouelkomsan and Emil Blomquist for important input and discussions.