Magnetism and Quantum Melting in Moiré-Material Wigner Crystals

Our starting point is the continuum model description of TMD heterobilayers with type-I or type-II band alignment [1]. In this case the topmost moiré band is concentrated in one of the layers, which we refer to as the active layer, depicted in red in Fig. 1(a)-(b). The influence of the second layer, shown in blue in Fig. 1(a)-(b), is responsible for a moiré potential that affects charge carriers in the active layer. The moiré pattern can be induced by a small twist angle $\theta$ or a lattice mismatch $\delta$ between the two layers. The moiré lattice constant is given by $a_{M}=a_{0}/(\theta^{2}+\delta^{2})^{1/2}$, with $a_{0}$ the lattice constant of the active TMD layer. To date most heterobilayer experiments have focused on unrotated WSe₂/WS₂ with a moiré lattice constant of $a_{M}\approx 8.2$ nm, or MoSe₂/WS₂, with a moiré lattice constant of $a_{M}\approx 7.5$ nm. Because the moiré lattice constant reaches a maximum, the system is expected to be less sensitive to twist angle disorder at zero twist angle.

Valley and spin are locked in TMD heterobilayers and the valley (or spin) projected continuum Hamiltonian is given by [1]

The first term in Eq. (1), which corresponds to the kinetic energy of carriers in the top moiré band, is diagonal in momentum space and the second term, the moiré potential, is diagonal in coordinate space. The moiré potential depends on only two parameters $(V_{m},\psi)$ because [1] of the system’s $C_{3}$ symmetry. The phase $\psi$ fixes the geometry of the moiré superlattice, which we take to be triangular as it is in most of the TMD heterostructures, and the strength of the moiré potential $V_{m}$ can be related to the experimentally tunable displacement field. For concreteness, throughout this work we take the effective mass $m^{*}=0.45\,m_{0}$ and $\psi=45^{\circ}$, corresponding to WSe₂/WS₂ [22]. We take the modulation potential strength $V_{m}$ as an experimentally controllable parameter since it has been demonstrated to be sensitive to the displacement field $D$. The form we have used for the moiré modulation parameter assumes that it varies smoothly with position on the moiré scale; higher harmonics in the plane-wave expansion are more important in longer period moirés in which relation relative to rigidly twisted bilayers is stronger.

An example of the bandstructure obtained from diagonalizing Eq. (1) is shown in Fig. 1(c). We label the band energies of Eq. (1) as $\epsilon_{{\bm{k}}}^{n}$, while the eigenstates can be written in a plane wave expansion as $\ket{n,\psi_{\bm{k}}}=\sum_{{\bm{G}}}z^{n}_{{\bm{k}},{\bm{G}}}\ket{{\bm{k}}+{\bm{G}}}$, where $n$ is a band index and ${\bm{G}}$ are reciprocal lattice vectors. In order to study the emerging many-body phases we consider the interacting Hamiltonian projected to the topmost moiré valence band (hence we omit band index)

where the Coulomb matrix elements are given by

The summation involving the term in brackets results from the projection of interactions to a single band and can be rewritten in terms of the form factors, as described elsewhere [23]. In Eq. (4), $q=|{\bm{k}}_{i}+{\bm{G}}_{i}-{\bm{k}}_{k}-\bm{G}_{k}|$ is the momentum transfer, $A$ is the area of the system, $\varepsilon$ the effective dielectric constant, and total momentum conservation is implicit.

In this article we take the same approach as in previous work [13, 12], utilizing exact diagonalization (ED) to solve the many-body Hamiltonian. We present our results in phase diagrams that depend on dimensionless parameters obtained by taking ratios of the relevant energy scales of the system. In particular, three energy scales can be identified; the moiré potential depth $V_{m}$, the interaction strength $U_{M}=e^{2}/(\varepsilon a_{M})$, and the kinetic energy scale $W_{M}=\hbar^{2}/m^{*}a_{M}^{2}$. By varying the two ratios of these three scales, we can simulate any heterobilayer as long as the moiré period is much larger than the microscopic lattice constant. This ensures that our conclusions apply for arbitrary heterobilayers as long as their low energy physics is captured by the continuum model, Eq. (1) and that the shape of the moiré potential is similar to that of a triangular lattice model. The energy gap to the remote moiré bands $\Delta_{R}$ can, in principle, be viewed as another relevant energy scale. We take parameter values such that this scale is always larger than all other scales involved, as can be seen in Fig. 1(d), justifying the single band projection of Eq. (3). As we have pointed out above, however, interactions renormalize bands more at higher electron densities. For this reason the single-band approximation should be treated with caution for fillings $\nu>1$, where remote band mixing is often relevant.

