Singlet, Triplet and Pair Density Wave Superconductivity in the Doped Triangular-Lattice Moiré System

Introduction.—Moiré bilayer systems have attracted great attention over the last few years due to their high tunability and capacity to host a wealth of exotic states of matter Andrei et al. (2021); Balents et al. (2020); Kennes et al. (2021). Since the discovery of superconductivity (SC) and Mott insulating phase in magic-angle twisted bilayer graphene (TBG) Cao et al. (2018a, b), other Moiré systems have been realized and are under active studies Liu et al. (2020); Chen et al. (2019), including twisted bilayer transition metal dichalcogenides (TMDs) Zhang et al. (2020a); Shabani et al. (2021); Weston et al. (2020); Devakul et al. (2021); An et al. (2020); Naik and Jain (2018); Zhang et al. (2021); Regan et al. (2020); Schrade and Fu (2019); Zhang et al. (2020b). Compared to TBG, twisted bilayer TMDs have the advantages of accommodating flat Moiré bands over a much wider range of twist angles and fewer low-energy degrees of freedom, allowing for a simpler lattice model description Wu et al. (2018, 2019); Pan et al. (2020). Strong correlation effects such as correlated insulating phase Wang et al. (2020), metal-insulator transition Li et al. (2021a); Ghiotto et al. (2021), stripe phase Jin et al. (2021) and quantum anomalous Hall effect Li et al. (2021b) have recently been observed in these systems.

Twisted TMD bilayers can be classified into hetero- and homo-bilayers according to whether the two layers are made of the same or different materials. The low-energy electronic degrees of freedom in the former are believed to be described by a generalized triangular-lattice Hubbard model with pseudo-spin SU(2) rotation symmetry Wu et al. (2018); Tang et al. (2020), whereas in the latter the spin SU(2) symmetry is broken into U(1) by a vertical displacement field due to spin-valley locking and inversion symmetry breaking, and consequently the electron hopping acquires a spin-dependent phase $\theta$ Wu et al. (2019); Pan et al. (2020); Schrade and Fu (2019); Wang et al. (2023). Note that the standard Hubbard and t-J models on triangular lattices, i.e. $\theta=0$, have exhibited a rich phenomenology enhanced by further-neighbor couplings due to the complex interplay between geometric frustration, quantum fluctuations and hole dynamics Raghu et al. (2010); Jiang (2021); Zhu and Chen (2023); Huang et al. (2023); Wang et al. (2004); Baskaran (2003); Motrunich and Lee (2004); Kumar and Shastry (2003); Chen et al. (2013); Venderley and Kim (2019a); Gannot et al. (2020); Peng et al. (2021). The hopping phase $\theta$ is shown to be widely tunable by the displacement field and thus may serve as a novel control knob of the many-body ground states of twisted TMD homobilayers. The magnetic and superconducting phases under the variation of both carrier density and $\theta$ of the U(1) Hubbard model and/or its closely related t-J model (for strong Hubbard U limit) at/near half-filling have been explored through mean-field calculations, renormalization group analysis, quantum cluster methods and Gutzwiller approximation Zang et al. (2021, 2022); Pan et al. (2020); Zhou and Zhang (2023); Wu et al. (2023a); Bélanger et al. (2022); Zegrodnik and Biborski (2023). However, these methods generally are not accurate in treating the strong electronic correlations present in the model Qin et al. (2022). Here we implement density matrix renormalization group (DMRG) White (1992) to accurately capture the ground states on quasi-1D few-leg cylinders, and thus reveal the different ordering tendencies at play and gain some insights into the physics at the 2D limit Stoudenmire and White (2012); Arovas et al. (2022). Particularly, DMRG has been applied onto a three-leg cylindrical Moiré Hubbard model but only weak SC correlations were observed Wietek et al. (2022). The effective spin-model derived at strong U and half-filling limit was also considered for exploring quantum spin liquid Kiese et al. (2022).

