Kanamori-Moiré-Hubbard model for transition metal dichalcogenide homobilayers

Introduction.— Transition metal dichalgenide (TMD) moiré materials provide unprecedented platforms to study the effect of electronic correlations on flat moiré bands [1, 5, 2, 3, 4]. A variety of low-energy Hamiltonians can be realized in these TMD moiré materials [6]. For example, the $\mathrm{WSe_{2}/WS_{2}}$ heterobilayer simulates the one-orbital triangular lattice Hubbard model [7, 8, 9, 10], while the AB-stacked $\mathrm{MoTe_{2}/WSe_{2}}$ leads to non-trivial moiré bands demonstrating quantum anomalous Hall effect [11]. In addition, recent ab-initio and continuum model calculations have shown that twisted $\Gamma$-valley homobilayers, such as $\mathrm{MoS_{2}}$, $\mathrm{MoSe_{2}}$, and $\mathrm{WS_{2}}$, produce two valence moiré bands with Dirac cone mimicking a honeycomb lattice, while the next set of lower energy four moiré bands simulates the two-orbital asymmetric $p_{x}$-$p_{y}$ honeycomb lattice model [12, 13, 14, 15]. Moreover, surprisingly in recent ARPES experiments $\Gamma$-valley moiré bands have been observed in the twisted $\mathrm{WSe_{2}}$ [16, 17], rendering it also a candidate material to realize the two-orbital honeycomb lattice model. These findings opens up an exciting avenue to simulate multi-orbital Hubbard-like models in TMD moiré materials.

The Kanamori-Hubbard (KH) model [18, 19] has been extensively studied for many conventional materials where multiple orbitals are active, as in iron based superconductors, iridates, manganites. etc. [20, 21, 22, 23, 24]. The moiré potential is shallower than the ionic potential present in conventional materials, leading to relatively broader Wannier functions in moiré materials and making non-local correlations important [25], which are typically ignored in the often used KH model. This suggests that the theoretical studies of twisted $\Gamma$-valley homobilayers require a model going beyond the standard KH model. In this publication, for the first time we provide a Kanamori-Moiré-Hubbard (KMH) model which includes non-local correlations, where the interaction parameters are calculated using the well-localized and accurate Wannier functions of the twisted $\mathrm{MoS_{2}}$ bilayer [26, 27, 28]. The importance of the KHM model is depicted by discussing the effective dielectric constant $\epsilon$ vs the twist angle $\theta$ phase diagrams for the half-filled KMH model, unveiling surprising results which definitely cannot be captured by the standard KH model. It is interesting to note that the relevance of the non-local correlations in the flat moiré bands of twisted bilayer graphene (TBG) has also been discussed [29, 30, 31], so we believe the KMH model can also be used for TBG, but only near magic angles [32, 33, 34, 35] unlike in TMD bilayers where the flat bands are present in a larger range of twist angles.