The model we study has orbital and spin degrees of freedom, discrete triangular lattice translational symmetry, $SU(2)$ spin-rotational invariance, and no spin-orbit coupling. The Hilbert space can be divided into smaller subspaces with total momentum ${\bm{K}}$ with discrete values determined by the number of moiré unit cells $M$, total spin ${\bm{S}}$, and azimuthal spin $S^{z}$. The basis is constructed in an occupation number representation, distributing particles among single-particle states labeled by $S^{z}$ quantum number and quasi-particle $(k_{x},k_{y})$ momentum. The total number of possible configurations $N_{conf}$ for particles distributed on $M$ single particle states with $N_{\uparrow}$ spins up and $N_{\downarrow}$ spins down is determined by a product of binomial coefficients, $N_{conf}=\binom{M}{N_{\uparrow}}\cdot\binom{M}{N_{\uparrow}}$. The many-body Hamiltonian projected to a given total momentum ${\bm{K}}$ is diagonalized in $S^{z}$ subspaces. We do not rotate the Hamiltonian matrix to a ${\bm{S}}$ basis as this is an additional computational cost, and instead determine the ground state total total spin $S$ by identifying multiplets from the $S^{z}$-dependent energy eigenvalues. The total spin $S$ assignments have been confirmed by calculating the spin structure factor $\mathcal{S}({\bm{q}}=0)$. For a given momentum, the $S$-multiplet structure implies that the largest subspace corresponds to the lowest possible $S^{z}$ which contains states with all values of total spin ${\bf S}$.

All of our exact diagonalization calculations are limited only by the maximal matrix size of a given subspace, with the largest subspace corresponding to the lowest possible $S^{z}$. The results presented here correspond to systems containing $M=9,12$ and $16$ moiré unit cells. Despite the limited system sizes that can be reached with the exact diagonalization method, important information can still be extracted concerning the behaviors of charge gaps and the nature of the magnetic order of insulating states. Numerical non-perturbative approaches, like the one taken in this paper or DMRG, are particularly important in the quantum melting regime, where Hartree-Fock is known to overestimate the stability of the insulating phase.

We can identify three main regimes in our phase diagram, Fig. 2(a)

1. 1.

   Atomic limit – $W/V_{m}\ll 1$. This limit corresponds to a perfectly flat band with $t=0$, equivalent to a classical electron gas on a lattice created by the strong moiré potential. Electrons are expected to localize on a subset of superlattice sites so as to minimize on-site and near-neighbor interactions, giving rise to generalized Wigner crystal states. Classical Monte Carlo [17, 7] and Hartree-Fock [25] calculations can address the specific form of the charge order, which depends on the filling. We see in Fig. 2(b) that the charge gaps at all multiples of $1/3$ and $1/4$ are equal in the atomic limit. The many-body ground state manifold contains a large number of nearly degenerate states, corresponding to different spin states on the same sublattice of occupied sites. As the atomic limit is approached, the spacings between these levels become too small to allow numerical determination of the ground state magnetic order. The small spacing of these levels implies that the full entropy of the spin subspace will be realized at a low temperature.

2. 2.

   Melting regime – $W/V_{m}\sim 0.05-0.1$. In the intermediate regime, the competition between electron localization due to long-range Coulomb interaction and quantum fluctuations controls the generalized Wigner crystal melting. Quantum fluctuations of charge distributions are enabled by the hopping term, $t$. When $t$ increases from the atomic limit, interactions between spins located on different sites strengthen. The ground state spin manifold then broadens sufficiently to allow the magnetic properties of some crystal states to be addressed. Interestingly, depending on the dielectric constant value $\varepsilon$ and the filling fraction $\nu$, we can obtain either ferromagnetic or antiferromagnetic states.