In this work we study SC of the lightly doped triangular lattice U(1) Moiré t-J model on a four-leg cylinder through DMRG calculations. By varying $\theta$ in the region of $(0,\frac{2\pi}{3})$, we identify two conjugated superconducting phases as shown in Fig. 1(b): (i) Mixed spin singlet $d$-wave and triplet $p$-wave SC coexisting with spin, charge and pair density waves (PDW); (ii) Ferromagnetic triplet $p$-wave SC coexisting with charge density wave (CDW). These two phases are related by a combined operation of spin flip and local gauge transformation, up to a change of the boundary condition. Their pairing correlations decay algebraically with the Luttinger exponents smaller or around two, demonstrating a robust quasi-long-range SC order Gong et al. (2021); Jiang and Kivelson (2021). Particularly, distinct from other SC phases on the triangular-lattice t-J model Jiang (2021); Huang et al. (2023); Zhu and Chen (2022), PDW is a novel SC state where Cooper pairs carry finite center-of-mass momentum Agterberg et al. (2020), which are not commonly realized in microscopic models Berg et al. (2010); Wu et al. (2023b, c); Huang et al. (2022); Jaefari and Fradkin (2012); Cho et al. (2012); Lee (2014); Soto-Garrido and Fradkin (2014); Venderley and Kim (2019b); Shaffer and Santos (2023). The plethora of interesting phases found in our calculations could motivate future experimental endeavour in search of novel SC in twisted TMD homobilayers.

Model and Method.— The Moiré t-J model is defined as

where $\sigma=\pm$ represents spin up/down, $c^{\dagger}_{i\sigma}$ and $c_{i\sigma}$ are the creation and annihilation operators for the electron with spin $\sigma$ at the site $i$, $\langle ij\rangle$ denote nearest neighbors whose locations satisfy $\bm{r}_{j}-\bm{r}_{i}\in\{\bm{e}_{a},-\bm{e}_{b},\bm{e}_{c}\}$ (see Fig. 1(a)), $\hat{S}^{z}_{i},\hat{S}^{+}_{i},\hat{S}^{-}_{i}$ are the spin-$\frac{1}{2}$ $\hat{z}$ component, raising and lowering operators at site $i$ respectively, and $\hat{n}_{i}=\sum_{\sigma}\hat{c}^{\dagger}_{i\sigma}\hat{c}_{i\sigma}$ is the electron number operator. Double occupancy is prohibited. The hopping phase $\theta$ produces a flux of $\pm 3\theta$ at each triangular plaquette, and a gauge transformation connects two models differing in the fluxes by $2\pi$. We therefore focus on the region of $0<\theta<2\pi/3$. In the present study, we set the hole doping level $\delta=1/12$, and choose $J=1$ and $t=3$, corresponding to a realistic situation of $U/t=12$ Pan et al. (2020).

To obtain the ground state, we employ DMRG simulation with U(1)$\times$U(1) symmetry corresponding to charge and spin conservation on a cylindrical system with periodic boundary condition (PBC) along the circumferential ($\bm{e}_{b}$ or $y$-) direction and open boundary condition along the axial ($\bm{e}_{a}$ or $x$-) direction. The number of lattice sites is given by $N=L_{x}\times L_{y}$, where $L_{x}$ and $L_{y}$ are the number of sites along $x$- and $y$-direction respectively and are set as $L_{y}=4$ and $L_{x}=36$ in the main text. The corresponding geometry is called YC$L_{y}$ Yan et al. (2011). The doping level is defined by $\delta=1-N_{e}/N$ and we consider the zero total spin-z sector: $\sum_{i}\hat{S}^{z}_{i}=0$, which hosts the ground state as verified in Sec. A of the supplemental materials (SM) SM . In DMRG, the number of Schmidt states kept for representing the reduced density matrix on either side of the system under bipartition is called “bond dimension” $M$ White (1992). The calculations improve with the increase of $M$ and become exact for a sufficiently large $M$.

Coexisting Singlet, Triplet and PDW SC (STPSC).—The SC order is examined by the spin-singlet and triplet pairing correlation functions $P_{\alpha\beta}^{s}(r)$ and $P_{\alpha\beta}^{t_{n}}(r)$ defined by

where the reference point $\bm{r}_{0}\equiv(x_{0},y_{0})=(L_{x}/4,L_{y})$ and the pairing operators $\hat{\Delta}^{s}_{\alpha}(\bm{r}_{1})$ and $\hat{\Delta}^{t_{n}}_{\alpha}(\bm{r}_{1})$ are defined on the bond along $\bm{e}_{\alpha}$ ($\alpha=a,b,c$) at site $\bm{r}_{1}$:

Here $\hat{\Delta}_{\alpha}^{t_{n}}$ corresponds to the triplet pairing with total spin-z $S_{z}=n$.

Fig. 2(a) and (c) show two dominant pairing components: $b$-bond singlet pairing $P^{s}_{bb}(r)$ and opposite-spin-z ($S_{z}=0$) triplet pairing $P^{t_{0}}_{bb}(r)$ for $\theta=\pi/12$ in the STPSC phase. Both exhibit power-law decay $P^{s(t_{0})}_{bb}(r)\sim r^{-K^{s(t_{0})}_{SC}}$ with the Luttinger exponents $K^{s(t_{0})}_{SC}\approx 0.3$, suggesting strongly diverging SC susceptibilities $\chi\sim T^{-(2-K_{SC})}$ as the temperature $T\rightarrow 0$ Arrigoni et al. (2004). Note also that slow power-law decays are already exhibited by the largest-$M$ results with exponents around 0.97. The singlet pairing component is larger in amplitude than the triplet one, and they exhibit $d$-wave and $p$-wave symmetry respectively Raghu et al. (2010); Hsu et al. (2017); Venderley and Kim (2019a). The mixing of singlet and triplet pairings are permitted by the absence of the inversion and spin SU(2) symmetry Yip (2014). In particular, the absence of inversion center allows the mixing of parity-odd $p$-wave and parity-even $d$-wave basis functions in the irreducible representation $E$ of the symmetry group $C_{3v}$ of the system Hsu et al. (2017). The charge density correlation function $D(r)\equiv\langle\hat{n}(\bm{r}_{0})\hat{n}(\bm{r}_{0}+r\bm{e}_{x})\rangle-\langle\hat{n}(\bm{r}_{0})\rangle\langle\hat{n}(\bm{r}_{0}+r\bm{e}_{x})\rangle$ in Fig. 2(b) decays algebraically with a relatively larger exponent (around 0.86), suggesting weaker charge density modulations coexisting with stronger SC. Correspondingly we observe a charge stripe order with two holes per stripe in the inset. For comparison, Fig. 2(d) presents also the in-plane spin-spin correlations $S_{xy}(r)$ defined by

and the Green’s function $G(r)\equiv\sum_{\sigma}\langle\hat{c}^{\dagger}_{\bm{r}_{0},\sigma}\hat{c}_{\bm{r}_{0}+r\bm{e}_{x},\sigma}\rangle$. The in-plane spin correlation is the strongest among all correlations, characterizing a robust spin density wave order inherited from the 2D in-plane $120^{\circ}$ Néel order at half filling based on the spin structure factor calculations Wu et al. (2019); Zang et al. (2021); SM . The Green’s function squared $\lvert G(r)\rvert^{2}$ is much weaker than the main pairing correlations, confirming the dominance of two-electron pairing over single-electron tunnelings.

Moreover, in the $S_{z}=\pm 1$ triplet pairing components, we observe quasi-long-range PDW orders with a Luttinger exponent around 0.58 in Fig. 3(a). The PDW wavevector $\bm{k_{\text{pdw}}}$ can be determined by the variation of the phase of the pairing correlation under displacement along both $\bm{e_{a}}$ and $\bm{e_{b}}$. Specifically,

characterizes spatial variation of the phase of the $b$-bond triplet pairing order parameters. In Fig. 3(b), $\bm{k^{{\pm 1}}_{\text{pdw}}}$ is determined to be $\pm\bm{K^{\prime}}$, which are the nearest accessible wavevectors to the Brillouin zone corners $\pm\bm{K}$ in the YC4 geometry. The same PDW wavevectors are identified for $a$- and $c$-bond. Note that a PDW ground state with $\bm{k^{\pm 1}_{\text{pdw}}}=\mp\bm{K}$ was also predicted for the Moiré Hubbard model at $\theta=\pi/3$ by perturbative renormalization group analysis in the weak coupling regime Wu et al. (2023a, b).