Wannierization and tight-binding model.— We calculate the moiré bands structure and the Bloch states using the continuum moiré Hamiltonian $H=-\hbar k^{2}/2m^{*}+\Delta({\bf{r}})$. The moiré potential $\Delta(\bf{r})$ is defined as $\Delta({\bf{r}})=\sum_{s}\sum_{j=1}^{6}V_{s}e^{i({\mathbf{g}}_{j}^{s}.{\mathbf{r}}+\phi_{s})}$, where ${\bf g}_{j}^{s}$ are the moiré reciprocal lattice vectors connecting to $s$-th nearest neighbour. The model parameters $\{V_{s},\phi_{s}\}$ are fixed following earlier studies [12], considering the $\mathrm{MoS_{2}}$ homobilayer, so that the band structure obtained from continuum model and ab-initio matches very well. All of our predictions are also valid for other $\Gamma$ valley homobilayers like $\mathrm{MoSe_{2}}$ and $\mathrm{WS_{2}}$. The valence bands closest to the chemical potential can be described by a one-orbital tight-binding model on a honeycomb lattice, see [12, 13]. Here, we focus on the second-set of 4 composite valence bands, which can be described by a two-orbital $p_{x}$-$p_{y}$ tight-binding model on the honeycomb lattice. Until now, the Wannier functions have not been calculated for these set of bands. We perform Wannierization, using projection technique [36, 37], to obtain 4 well-localized Wannier functions, two on each sublattice namely $A$ and $B$, see Fig. 1(k) (for details see supplementary [38]). The calculated Wannier fuctions have nodes at the moiré sites and a pair of lobes like in the $p$-orbitals of the hydrogen atom, as shown in Fig. 1(c-f) and Fig. 1(g-j) for twist angles $1^{\circ}$ and $2.5^{\circ}$, respectively. We noticed that $\Psi_{A(B)p_{x}}({\bf{r}})$ cannot be obtained by a $90^{\circ}$ rotation of $\Psi_{A(B)p_{y}}({\bf{r}})$ unlike in the ideal $p_{x}$-$p_{y}$ orbitals, which follows from the absence of full rotational symmetry in the moiré potential. Moreover, we found $\Psi_{Bp_{x(y)}}({\bf{r}})=-\Psi_{Ap_{x(y)}}(-{\bf r})$ due the inversion symmetry of the moiré potential on two sublattices given by $\Delta({\bf r}-{\bf R}_{A})$=$\Delta(-{\bf r}-{\bf R}_{B})$.

Using the above Wannier functions, we calculated the hopping parameters for the two-orbital tight-binding model on the honeycomb lattice, up to third nearest-neighbour using $t_{{\bf j}-{\bf i}}^{SS^{{}^{\prime}}}(\mu,\nu)=\langle\Psi_{S\mu}^{\bf j}|H|\Psi_{S^{{}^{\prime}}\nu}^{\bf i}\rangle$ (for details see Supplementary Material [38]), where the $\{{\bf i},{\bf j}\}$, $\{S,S^{{}^{\prime}}\}$, and $\{\mu,\nu\}$ indices denotes unit-cell, sublattice, and orbitals ($p_{x}$ or $p_{y}$), respectively. We write the kinetic energy as $H_{\textrm{K.E.}}=\sum_{\bf{i}\sigma}K_{\bf{i}\sigma}^{1}+K_{\bf{i}\sigma}^{2}+K_{\bf{i}\sigma}^{3}$, where the terms $K_{\bf{i}\sigma}^{n}$ consists of hoppings between the $n^{th}$ nearest neighbour sites in the honeycomb lattice. The hopping connections up to the 3rd nearest-neighbour are pictorially shown in Fig. 1(k). $K_{\bf{i}\sigma}^{1}$ is presented below:

$$K_{\bf{i}\sigma}^{1}=\sum_{\begin{subarray}{c}\nu,\mu\in\{p_{x},p_{y}\}\\ \bf{r}\in\{0,-\bf{a}_{2},\bf{a}_{1}-\bf{a}_{2}\}\end{subarray}}t_{\bf{r}}^{BA}(\mu,\nu)c_{{\bf{i}+\bf{r}}B_{\mu}\sigma}^{\dagger}c_{{\bf{i}}A_{\nu}\sigma}+h.c.$$

(1)