3. 3.

   Electron gas limit – $W/V_{m}\gg 1$: This is the limit of weak moiré modulation. The system will resemble a 2DEG and the ground state is determined by the electron density. The question of how Wigner crystallization in the 2DEG is modified by adding an underlying lattice as a small perturbation is not a trivial one and the properties expected for electron crystals in the 2DEG are significantly modified. We obtain vanishing charge gaps at all fillings corresponding to Fermi liquid states as can be seen in Fig. 2(a) (in Appendix B we show examples of occupation distributions where the Fermi surfaces can be identified). This is a limit where Hartree-Fock calculations do not give reliable results [26].

In focusing on filling fractions $\nu=1/3$ and $\nu=2/3$, we will mainly discuss results for a system containing $M=9$ moiré unit cells in the main text and present some additional results for $M=12$ in Appendix B. Despite the small system size, the $M=9$ geometry captures the charge distribution observed in experiments [9, 27], hence it can give us some insight into how increasing hopping amplitudes from the localized limit initially determines the magnetic ground state and then, ultimately, drives melting.

In order to determine the nature of the magnetic ground states, translational symmetry breaking can be tested by evaluating the static spin-spin correlation function

and using it to calculate a static spin structure factor defined as

In these equations $\bm{q}$ is an extended zone wavevector, $A$ is the sample area, and the expectation values are taken in the many-body ground state. Finite size peaks in this quantity at wavevectors that are not reciprocal lattice vectors indicate broken translational symmetry. The inverse Fourier transform of $\mathcal{S}({\bm{q}})$,

is used below to characterize the correlations between spins at positions separated by ${\bm{r}}$.

Additional insight is provided by the number of finite-size many-body eigenvalues in different energy ranges. The Hilbert space can be divided into a singly-occupied subspace (low energy sector) and a doubly-occupied subspace (high energy sector) analogous to the lower and upper Hubbard bands for half-filling. Since the spectrum connects adiabatically to the atomic limit, the number of states with energy smaller than $\sim U_{0}$ is equal to the dimension of the full single-occupancy Hilbert space. For filling $\nu=N/M$, the subspace of the single-occupancy Hilbert space with $N_{\uparrow}$ spins up and $N_{\downarrow}$ spins down has dimension given by

By considering all possible configurations $(N_{\uparrow},N_{\downarrow})$ we obtain the total dimension of the single occupancy, low-energy, sector.

In Fig. 3(a) we illustrate the spatial configuration of the ground state for $\nu=1/3$. In the atomic limit ($t=0$) all spin configurations are degenerate and the ground state charge distribution corresponds to a triangular lattice with periodicity $\sqrt{3}\,a_{M}$ (represented as green sites). This is the only charge configuration that avoids the $U_{1}$ nearest neighbor interaction scale. This charge order has been measured by STM [9, 27]. Other theoretical studies [18, 19, 20] have also found the charge configuration of Fig. 3(a) in the highly localized limit.

Starting from the atomic limit, there is a regime of finite $t$ in which the charge order pattern and the insulating state gap survive, but quantum fluctuations induce interactions between localized spins that lift the large ground state spin degeneracy. From our ED results we obtain a ground state with minimum total spin, indicating an antiferromagnetic state. The magnetic configuration in this regime is also shown in Fig. 3(a). We confirm that spins form a $120^{\circ}$-Néel antiferromagnetic state by evaluating the structure factor $\mathcal{S({\bf q})}$ and its inverse Fourier transform $\xi({\bm{r}})$, presented in Fig. 3(d) and (e) respectively. We observe structure factor peaks at the middle inner points ${\bm{\vartheta}}$ of the Brillouin zone, indicating a magnetic unit cell with 9 sites. From the correlation function $\xi({\bm{r}})$ we see that all first neighbors of an occupied site are empty, while the spin orientations on the second neighboring sites form an angle larger than $\pi/2$ with respect to the occupied site, which translates into negative values. In the limit of classical spins, for this triangular Néel state, evaluation of the structure factor at high symmetry points results in a finite value only at ${\bm{\vartheta}}$

When quantum fluctuations are not too strong, these classical estimates are still approximately valid [28]. Therefore the appearance of peaks at ${\bm{\vartheta}}$ in Fig. 3(d) is in agreement with the indicated spin configuration.