Ferromagnetic Triplet SC (FMTSC).—In the FMTSC phase, we find the dominant pairing channel to be a $p$-wave spin triplet. In Fig. 4(a) and (c), both $P^{t_{0}}_{aa}$ and $P^{t_{1}}_{aa}$ are non-oscillatory, in accordance with uniform SC order in the bulk of the system, and decay algebraically with exponents slightly larger than 2. An accompanying CDW order is confirmed in Fig. 4(c) by both the quasi-long-range density correlation ($\sim r^{-1.75}$) and charge stripes in the electron density profile (one hole per stripe). In Fig. 4(d), a robust in-plane ferromagnetic spin correlation is observed in reminiscence of the parent ferromagnetic order Wu et al. (2019); Zang et al. (2021), with the total spin $S\approx\sqrt{\langle\hat{S}^{2}\rangle}\approx 0.326N_{e}$. The singlet paring is shown much weaker than the triplet ones as the triplet pairing is favored by ferromagnetism. The opposite-spin-z triplet pairing correlation $P^{t_{0}}_{aa}$ has stronger amplitude and slower decay rate than those of the same-spin-z component $P^{t_{1}}_{aa}$ because the ferromagnetic order is in-plane.

Discussion and Summary—The FMTSC and STPSC phases are related by a spin-flip operation followed by a local gauge transformation Zhou and Zhang (2023) as demonstrated in the SM. Particularly, the uniform z-spin-polarized triplet pairing order at $\theta$ in the FMTSC region is conjugated to the PDW order with $\bm{k^{\pm 1}_{\text{pdw}}}=\pm\bm{K}$ at $(2\pi/3-\theta)$ in the STPSC region:

with $\beta_{\alpha}=\bm{e_{\alpha}}\cdot(\bm{b_{1}-b_{2}})/3$, where $\bm{b_{1,2}}$ are the reciprocal wavevectors conjugated to $\bm{e_{a}}$ and $\bm{e_{b}}$ respectively. This is consistent with our observations at $\theta=7\pi/12$ (Fig. 4(c)) and its conjugated partner $\theta=\pi/12$ (Fig. 3), albeit with a different flux (y-boundary phase) into the 4-leg cylinder. Moreover,

which means that the singlet and opposite-spin triplet pairing components are superposed to produce their counterparts in the conjugated phase. Since the singlet pairing component at $\theta=7\pi/12$ is found negligible compared to the triplet components, and $\beta_{\alpha}=2\pi/3$ ($\alpha=a,b$) or $4\pi/3$ ($\alpha=c$), one has $\absolutevalue{\Delta^{s}_{\alpha}(\theta=\pi/12)}\approx\sqrt{3}\absolutevalue{\Delta^{t_{0}}_{\alpha}(\theta=\pi/12)}$ according to Eq. 6, which explains the larger magnitude of the spin singlet pairing than that of the triplet and the same power-law exponents in Fig. 2 (a) and (c).

However, the pairing correlations at $\theta=\pi/12$ has much stronger magnitude (over one order of magnitude larger) and slower decay rate compared to those at $\theta=7\pi/12$ ($K_{\text{SC}}\approx 0.29$ vs. $K_{\text{SC}}\approx 2.27$). This in addition to the difference in charge distributions (two vs. one holes per stripe) is caused by the change of the boundary condition: the periodic boundary condition at $\theta=\pi/12$

turns into a twisted boundary condition Gannot and Kivelson (2023) at $\theta=7\pi/12$

after the gauge transformation, corresponding to inserting a magnetic flux of $\pm 2\pi L_{y}/3$ through the interior of the cylinder for electrons. The spin structure factor of the $120^{\circ}$ Néel order for $0<\theta<\pi/3$ is peaked at $\pm\bm{K}$, which are not resolved in the 4-leg cylinder under PBC, whereas for $\pi/3<\theta<2\pi/3$ the system is ferromagnetic with the peak at the system-supported momentum $\bm{\Gamma}$. Therefore, the former regime is more frustrated than the latter in the YC4 geometry and this might result in stronger SC. The sensitivity of SC to boundary conditions reveal finite-size effects in our four-leg system, so we also study a different cylinder geometry XC4 Szasz et al. (2020) (Sec. H in SM) as well as a YC3 system with $N=40\times 3$ (Sec. F in SM). Both systems preserve the PBC under local gauge transformation and support $\bm{\Gamma}$ and $\bm{\pm K}$ in the Brillouin zone, therefore introducing no frustration. In the XC4 geometry, we again obtain the STPSC and FMTSC phases and their SC correlations now have similar amplitudes and decay with close exponents ($\approx 2$), consistent with Eq. 6. In the YC3 cylinder at $\theta=7\pi/12$, the Luttinger exponents for SC ($K^{t_{0}}_{\text{sc}}\approx 2.28$) is nearly identical to that of the YC4 cylinder ($K^{t_{0}}_{\text{sc}}\approx 2.27$). The observation of quasi-long-range SC order at different boundary conditions, cylinder geometries and sizes is positive evidence for the existence of SC in the 2D limit Szasz et al. (2020).