The $K_{\bf{i}\sigma}^{2(3)}$ terms can be written similarly, as shown in supplementary [38]. The 1st nearest-neighbour hopping term $K_{\bf{i}\sigma}^{1}$, shown in eq. 1, depends on three 2$\times$2 matrices namely $\{t_{\bf{0}}^{BA},t_{-\bf{a}_{2}}^{BA}$, $t_{\bf{a}_{1}-\bf{a}_{2}}^{BA}\}$. Similarly six 2$\times$2 matrices $\{t_{\bf{a}_{1}}^{AA},t_{-\bf{a}_{2}}^{AA},t_{\bf{a}_{1}-\bf{a}_{2}}^{AA},t_{\bf{a}_{1}}^{BB},t_{-\bf{a}_{2}}^{BB},t_{\bf{a}_{1}-\bf{a}_{2}}^{BB}\}$ and three 2$\times$2 matrices $\{t_{\bf{a}_{1}}^{AB},t_{\bf{a}_{1}}^{BA},t_{\bf{a}_{1}-2\bf{a}_{2}}^{BA}\}$ are required for the 2nd and 3rd nearest-neighbour hoppings, respectively. All of these 12 matrices are dependent on $\theta$. We found a good match between the band-structure calculated using the above tight-binding model and the continuum model, as shown in Fig. 1(a,b), suggesting that we have accurate Wannier functions. We noticed that for $\theta\lessapprox 1.2$, only nearest-neighbour hoppings are enough to obtain the correct band-structure, as shown in Fig. 1(a) for $\theta=1.0$. However, for larger $\theta$ longer-range hoppings are required to reproduce the continuum model results (see Fig. 1(b), for $\theta=2.5$). We show the evolution of the some dominant hopping parameters in Fig. 1(l), depicting the exponential fast growth of hoppings with $\theta$.

Interaction parameters and Kanamori-Moiré-Hubbard model.— Now we will derive the Coulomb interaction between the fermions in the Wannier states discussed above. The generic interaction term can be written as:

$$H_{\textrm{Int}}=1/2\sum_{\begin{subarray}{c}\mathbf{i},\mathbf{j},\mathbf{k},\mathbf{l},\\ \alpha,\beta,\gamma,\delta,\\ \sigma,\sigma{{}^{\prime}}\end{subarray}}V_{\mathbf{ijkl}}^{\alpha\beta\gamma\delta}c_{\bf{i}\alpha\sigma}^{\dagger}c_{\bf{j}\beta\sigma{{}^{\prime}}}^{\dagger}c_{\bf{l}\delta\sigma{{}^{\prime}}}c_{\bf{k}\gamma\sigma},$$

(2)

where $V_{\mathbf{ijkl}}^{\alpha\beta\gamma\delta}=\langle\Psi^{\bf{i}}_{\alpha}\Psi^{\bf{j}}_{\beta}|V|\Psi^{\bf{k}}_{\gamma}\Psi^{\bf{l}}_{\delta}\rangle$ and $V=e^{2}/\epsilon|{\bf r}_{1}-{\bf r}_{2}|$. $\epsilon$ is produced by the surrounding dielectric enviroment, such as nearby h-BN layers. The exact value of $\epsilon$ is not known so we keep it as a free parameter. The $\{\alpha,\beta,\gamma,\delta\}$ indices represent the sublattice $S$ and the orbital $\mu$ via $\alpha=2S+\mu=S_{\mu}$, where the sublattice $A(B)=0(1)$ and the orbital $p_{x}(p_{y})=0(1)$.