The many-body spectrum separates into a high-energy sector associated with double occupations and characteristic energy $U_{0}$ and a low-energy sector of configurations with no double occupations. For $N=3$ particles on $M=9$ sites, the full Hilbert space has dimension 816 and the low-energy sector of single-occupancy consists of 672 states. Fig. 3(c) shows the many-body energies as a function of $W/V_{m}$, obtained by ED, corresponding to the low-energy single-occupancy Hilbert space, that in turn separates into three bands.

The lowest band, or ground state manifold, contains 24 states that are degenerate in the atomic limit. This number can be understood by noting that there are three inequivalent Wigner crystal configurations at $\nu=1/3$–filling, corresponding to choosing one of the three sites to occupy within the unit cell. For each of these states there is a multiplicity of $2^{3}$ when the spin degree of freedom is taken into account. As $t$ is increased (right limit of the plot) the degeneracy is lifted. In contrast to the half-filled case, the $\nu=1/3$ state has a branch of excitonic states with energies that lie below the charge gap, indicated by a dark green line with dots, as can be observed in our spectrum. The number of these particle-hole excitations is $N_{ex}=432=24\times 3\times 6$, corresponding to the product of the number of states in the ground state manifold, the number of electrons, and the number of neighboring empty sites around each filled site. Finally, the third band in Fig. 3(c) is the manifold of itinerant charged excitations. Its minimum coincides with the charge gap $\Delta_{c}$ in the atomic limit. This band has 216 states for the system size considered, corresponding to configurations where the three particles occupy neighboring sites. The multiplicity of the itinerant charged excitation branch grows most quickly as the system size is increased.

From the mapping of the TMD continuum model to an extended triangular Hubbard model [1, 12], we have derived an effective spin model that describes the magnetic ground state we observe. The Heisenberg interaction between localized spins in the configuration shown in Fig. 3(a), up to order $t_{1}^{4}/U_{0}U_{1}^{2}$, is

where $t_{1}$ and $t_{2}$ stand for the first and second nearest-neighbor hopping parameters, $X_{2}$ is the second-neighbor direct exchange and we have defined an effective second-neighbor hopping $t_{\text{eff}}=t_{2}-t_{1}^{2}/U_{1}$. Fig. 3(a)-(b) show a schematic of the virtual state processes that contribute to $t_{\text{eff}}$. This virtual second-neighbor hop process does give a contribution to the effective spin model, Eq. (10), that is similar to the analogue half-filling $t_{1}$ process.

Hopping of an arbitrary particle to a neighboring site by $t_{1}$, increases the energy by $2U_{1}$ (we neglect the long range part of Coulomb interaction in this analysis, for simplicity). Note that terms of order $t_{1}^{2}$ do not influence spin interactions, but repeated $t_{1}$ and $t_{2}$ hops without double-occupying a site yield the second and third terms in Eq. (10). The value of $J_{1}(1/3)$ as a function of interaction strength is plotted in Fig. 5(a) below.

The spatial configuration of the ground state for $N=6$ particles on $M=9$ sites (for the same parameters as those of Fig. 3) is shown in Fig. 4(a). In the atomic limit the charge distribution forms a honeycomb lattice, in agreement with recent STM measurements [9]. When quantum fluctuations are turned on, the ground state has antiferromagnetic order, consistent with the measured Curie-Weiss (CW) temperature indicator [14] on a similar heterobilayer system.

The structure factor, shown in 4(d), has peaks at ${\bm{\kappa}}$, indicating a magnetic unit cell with three sites, while also showing smaller non-zero values for the middle points ${\bm{\vartheta}}$. The function $\xi({\bm{r}})$ is illustrated in Fig. 4(e). The second neighbor sites have a positive value indicating same spin orientation, while the six nearest neighbors of an occupied site have negative values. We interpret these small negative values as the average of three occupied sites with opposite spin orientation and three empty sites. These spin correlation function results seem to indicate that at $\nu=2/3$–filling, there is greater virtual occupation of empty sites compared with the $\nu=1/3$–filling case. A classical analysis of the structure factor for the spin configuration in Fig. 4(a) yields the following values for high-symmetry points

allowing us to conclude in favor of the state in Fig. 4(a) when quantum fluctuations are not strong.