In contrast with the topological SC phases reported in the mean field and perturbative renormalization group studies of the doped TMD homobilayer Zhou and Zhang (2023); Wu et al. (2023a) or monolayer Hsu et al. (2017); Yuan et al. (2014), both the $d$- and $p$-wave SC phases found here are topologically trivial as the nearest-neighbor pairings acquire a phase of either 0 or $\pi$ after a $\pi/3$ rotation, instead of the nontrivial phases of $\pm\pi/3$ and $\pm 2\pi/3$ for $p\pm ip$ and $d\pm id$-wave topological SC phases Huang and Sheng (2022); Jiang and Jiang (2020); Huang et al. (2023). Furthermore, the SC phase here is distinct from the Ising SC found in electron-doped TMD monolayers Zhou et al. (2016); Lu et al. (2015); Xi et al. (2016); Saito et al. (2016) in that the former arises from hole doping the parent in-plane magnetic Mott insulator at strong electronic couplings whereas the latter the pinning of the electron spins in the Cooper pairs to the out-of-plane directions by the Ising spin-orbit interaction at weak electronic couplings. Finally, the $\theta=\pi/6$ case was also studied in Ref. Wietek et al. (2022), but a rather large power-law decay exponent ($\approx 3.34$) was found, so only weak SC was claimed there. Consistently we find that $\theta=\pi/6$ is located at the boundary of the SC region in Fig. 1, and its conjugated pair $\theta=\pi/2$ exhibits no clear signature of SC possibly because of less frustration.

In summary, we perform large-scale DMRG simulations of the Moiré t-J model on four-leg cylinders at small hole doping. By varying the spin-dependent hopping phase induced by the out-of-plane electric field, we identify two conjugated SC phases, one of which is characterized by the coexistence of singlet $d$-wave, triplet $p$-wave SC and PDW, and the other ferromagnetic triplet SC. Our study supports twisted TMDs as a highly tunable platform for realizing exotic SC phases.

Data Availability.— The ITensor DMRG code and the data for all the figures in the main text and SM can be accessed by https://github.com/cfengno1/Moire-t-J-Model.

Acknowledgments.— We thank Yuchi He for useful comments. This work was supported by the U.S. Department of Energy, Office of Basic Energy Sciences under Grant No. DE-FG02-06ER46305. ITensor library Fishman et al. (2022) is used in this work for all DMRG calculations.

Supplemental Materials

In the Supplemental Materials, we provide additional results to support the claims made in the main text. In Sec. .1, we verify that the ground state satisfies $\sum_{i}\hat{S}_{i}^{z}=0$. In Sec. .2, we present the relative signs for pairings along different bonds, from which one can deduce the pairing symmetries. In Sec. .3, the spin structure factors for $\theta=\pi/12$ and $7\pi/12$ are given to show the underlying in-plane $120^{\circ}$ Néel and ferromagnetic orders respectively. In Sec. .4, the Green’s functions $G(r)$ are fitted by exponential decay. Sec. .5 shows the transformation between the STPSC and FMTSC phases. Sec. .6 shows the power-law decay of the pairing and CDW correlations on a width-3 cylinder. Sec. .7 shows the power-law fittings for the pairing correlations in the gray area (non-SC regime) of the phase diagram Fig. 1 in the main text, where large Luttinger exponents are found. Finally, in Sec. .8 the spin structure factors and different correlation functions of the XC4 cylinder are shown, exhibiting the same phases as those in the YC4 cylinder.