In the present work, for simplicity, we limit the non-local Coulomb interactions only up to nearest-neighbour sites of the honeycomb lattice. A priori, the longer range interactions are not expected to be very relevant at and near half-filling [39]. To study Wigner crystals at fractional fillings, the approximate longer range interactions can be easily included by assuming the $(\frac{1}{|r|}-\frac{1}{\sqrt{r^{2}+d^{2}}})$ functional form, where $d$ is the screening length [40, 41, 42]. The Coulomb interaction term which includes up to nearest-neighbour interactions can be divided into three parts, $H_{\textrm{Int}}=\sum_{\bf{i}}H_{\bf{i}}+H_{\bf{i},{\bf{i}\textrm{-}}{\bf{a}_{2}}}^{AB}+H_{\bf{i},{\bf{i}\textrm{+}\bf{a}_{1}}{\textrm{-}}{\bf{a}_{2}}}^{AB}$, where $\bf{i}$ is the unit cell index and $\bf{a_{(1,2)}}$ are the Bravais lattice vectors. The first part $H_{\bf{i}}$ consists of all the Coulomb interactions possible within the unit cell $\bf{i}$, including both local and nearest-neighbour interactions given by $V_{{\bf{iiii}}}^{\alpha\beta\gamma\delta}$ (total $4^{4}$ terms). The second $H_{\bf{i},{\bf{i}\textrm{-}}{\bf{a}_{2}}}^{AB}$ and third parts $H_{\bf{i},{\bf{i}\textrm{+}\bf{a}_{1}}{\textrm{-}}{\bf{a}_{2}}}^{AB}$ contains the Coulomb interactions between the nearest-neighbour sites belonging to different unit cells. Now we will discuss the $H_{\bf{i}}$ term in detail; the other two terms are very similar and shown in the supplementary. $H_{\bf{i}}$ is shown in Eq. 3, where ${\bf{S}}_{{\bf{i}}\alpha}$=$\frac{1}{2}\sum_{s,s^{{}^{\prime}}}c_{{\bf i}\alpha s}^{\dagger}{\bm{\tau}}_{ss^{{}^{\prime}}}c_{{\bf i}\alpha s^{{}^{\prime}}}$ represent the spin at unit cell ${\bf i}$, orbital $\mu$=mod$(\alpha,2)$, and sublattice=$(\alpha-\mu)/2$. The pair anhilation operator is defined as $P_{\bf{i}\alpha}=c_{\bf{i}\alpha\downarrow}c_{\bf{i}\alpha\uparrow}$. $s=1(-1)$ for $\sigma=\uparrow(\downarrow)$, and the set $\mathbb{S}=\{\{0123\},\{0132\},\{0213\}\}$.

$$H_{\bf{i}}=U_{0}\sum_{\alpha}n_{\bf{i}\alpha\uparrow}n_{\bf{i}\alpha\downarrow}+\sum_{\alpha<\beta}(U_{\alpha\beta}-\frac{J_{\alpha\beta}}{2})n_{\bf{i}\alpha}n_{\bf{i}\beta}\\ -2\sum_{\alpha<\beta}J_{\alpha\beta}{\bf S}_{\bf{i}\alpha}\cdot{\bf S}_{\bf{i}\beta}+\sum_{\alpha<\beta}J_{\alpha\beta}(P_{\bf{i}\alpha}^{\dagger}P_{\bf{i}\beta}+h.c.)\\ +1/2\sum_{\begin{subarray}{c}\sigma,\sigma{{}^{\prime}},\alpha\neq\beta\neq\gamma\end{subarray}}(A_{\beta\alpha\gamma}-\delta_{\sigma\sigma{{}^{\prime}}}\tilde{J}_{\beta\alpha\gamma})(c_{\bf{i}\alpha\sigma}^{\dagger}c_{\bf{i}\gamma\sigma}n_{\bf{i}\beta,\sigma{{}^{\prime}}}+h.c.)\\ +\sum_{\sigma,\alpha\neq\beta}\tilde{A}_{\alpha\beta}(c_{\bf{i}\alpha\sigma}^{\dagger}c_{\bf{i}\beta\sigma}n_{\bf{i}\beta\bar{\sigma}}+h.c.)\\ -\sum_{\alpha\neq\beta\neq\gamma}\tilde{J}_{\alpha\gamma\beta}(S_{\bf{i}\alpha}^{+}c_{\bf{i}\beta\downarrow}^{\dagger}c_{\bf{i}\gamma\uparrow}+h.c.)\\ +1/2\sum_{\begin{subarray}{c}\sigma,\alpha\neq\beta\neq\gamma\end{subarray}}\tilde{J}_{\alpha\gamma\beta}s(P_{\bf{i}\alpha}^{\dagger}c_{\bf{i}\gamma\bar{\sigma}}c_{\bf{i}\beta\sigma}+h.c.)\\ +\sum_{\begin{subarray}{c}\sigma,\sigma{{}^{\prime}},\\ \{\alpha\beta\gamma\delta\}\in{\mathbb{S}}\end{subarray}}T_{\alpha\beta\gamma\delta}(c_{\bf{i}\alpha\sigma}^{\dagger}c_{\bf{i}\gamma\sigma}(c_{\bf{i}\beta\sigma{{}^{\prime}}}^{\dagger}c_{\bf{i}\delta\sigma{{}^{\prime}}}+c_{\bf{i}\delta\sigma{{}^{\prime}}}^{\dagger}c_{\bf{i}\beta\sigma{{}^{\prime}}})+h.c.)$$