For $N=6$ particles on $M=9$ sites, the full Hilbert space has dimension 18564, while the dimension of the single-occupancy Hilbert space is 5376. Fig. 4(c) shows the low-energy many-body spectrum, which separates into three bands as in the $\nu=1/3$ case, with a higher-energy charged sector, and a mid-energy region containing excitonic states with bound electron-hole pairs, indicated by red arrows in Fig. 4(a)-(b), and the ground state manifold. The ground state band has 192 states, corresponding to three inequivalent Wigner crystal configurations (the conjugate configurations of the $\nu=1/3$ case) and a multiplicity of $2^{6}$ for each, when the spin degree of freedom is accounted for. The electron-hole bound pair sector has $N_{ex}=3456=192\times 6\times 3$ states, which result from the product of the number of states in the ground state manifold, the number of electrons and the number of empty sites where each electron can hop. The higher band contains the charged excitations and has 1728 states, which correspond to the conjugates of the charge distributions forming the higher band in the $\nu=1/3$ case, times the spin multiplicity $2^{6}$. In this case the Heisenberg coupling constant of the effective spin honeycomb model for the configuration in Fig. 4(a), at lowest order, is

where $X_{1}$ is nearest neighbor direct exchange. Virtual hopping to a double occupied site by $t_{1}$ is the main process responsible for antiferromagnetism, the first term in Eq. (12), shown by blue arrows in the cartoon of Fig. 4(a)-(b). The value of $J_{1}(2/3)$ as a function of interaction strength is plotted in Fig. 5(b).

The magnetic honeycomb pattern found at $\nu=2/3$, shown in Fig. 4(a), is sensitive to model parameters. If the value of the dielectric constant $\varepsilon$ is increased, the ground state transits from antiferromagnetic to ferromagnetic. The dielectric function increase can be accomplished experimentally by modifying the substrate or varying the distance from active device to electrical gates. This change of properties is illustrated in Fig. 5, where we present structure factors for two values of $\varepsilon$ for both fillings $\nu=1/3$ and $2/3$. (We also include results corresponding to $M=12$.) For $\nu=1/3$, shown in Fig. 5(a), the main peaks are at the middle points, ${\bm{\vartheta}}$, for the two system sizes and for both values of $\varepsilon$, in agreement with the state of Fig. 3(a). On the other hand, structure factors for $\nu=2/3$ change qualitatively from $\varepsilon=20$ to $\varepsilon=10$, as seen in Fig. 5(b). In particular, we observe peaks at ${\bm{\gamma}}$ in the structure factors for $\varepsilon=10$ in both system sizes, characteristic of a ferromagnetic state. Smaller peaks at ${\bm{\kappa}}$ can also be observed, indicating that the system still breaks the moiré lattice translation symmetry. For $\varepsilon=20$, the main peaks in the structure factor are at ${\bm{\kappa}}$, in agreement with the state depicted in Fig. 4(a).

The behavior at $\nu=2/3$ can be qualitatively understood in terms of the effective spin model of Eq. (12). As the value of $\varepsilon$ is decreased, the contribution from direct exchange starts to dominate, eventually turning the coupling constant negative and driving the ground state ferromagnetic. We show in Fig. 5(b) that $J_{1}(\nu=2/3)$ changes sign as $\varepsilon$ decreases. For the other filling, $\nu=1/3$, the effective spin coupling $J_{1}(\nu=1/3)$ remains positive for all values of $\varepsilon$ considered, as seen in Fig. 5(a). However, because the main processes determining the magnetic properties have a much smaller energy scale for $J_{1}(\nu=1/3)$, it is almost one order of magnitude smaller than $J_{1}(\nu=2/3)$ at the same $V_{m}$. This trend agrees with a recent experiment by the Cornell group [14], where the Curie-Weiss temperature measurements for $\nu=2/3$–filling clearly confirmed antiferromagnetic order, while for $\nu=1/3$–filling the result is inconclusive due to the very small values obtained. The experiment also found that the Curie-Weiss temperature vs. $\nu$ has a strong local minimum at $\nu=2/3$, suggesting a competition between super-exchange and direct exchange. Our simple finite size calculations allow us to reproduce this phenomenology. We note that in [14] a similar effective honeycomb spin model with positive first-neighbor coupling but negative second-neighbor coupling was proposed, and it can also give rise to a ferromagnetic groundstate.