(3)

Equation 3 encompasses all $4^{4}$ intra-unit cell interaction terms. The first four terms look similar to the conventional multiorbital Hubbard model, but here they capture the non-local interactions as well. This is the first time such a model is shown.

The first term is the standard onsite intra-orbital Hubbard repulsion, where $U_{0}=V_{\bf{iiii}}^{\alpha\alpha\alpha\alpha}$ (same for all $\alpha$’s). The second term incorporates the onsite inter-orbital density-density repulsions via parameters {$U_{01}$, $U_{23}$, $J_{01}$, $J_{23}$} and the non-local orbital resolved repulsions via parameters like $U_{02}$, $J_{02}$, etc., where $U_{\alpha\beta}=V_{\bf{iiii}}^{\alpha\beta\alpha\beta}$ and $J_{\alpha\beta}=V_{\bf{iiii}}^{\alpha\beta\beta\alpha}$. The well known local Hund’s coupling is present in the third term via the dominant $J_{01}$ and $J_{23}$ parameters; this term also includes the non-local ferromagnetic direct exchange terms ($J_{02}$,$J_{13}$, $J_{03}$, $J_{12}$). The fourth term incorporate the onsite inter-orbital and non-local pair hopping terms. We also found interaction assisted hoppings (term-5 and term-6), spin-flip hopping accompanied with local spin flip (term-7), and scattering of doublon to different states (term-8) quantified by ($A_{\beta\alpha\gamma}$, $\tilde{A}_{\alpha\beta}$, $\tilde{J}_{\beta\alpha\gamma}$). The remaining interactions are present in term-9.

We show the interaction parameters of first 6 terms as a function of $\theta$ in Fig. 2. The density-density terms are dominant interactions, see Fig. 2(a). The onsite intraorbital repulsion ($U_{0}$) suggests that $\epsilon U_{0}/W$ can be of order of 10 to 1000 in real materials, depending on $\theta$, where $W$ is the non-interacting bandwidth. For example, $U_{0}/W$ is about 1200$\epsilon^{-1}$ and 25$\epsilon^{-1}$ for $\theta$=$1^{\circ}$ and $\theta$=$2.5^{\circ}$, respectively. The local Hund’s coupling and the non-local ferromagnetic direct exchange is shown in Fig. 2(b). Fig 2(c,d) displays the interaction assisted hoppings vs. $\theta$. The rest of the interaction parameters are relatively smaller, and shown in the supplementary. We call the total Hamiltonian $H=H_{\textrm{K.E.}}+H_{\textrm{Int}}$ the Kanamori-Moiré-Hubbard (KMH) model because of the presence of non-local interaction terms, which are ignored in the standard KH model. These non-local correlations can lead to unexpected results, as shown in the next sextion. It should be noted that the KMH model shown here has larger scope and can be also used for magic-angle TBG and future moiré materials addressing multiorbital physics on honeycomb lattice (only the values of hopping and interaction parameters will depend on the specific material).

Numerical results at half-filling.— We create $\epsilon$ vs $\theta$ phase diagrams to investigate the physics of the KMH model at half-filling $n$=$N/L$=$2$, where $N$ is the total number of fermions and $L$=$L_{1}$$\times$$L_{2}$ the total number of unit cells. We studied 6$\times$6 and 12$\times$12 system sizes using the unrestricted Hartree-Fock technique. We choose a broad range of $\epsilon$ $\in[1,80]$ as it can be tuned by changing the distance with the nearby metallic gate. Moreover, $\epsilon$ will be enhanced by the charge-fluctuations between the moiré bands considered here and other remote moiré bands. The $\epsilon$ vs $\theta$ phase diagram for the KMH model is shown in Fig. 3(a). Surprisingly, in addition to the expected antiferromagnetic (AFM) state, we have unveiled two new states not anticipated to be stable: the S=1 ferromagnetic (FM) state for $\theta<1.75$ and the charge density wave (CDW) state for $\theta\geq 1.75$. The non-local density-density repulsion plays the key role to stabilize the CDW state. The competition between the non-local FM direct exchange ($\propto\epsilon^{-1}(J_{02}+J_{13}+2J_{03})$) and the AFM superexchange ($\propto(\epsilon t^{2})/(U_{0}+J_{01})$) leads to the transition from FM to AFM state as $\epsilon$ increases. We found that the AFM state is present only for $\epsilon>20$ with local moment $\mathrm{S}<1$. See Fig. 3(c,d,e) for the pictorial represention of the FM, AFM, and CDW states.

We also used the simplified KMH model, only keeping the first 4 terms in Eq. 3, and found all three phases are present nearly in the same region of the phase diagram (see Fig. 3(b)), suggesting that the FM direct exchange and the density-density repulsion are the most important non-local interactions for the half-filled KMH model.

To investigate the effect of the quantum fluctuations, we used the Lanczos technique and studied a small $2\times 2$ cluster with periodic boundary conditions (Fig. 4(b)), using the simplified KMH model. The phase diagram is shown in Fig. 4(a). We again found the FM, AFM, and CDW states in the same region of the phase diagram. Fig. 4(c) shows the spin-spin correlation, with respect to site=1, $\langle{\bf S}_{1}\cdot{\bf S}_{j}\rangle$ for $\theta$=$1.5^{\circ}$, depicting strong FM correlations for $\epsilon$=10, and AFM correlations for $\epsilon\in\{30,50,80\}$. Fig. 4(e) and Fig. 4(f) show $\langle{\bf S}_{1}\cdot{\bf S}_{j}\rangle$ and density-density correlations $N_{1j}=\langle n_{1}n_{j}\rangle-\langle n_{1}\rangle\langle n_{j}\rangle$, respectively, at fixed $\theta$=2.0 depicting the growth of AFM correlations and suppression in the CDW as $\epsilon$ increases. This indicates a smooth transition from the CDW phase to the AFM phase. We believe larger systems are required to confirm whether it is a 2nd-order phase transition or a crossover. The FM to AFM or the FM to CDW are first order transitions because the total spin suddenly changes from $2L$=$8$ to $0$. The averaged local moment ${\bf S}^{2}_{L}=(1/4L)\sum_{i\alpha}\langle{\bf S}_{i\alpha}^{2}\rangle$ as a function of $\epsilon$ is shown in Fig. 4(d). We found, for $\theta<1.75$ that ${\bf S}^{2}_{L}$ decreases as $\epsilon$ is increased while the system transits from the S=1 FM to AFM phases, whereas for $\theta\geq 1.75$, ${\bf S}^{2}_{L}$ grows with $\epsilon$ developing AFM correlations with weak CDW.

Conclusions.— We showed that the twisted $\Gamma$-valley TMD bilayers contains physics beyond the conventional multi-orbital Hubbard model. We provide a KMH model which can be used to theoretically study the multi-orbital physics of TMD bilayers. Using our numerical studies at half-filled moiré bands we show that the non-local direct-exchange terms and density-density interactions can lead to S=1 FM insulators and CDW states, respectively, depending on $\epsilon$ and $\theta$. The AFM state can also be obtained but at large $\epsilon>20$. Our theoretical prediction of a S=1 FM insulator can be verified by measuring the magnetic susceptibility and Weiss constant in real materials [7, 43], and the charge ordered state can be observed using high-resolution scanning tunneling experiments [44]. The KMH model can also be used for further theoretical investigations like doping near half-filled correlated insulators and for studying Mott-Wigner crystals at fractional fillings by including longer range density-density interactions.

N. Kaushal and E. Dagotto were supported by the US Department of Energy, Office of Science, Basic Energy Sciences, Materials Science and Engineering Division.