Mixed singlet-triplet superconducting state within the moir\'e $t$-$J$-$U$ model as applied to the description of twisted WSe$_2$ bilayer

TL;DR

This paper investigates a moiré $t$-$J$-$U$ model for twisted WSe$_2$ bilayer, revealing a mixed singlet-triplet superconducting state with two domes in the phase diagram and potential tunability via gate voltage.

Contribution

It introduces a detailed analysis of the mixed singlet-triplet superconducting state in twisted WSe$_2$ using the $t$-$J$-$U$ model and compares Gutzwiller and Hartree-Fock methods.

Findings

Reproduces two superconducting domes consistent with experiments.

Identifies a highly unconventional mixed $d+id$ and $p-ip$ pairing symmetry.

Suggests possible transition between pure singlet and triplet pairing via gate voltage.

Abstract

We analyze an analog of the $t$-$J$-$U$ model as applied to the description of a single moir\'e flat band of twisted WSe$_2$ bilayer. To take into account the correlation effects induced by a significant strength of the Coulomb repulsion, we use the Gutzwiller approach and compare it with the results obtained by the Hartree-Fock method. We discuss in detail the graduate appearance of a two dome structure of the superconducting state in the phase diagram by systematically increasing the Coulomb repulsion integral, $U$. The two superconducting domes residing on both sides of a Mott insulating state can be reproduced for a realistic parameter range in agreement with the available experimental data. According to our analysis the paired state has a highly unconventional character with a mixed $d+id$ (singlet) and $p-ip$ (triplet) symmetry. Both components of the mixed paired state are of comparable amplitudes. However, as shown here, a transition between pure singlet and pure triplet pairing should be possible in the considered system by tuning the gate voltage, which controls the magnitude of the valley-dependent spin-splitting in the system.

Tables

Table 1. Table 1: The six possible values of the symmetry factor, M 𝑀 M , and the corresponding values of the parity factor p 𝑝 p (second column), appearing in Eq. ( 12 ). The corresponding six symmetries of the superconducting gap together with their parities are provided in the third and fourth columns. The Cooper pair spin state compatible with a given symmetry is provided in the fifth column.

M	p	gap symmetry	parity	spin state
0	0	$extended$ $s$	even	singlet
1	1	$p_x+i{p}_y$	odd	triplet
2	0	$d_{x^2-y^2}+i{d}_{xy}$	even	singlet
3	1	$f$	odd	triplet
4	0	$d_{x^2-y^2}-i{d}_{xy}$	even	singlet
5	1	$p_x-i{p}_y$	odd	triplet

Equations

$$\hat{H}={\hat{H}}_0+{\hat{H}}_J+{\hat{H}}_U,$$

$${\hat{H}}_0=\sum _{<ij>\sigma }{t}_{ij\sigma }{\hat{c}}_{i\sigma }^{†}{\hat{c}}_{j\sigma },$$

$$t_{ij\sigma }=|t|e^{i\sigma {\varphi }_{ij}},$$

$$\begin{split}\hat{H}_{J}&=J\sideset{}{{}^{\prime}}{\sum}_{<ij>}\bigg{(}\hat{S}^{z}_{i}\hat{S}^{z}_{j}+\cos{(2\phi_{ij})}\sum_{\alpha=x,y}\hat{S}^{\alpha}_{i}\hat{S}^{\alpha}_{j}\\ &+\sin{(2\phi_{ij})}(\mathbf{\hat{S}}_{i}\times\mathbf{\hat{S}}_{j})\cdot\hat{z}\bigg{)},\end{split}$$

$${\hat{H}}_U=U\sum _i{\hat{n}}_{i\uparrow }{\hat{n}}_{i\downarrow },$$

$$|{\Psi }_G⟩=\hat{P}|{\Psi }_0⟩,$$

$$\hat{P}=\prod _i\sum _{\Gamma }{\lambda }_{i,\Gamma }|\Gamma {⟩}_{ii}⟨\Gamma |,$$

$$\begin{array}{rl}{S}_{ij}^{\sigma \bar{\sigma }}=⟨{\hat{c}}_{i\sigma }{\hat{c}}_{j\bar{\sigma }}& {⟩}_0,P_{ij\sigma }=⟨{\hat{c}}_{i\sigma }^{†}{\hat{c}}_{j\sigma }{⟩}_0,\\ & n_{i\sigma }=⟨{\hat{n}}_{i\sigma }{⟩}_0,\end{array}$$

$$E_G=\frac{⟨{\Psi }_G|\hat{H}|{\Psi }_G⟩}{N⟨{\Psi }_G|{\Psi }_G⟩}=\frac{1}{N}⟨\hat{H}{⟩}_G,$$

$$\begin{array}{rl}{\Delta }_{ij}^{\sigma \bar{\sigma }}& =⟨{\hat{c}}_{i\sigma }{\hat{c}}_{j\bar{\sigma }}{⟩}_G=q^2⟨{\hat{c}}_{i\sigma }{\hat{c}}_{j\bar{\sigma }}{⟩}_0,\\ {\tilde{P}}_{ij\sigma }& =⟨{\hat{c}}_{i\sigma }^{†}{\hat{c}}_{j\sigma }{⟩}_G=q^2⟨{\hat{c}}_{i\sigma }^{†}{\hat{c}}_{j\sigma }{⟩}_0,\end{array}$$

$${\Delta }_{\theta }^{\sigma \bar{\sigma }}=⟨{\hat{c}}_{i\sigma }{\hat{c}}_{j(\theta )\bar{\sigma }}{⟩}_G,$$

$${\Delta }_{M,p}^{\sigma \bar{\sigma }}=\frac{i^p}{6}\sum _{\theta }{e}^{-iM\theta }{\Delta }_{\theta }^{\sigma \bar{\sigma }},$$

$$\begin{array}{rl}{\Delta }_{M,p}^s& =({\Delta }_{M,p}^{\uparrow \downarrow }-{\Delta }_{M,p}^{\downarrow \uparrow })/\sqrt{2},\\ {\Delta }_{M,p}^t& =({\Delta }_{M,p}^{\uparrow \downarrow }+{\Delta }_{M,p}^{\downarrow \uparrow })/\sqrt{2}.\end{array}$$

$${\hat{P}}_i^2=1+x{\hat{d}}_i^{HF},$$

$$\begin{array}{rl}{\lambda }_d^2& =1+x(1-n_s{)}^2,\\ {\lambda }_s^2& =1-x{n}_s(1-n_s{)}^2,\\ {\lambda }_{\varnothing }^2& =1+x{n}_s^2,\end{array}$$

$$⟨{\hat{H}}_0{⟩}_G=\sum _{ij}{q}^2{t}_{ij\sigma }⟨{\hat{c}}_{i\sigma }^{†}{\hat{c}}_{j\sigma }{⟩}_0,$$

$$\begin{split}\langle\hat{H}_{J}\rangle_{G}&=\lambda^{4}_{s}\;J\sum_{<ij>}\bigg{(}\frac{1}{2}\;\sum_{\sigma}\;e^{i\sigma 2\phi_{ij}}\langle\hat{c}^{\dagger}_{i\sigma}\hat{c}_{i\bar{\sigma}}\hat{c}^{\dagger}_{j\bar{\sigma}}\hat{c}_{j\sigma}\rangle_{0}\\ &+\frac{1}{4}\;\sum_{\sigma\sigma^{\prime}}\sigma\sigma^{\prime}\langle\hat{n}^{HF}_{i\sigma}\hat{n}^{HF}_{j\sigma^{\prime}}\rangle_{0}\bigg{)},\end{split}$$

$$⟨{\hat{H}}_U⟩=N{\lambda }_d^{2}U{n}_s^2,$$

$$\begin{split}\hat{\mathcal{H}}_{\textrm{eff}}&=\sum_{ij\sigma}t^{\textrm{eff}}_{ij\sigma}\hat{c}^{\dagger}_{i\sigma}\hat{c}_{j\sigma}-\mu^{\textrm{eff}}\sum_{i\sigma}\hat{n}_{i\sigma}\\ &+\sideset{}{{}^{\prime}}{\sum}_{<ij>\sigma}\big{(}\Delta^{\textrm{eff}}_{ij\sigma\bar{\sigma}}\hat{c}_{j\sigma}\hat{c}_{i\bar{\sigma}}+H.c.\big{)},\end{split}$$

$$t_{ij\sigma }^{\text{eff}}\equiv \frac{\partial \mathcal{F}}{\partial P_{ij\sigma }},({\Delta }_{ij\sigma \bar{\sigma }}^{\text{eff}}{)}^{\ast }\equiv \frac{\partial \mathcal{F}}{\partial S_{ji}^{\sigma \bar{\sigma }}},{\mu }^{\mathrm{eff}}\equiv -\frac{\partial \mathcal{F}}{\partial n_s},$$

$${\Delta }^{\sigma {\sigma }^{\prime }}(\mathbf{k})=\sum _{\theta }{\Delta }_{\theta }^{\sigma {\sigma }^{\prime }}{e}^{i\mathbf{k}{\mathbf{R}}_{ij(\theta )}},$$

$$\begin{split}\Delta^{s}(\mathbf{k})&=\big{(}\Delta^{\uparrow\downarrow}(\mathbf{k})-\Delta^{\downarrow\uparrow}(\mathbf{k})\big{)}/\sqrt{2}\\ \Delta^{t}(\mathbf{k})&=\big{(}\Delta^{\uparrow\downarrow}(\mathbf{k})+\Delta^{\downarrow\uparrow}(\mathbf{k})\big{)}/\sqrt{2}\end{split}$$

$${\hat{H}}_J=J\sum _{ij}{\hat{\mathbf{S}}}_i{\hat{\mathbf{S}}}_j.$$

Full text

I Introduction

In recent years the twisted moiré structures have gathered a significant amount of interest due to their rich physics on one side and high degree of tunability on the other. The development in the field has been mostly triggered by the discovery of both correlated insulating state and unconventional superconductivity in the magic angle twisted bilayer graphene (TBG)Cao et al. (2018a, b). Since than, many other moiré structures have been intensively studied which are bilayer boron nitrideXian et al. (2019), multilayer grapheneZhang et al. (2021); Chen et al. (2019); Shen et al. (2020) as well as both homo- and hetero-bilayer transition metal dichalcogenides (TMD)Wang et al. (2020); Liao et al. (2020); Ulstrup et al. (2020); Ghiotto et al. (2021).

In moiré systems, for certain twist angles between the atomic layers, flat energy bands are created. This feature points to a significant role of inter-electronic interactions which is a situation similar to that known from the so-called strongly correlated electron systems such as copper based high-temperature superconductorsSpałek et al. (2022). In fact, both the well known cuprates and the novel moiré systems exhibit similar feature in the phase diagram which is the appearance of superconducting (SC) domes residing in close proximity of the correlated insulating stateCao et al. (2018b); Wang et al. (2020).

For the case of TBG the flat bands and the resulting correlation induced phenomena are fragile when it comes to deviations of the twist angle around the so-called magic values. On the other hand, in TMD moiré systems the effect is more robust and the flat bands are observed for larger spectrum of anglesWang et al. (2020); An et al. (2020); Ghiotto et al. (2021). What also distinguisches the TMD systems from those graphene based is the appearance of significant spin-orbit coupling.

As it has been shown both by the DFT calculations and continuum model based approach, the twisted homo-bilayer WSe2 (tWSe2) can be described by a single-band tight-binding model on a triangular lattice with spin-orbit coupled hoppings which reproduce the flat electronic band with a valley-dependant spin-splittingWang et al. (2020); Pan et al. (2020). In such an approach the so-called displacement field controlled by the bias voltage across the bilayer, tunes the spin-dependent phase of the hopping amplitudes. This allows for changing the magnitude of the spin-splitting in the system. In order to take into account the electron-electron interactions, such approach can be supplemented with the onsite Coulomb repulsion term leading to the Hubbard model. It has been estimated that the value of the Coulomb integral $U$ places tWSe2 in the moderate correlations regimeWang et al. (2020). However, the magnitude of electronic correlations in those systems can be controlled by the twist angle which changes the $U/W$ parameter, with $W$ being the bare band width. The Hubbard model has been studied in the context of both magnetically ordered states and $d+id$ spin-singlet supercondcutivity of tWSe2Zang et al. (2021); Bélanger et al. (2022); Zang et al. (2022); Pan et al. (2020); Pan and Das Sarma (2022). In such picture the SC pairing comes as a result of the electronic correlations which become significant for large enough value of the Hubbard $U$. Another canonical model which can be analyzed in the context of electron-electron correlation effects in TMD moiré bands is the $t$-$J$ modelWu et al. (2018); and Sheng (2022); Chen and Sheng (2023) within which the Cooper pairing comes naturally from the exchange interaction term $\sim J$. However, it should be noted that the situation in this context is not completely resolved so far and a well established description of the paired state in those systems has not been reached.

Here, we propose an analog of the $t$-$J$-$U$ model as applied to the description of a single moiré flat band of tWSe2. The $t$-$J$-$U$ hamiltonian has been initially proposed for the description of copper based high-TC superconductorsDaul et al. (2000); Basu et al. (2001); Xiang et al. (2009) and as shown by us in previous years it leads to good quantitative agreement between theory and experiment for selected fundamental features of the paired stateSpałek et al. (2017). This fact together with similarities between the cuprates and twisted moiré van der Waals structures constitutes the main motivation for the analysis of unconventional superconducting phase in tWSe2 with the use of the $t$-$J$-$U$ based approach. Due to the significant magnitude of the onsite Coulomb repulsion in such an analysis one has to take into account the electronic correlation effects beyond the level of a simple Hartree-Fock method (HF). Here, we use the Gutzwiller approach (GA), which is efficient and catches the basic effects related with electron-electron interactions. As we show the two superconducting domes residing on both sides of the insulating state can be reproduced with the use of the presented method in agreement with the experimental indicatorsWang et al. (2020). The character of the obtained superconducting state is highly unconventional with a mixed singlet-triplet Cooper pairing. Interestingly, according to our analysis the balance between the singlet and triplet contributions to the pairing could be controlled by the gate voltage in a real experimental setup. In particular, in proper gate voltage range, a transition between pure singlet and pure triplet pairing can be achieved within a single non-hybrid system.

II Model and method

Here we apply a moiré band $t$-$J$-$U$ hamiltonian of the form

[TABLE]

where $\hat{H}_{0}$ represents a tight-binding model with spin-orbit coupled hoppings on a triangular lattice. As shown in Refs. Wang et al., 2020; Pan et al., 2020, such single particle hamiltonian can well describe the top moiré band of tWSe2. Namely,

[TABLE]

where $t_{ij\sigma}$ are the hopping parameters between the nearest-neighboring sites of the form,

[TABLE]

where $\phi_{ij}=-\phi_{ji}$ due to hermiticity condition and $\sigma=+1,-1$ represents the spin up and down states associated with $+K$ and $-K$ valleys, respectively. The hopping parameters fulfill the $C_{3}$ rotational symmetry and reconstruct the spin-splitted Fermi surfaces. This is a manifestation of the fact that the system is non-centrosymmetric and spin-orbit coupling appears. In order to simplify the situation, here we limit to the nearest-neighbor hopping terms (summation over $<ij>$) what is justified by the fact that the longer-ranged hoppings are at least one order of magnitude smallerWang et al. (2020); Pan et al. (2020). In the further analysis we take the values of hoppings provided in Ref. Wang et al., 2020 for the twist angle $5.08^{\circ}$ which is close to the one for which the signatures of superconductivity have been identified experimentally. As shown in Fig. 1, by changing the value of the so-called displacement field one can tune both $|t|$ and $\phi$, where the latter determines the magnitude of the valley-dependent spin-splitting. In real experimental setup the displacement field is controlled by the bias voltage across the bilayer.

The intersite exchange interaction term appearing in Hamiltonian (1) has the following form

[TABLE]

where $\hat{S}^{\alpha}_{i}$ for $\alpha=x,y,z$ are the three components of the electron spin operator and $\phi_{ij}$ is the phase of the spin up hopping between $i$ and $j$ lattice sites. The summation above runs over the nearest neighbor sites only. The prime indicates that each bond between $i$ and $j$ is taken only once. Due to the appearance of the complex spin-dependant hoppings, the exchange term is slightly modified with respect to the one originally appearing within the $t$-$J$-$U$ model as applied for the cupratesSpałek et al. (2017). In particular, here the Dzyaloshinskii-Moriya (DM) term appears which, however, might be significantly suppressed for very low values of the displacement fields for which $\phi\approx 0$. In such a case the exchange interaction takes the standard form, $\sum_{ij}\mathbf{\hat{S}}_{i}\mathbf{\hat{S}}_{j}$.

The last term in Hamiltonian (1) corresponds to the onsite Coulomb interaction

[TABLE]

where $\hat{n}_{i\uparrow}$ ($\hat{n}_{i\downarrow}$) are the number of particle operators on given site with spin up (down), while $U$ corresponds to the Coulomb interaction strength.

In order to take into account the correlation effects induced by a significant value of the Hubbard $U$ we apply the zeroth order diagrammatic expansion of the Gutzwiller wave function (DE-GWF)Spałek et al. (2022) which is equivalent to the $Statistically$ $consistant$ $Gutzwiller$ $Approximation$ (SGA)Abram et al. (2013). As we show here, for the considered problem the zeroth order expansion provides a good enough approximation to reach a qualitative agreement with fundamental experimental observations for tWSe2. Within this method the many-body effects are captured by taking the variational correlated wave function in the following form

[TABLE]

where $|\Psi_{0}\rangle$ is the non-correlated (mean-field) state to be defined below and $\hat{P}$ is the Gutzwiller-like correlation operator of the following form

[TABLE]

where $|\Gamma\rangle$ are the four states from the local basis $|\emptyset\rangle_{i}$, $|\uparrow\rangle_{i}$, $|\downarrow\rangle_{i}$, $|\uparrow\downarrow\rangle_{i}$ and $\lambda_{i,\Gamma}$ are the corresponding variational parameters. The $|\Psi\rangle_{0}$ state is taken here as the uncorrelated superconducting state which is defined by the pairing and hopping mean-field parameters, as well as the electron occupation expectation value

[TABLE]

respectively. In the above $\langle\hat{o}\rangle_{0}$ stands for the expectation value of the $\hat{o}$ operator in the state $|\psi\rangle_{0}$, and $\bar{\sigma}$ represents spin orientation opposite to $\sigma$.

Here we focus on the unconventional features of the paired state, therefore, for simplicity we assume no magnetic or charge ordering. Hence, we take $\lambda_{i,\Gamma}\equiv\lambda_{\Gamma}$, $\lambda_{\uparrow}=\lambda_{\downarrow}=\lambda_{s}$, $n_{i\uparrow}=n_{i\downarrow}\equiv n_{s}$. Due to the appearance of the spin-orbit coupling in the model both singlet and triplet pairings are possible which means that we need to treat $S_{ij}^{\uparrow\downarrow}$ and $S_{ij}^{\downarrow\uparrow}$ separately allowing both for pure spin-singlet pairing ($S_{ij}^{\uparrow\downarrow}=-S_{ij}^{\downarrow\uparrow}$), pure spin-triplet pairing ($S_{ij}^{\uparrow\downarrow}=S_{ij}^{\downarrow\uparrow}$) as well as their mixture. As can be shown the equal spin state (ESP) with $S_{ij}^{\uparrow\uparrow}=S_{ij}^{\downarrow\downarrow}$ and $S_{ij}^{\sigma\bar{\sigma}}=0$ is equivalent to the $m_{s}=0$ spin triplet state with $S_{ij}^{\uparrow\downarrow}=S_{ij}^{\downarrow\uparrow}$ and $S_{ij}^{\sigma\sigma}=0$ and as such is taken into account in our analysisSmidman et al. (2017).

The main task within our approach is to calculate the ground-state energy per lattice site of the system in the correlated state within the zeroth order diagrammatic expansion

[TABLE]

where $N$ is the number of moiré lattice sites. As we show in more detail in Appendix A, the ground state energy in the correlated state can be expressed as a function of the variational parameters, $\lambda_{\Gamma}$, the mean-field parameters, $S_{ij}^{\sigma\bar{\sigma}}$ and $P_{ij\sigma}$, as well as $n_{s}$. In order to carry out the minimization procedure effectively with the self-consistency condition fulfilled we apply the Effective Hamiltonian SchemeKaczmarczyk (2015), which we show in some more detail in Appendix B as applied to the considered situation.

After determining the mean-field and the variational parameters one can calculate the correlated superconducting gap amplitude as well as the hopping expectation values in the correlated state

[TABLE]

where $q=\lambda_{s}(\lambda_{d}\;n_{s}+\lambda_{\emptyset}(1-n_{s}))$. Since both the hopping and the exchange interactions are taken into account to the nearest neighbors only, we also calculate the six nearest neighbor real space pairing amplitudes both for $\uparrow\downarrow$ and $\downarrow\uparrow$ spin configurations. We use the following notation

[TABLE]

where $j(\theta)$ is the nearest neighboring lattice site of site $i$, while $\theta$ corresponds to the angle between the positive half-$x$ axis and the $\mathbf{R}_{ij}=\mathbf{R}_{i}-\mathbf{R}_{j}$ vector. The angles $\theta\in\{0,\;\pi/3,\;2\pi/3,\;\pi,\;4\pi/3,\;5\pi/3\}$ correspond to six subsequent nearest neighbor lattice sites, $j$. Due to the homogeneity of the system, the $\Delta^{\sigma\bar{\sigma}}_{\theta}$ gap amplitude does not depend on the position of the $i$ lattice site. In order to identify the symmetry of the obtained paired state for a given spin configuration we transform the six nearest neighbor pairing amplitudes into a six symmetry resolved pairing amplitudes of the following form

[TABLE]

where the summation runs over the angles $\theta\in\{0,\;\pi/3,\;2\pi/3,\;\pi,\;4\pi/3,\;5\pi/3\}$, and $M$ is the symmetry factor, which takes integer values and corresponds to six possible pairing symmetries (cf. Table 1). The $p$ parameter corresponds to the parity of the particular symmetry. Namely, for even parity symmetries we have $p=0$ and for odd parities we have $p=1$.

Finally, the correlated gap amplitudes corresponding to subsequent symmetries can be transformed to those corresponding to singlet and triplet spin states of the Cooper pairs

[TABLE]

Only three of the considered symmetries are compatible with singlet ($d+id$-, $d-id$-, and $extended$ $s$-$wave$) which means that they always give zero for the triplet component ($\Delta^{t}_{M,p}=0$ for $M=0,2,4$). The remaining $p+ip$-, $p-ip$-, and $f$-$wave$ symmetries are compatible with the spin-triplet pairing, therefore, $\Delta^{s}_{M,p}=0$ for $M=1,3,5$. In our calculation scheme we determine all the twelve nearest neighbor gap amplitudes, $\Delta^{\uparrow\downarrow}_{\theta}$, $\Delta^{\downarrow\uparrow}_{\theta}$, in the correlated state $|\Psi_{G}\rangle$ without any constraints which would limit us only to one selected symmetry among those provided in Table 1. Next, we transform them to the final symmetry resolved singlet and triplet gap amplitudes by using Eqs. (12), (LABEL:eq:gap_symmetries_spin). In such a manner we can encompass both the pure superconducting symmetries as well as their mixtures. For the sake of clarity, in the following Sections, we use the symmetry names ($p\pm ip$, $d\pm id$, $f$ etc.) in the subscripts of the symmetry resolved superconducting gaps instead of the values of the $M$ and $p$ factors.

It should be noted that within the so-called diagrammatic expansion method (DE-GWF) one starts from the renormalized mean-field theory as the zeroth-order result and proceed with incorporating systematically the nonlocal correlations of increased range in higher orders. Therefore, within the higher order calculations the correlated gaps and hopping expectation values are no longer related with their non-correlated counterparts by a simple renormalization given by Eqs. (10) (cf. Ref. Spałek et al., 2017). Also, in such case the hopping and pairing between longer-range neighboring sites needs to be included which makes the analysis much more involved. This allows to take into account the correlation effects with better accuracy. Nevertheless, previous calculations carried out within the SGA approach (equivalent to the zeroth order DE-GWF) and the higher order DE-GWF method, both has lead to the stability of the $d$-$wave$ paired state, and the dome like behavior of the SC gap for the case of the single band $t$-$J$-$U$ model on the square-latticeAbram et al. (2013); Spałek et al. (2017). Therefore, one would expect that also here the inclusion of the higher order terms would not change the main result of the analysis, which is related to the symmetry of the gap and the general form of the phase diagram.

The code which was written to carry out the numerical calculations as well as the data behind the figures presented in the following Sections are available in the open repositoryZegrodnik (2023).

III Results

In the following we analyze the features of the unconventional superconducting state for the WSe2 bilayer with the twist angle $5.08^{\circ}$ and the displacement field $D=0.45$ V/nm for which the signatures of both the correlated insulating state and superconductivity have been identified experimentallyWang et al. (2020). The bare band width for the corresponding hopping parameters is $W\approx 90$ meV. It has been estimated the onsite Coulomb repulsion integral, $U$, is comparable to $W$ for the considered systemWang et al. (2020); Pan et al. (2020). Unless stated otherwise, we take $U=120$ meV which places the system close to the moderately correlated regime. Within the present approach we take the value of $J$ as determined from the formula $J=4|t|^{2}/U=2.54$ meV corresponding to the mapping of the Hubbard model to the spin Heisenberg model shown in Ref. Pan et al., 2020.

As one can see in Fig. 2, for the chosen model parameters the symmetry resolved gap amplitudes create two superconducting domes residing on both sides of the half-filled situation ($n=1$). Here, we define the band filling as the number of electrons in the narrow band per single moiré lattice site, $n=2n_{s}$. According to our calculations the paired state has a mixed singlet-triplet character with comparable magnitudes of both components of the pairing. The singlet (triplet) pairing component is realized by the $d+id$ ($p-ip$) symmetry of the superconducting gap. Apart from the two gap amplitudes shown in the Figure, all the remaining ones listed in Table 1 have turned out to be zero for the given range of model parameters. One should note that pairing in the mixed $d/p$ channel has been proposed to appear due to density waves fluctuations very recently according to functional renormalization group calculationsKlebl et al. (2022). Also, the singlet/triplet pairing has been analyzed within the $t$-$J$ model showing the dominant singlet $d$-wave contributionChen and Sheng (2023). The singlet-triplet mixing obtained here is a consequence of spin-orbit coupling, due to which parity is no longer a good quantum numberSmidman et al. (2017). Among the possible pairing symmetries the $d\pm id$ and $p\pm ip$ optimize the condensation energy due to the fact that they provide a fully gapped situation in contradiction to the $extended$ $s$-$wave$ and $f$-$wave$ symmetries, for which the line nodes would cross the Fermi surface, in the considered doping range (cf. Appendix C).

In Fig. 2 (b) we show that the nearest neighbor hopping expectation value $\tilde{P}$ is being significantly suppressed as one approaches the half-filled situation ($n=1$) from both sides, which indicates the appearance of the correlation induced insulating state (Mott insulator). The insulating behavior suppresses superconductivity in the close proximity of $n=1$ creating the double-dome structure seen in experimentsWang et al. (2020). The resulting two SC domes are not symmetric with respect to half-filling, which is caused by the appearance of the van Hove singularity for $n\approx 0.8$, which in turn enhances the pairing at the hole-doped side of the phase diagram [cf. Fig. 2 (c)]. One should note that the approach based on the onsite Gutzwiller-like correlator results in a trivial character of the Mott transition. For the complete description of this effect one should take into account more involved Jastrov and/or doublon-holon correlatorsCapello et al. (2005); Miyagawa and Yokoyama (2011).

In Fig. 3 we show the subsequent contributions to the total system energy as well as their sum, all as functions of band filling as obtained from the Gutzwiller approach and the Hartree-Fock method. The former takes into account the correlation effects induced by the Hubbard $U$ while the latter neglects them. As one can see the two methods deviate significantly close to half-filling where the correlations are most pronounced. Again, within the Gutzwiller approach one can see the signature of the insulating behaviour for $n=1$ where the kinetic energy approaches zero, together with the Coulomb energy as both the electron hopping and the double occupancies are being suppressed. It should be noted that the appearance of the insulating state is concomitant with a significant decrease of the exchange interaction energy as shown in Fig. 3 (b). This effect is caused by the fact that the $\lambda_{s}$ parameter, over which the expectation value of the $\sim J$ term is premultiplied [cf. Eq. (17)], increases close to half filling as shown in Fig. 6 (a) in Appendix A. One may understand this as an electron correlation induced increase of the renormalized exchange interaction integral $J^{\prime}=\lambda_{s}^{4}\;J$. As shown in Fig. 3 (d), the Gutzwiller method gives total energy lower than HF as it should be since the former is of better quality.

The appearance of the mixed singlet-triplet pairing seen in Fig. 2 is a consequence of the spin-orbit coupled hoppings and the resulting spin-splitting of the Fermi surface. It should be noted that with the use of the displacement field one can tune the phase of the hoppings, $\phi$ [cf. Fig. 1 (b)], which in turn determines the magnitude of the spin splitting as well as the strength of the DM interaction [cf Eq. 4]. The influence of the DM interaction on the paired state is briefly discussed in Appendix D. Since the spin-splitting plays an important role when it comes to the appearance of the mixed spin-triplet pairing, it is tempting to check if the balance between the singlet and triplet contributions can be modified by changing $\phi$. In Fig. 4 we show the $\phi$-dependance of the two SC gaps, $\Delta^{s}_{d+id}$ (singlet) and $\Delta^{t}_{p-ip}$ (triplet), for two selected band fillings. For the sake of simplicity in these calculations we keep $|t|$ constant and change only the phase of the hopping. As one can see by changing $\phi=\pi\rightarrow\pi/2\rightarrow 0$, we also carry out the transition between singlet$\rightarrow$triplet$\rightarrow$singlet pairing, respectively. Purely real (imaginary) hopping values result in purely singlet (triplet) pairing. In between those three distinguished situations a mixed singlet-triplet state appears. The presented data indicate the possibility of tuning the balance between the singlet and triplet contributions by changing the displacement field. The important question is if it is possible to tune the gate voltage in a real experimental setup so as to reach the displacement field values for which $\phi=\pi/2$ in order to generate pure spin-triplet pairing. This would allow for switching between singlet and triplet pairing in an $in$ $situ$ manner by simply adjusting the top and bottom gate voltage.

As the last part of our analysis we consider the calculated phase diagram in the $(U,n)$ plane. It can be seen from Fig. 5, that large enough $U$ is crucial in order to induce the paired state in the system. As the $U$ value is increased, first, a weak superconducting state appears in the range of band fillings which correspond to the van Hove singularity ($n\approx 0.8$) where the high values of density of states create favorable conditions for the pairing. When $U$ approaches the value of $W\approx 90$ meV, the pairing is strengthen around half-filling. This effect is caused by the exchange interaction enhancement which is due to $\lambda_{s}>1$, as discussed earlier (cf. Eq. 17). For $U\gtrsim W$, the SC state stability range becomes divided and two SC domes are created which is a direct result of the insulating state appearance at $n=1$. It should be noted that even though the electronic structure of the model and the features of the superconducting state obtained here are significantly different from those corresponding the copper based superconductors there are some clear similarities between the two systems. Namely, in both systems the two superconducting domes divided by a correlated insulating state appear. With this respect, the structure of the $(U,n)$-phase diagram calculated here resembles the one obtained by us previously for the case of a single-band model as applied for the description of the cupratesSpałek et al. (2017); Zegrodnik and Spałek (2017). However, it should be noted that also other approaches have been reported in recent years which represent different methods and points of view to the problem of pairing in single band models of copper based superconductorsGunnarsson et al. (2015); Qin et al. (2020); Dong et al. (2022a, b). In particular, the approach presented here does not take into account the dynamical screening effects which leads to frequency dependent $U$ parameter and may affect the strength of the pairing as well as the resulting critical temperatureNilsson et al. (2019).

IV Conclusions

We have applied the moiré $t$-$J$-$U$ model to the description of unconventional superconducting state in tWSe2. In order to take into account the electron-electron correlations we have used the Gutzwiller method. As we show our approach leads to the reconstruction of two superconducting domes residing on both sides of the correlation induced insulating state appearing at half-filling, which is in agreement with the experimental data shown in Ref. Wang et al., 2020. A similar structure of the phase diagram as the one obtained here has also been discussed by us recently in a quantum-dot square lattice system with the use of the variational Monte Carlo methodBiborski et al. (2021). According to our analysis the suppression of the electron hopping processes near the half-filled situation is concomitant with the enhanced role of the exchange interaction term as the corresponding renormalization factor becomes larger than unity, $\lambda_{s}>1$ (cf. Fig. 6 and Eq. 17). This effect also enhances the superconducting phase around the $n=1$ case, however, for large enough value of the Coulomb repulsion integral, $U\gtrsim W$, the appearance of the Mott insulating behaviour separates the SC state stability regime into two superconduting domes (cf. Fig. 5). It is worth mentioning that it was proposed very recently within a two-orbital model that the Mott-Hubbard behavior in tWSe2 can be switched into Hund correlated state by tuning the twist angleRyee and Wehling (2023). Such effect is not taken into account here as it requires a two-band extension of the considered model.

The character of the paired state obtained here is highly unconventional with a mixed spin singlet-spin triplet pairing and non-zero values of the the corresponding $d+id$ and $p-ip$ components of the gap symmetry, respectively. As we show the balance between the singlet and triplet contributions can be tuned by adjusting the gate voltage which determines the valley-dependent spin splitting of the electronic structure (cf. Fig. 4). In particular, if the hopping phase $\phi=\pi/2$ can be reached experimentally, the transition between pure spin-singlet and pure spin-triplet pairing can be realized in the system.

It should be noted that various magnetic and charge ordered states have been analyzed by a number of methods within the Hubbard and Heisenberg models as applied to the description of tWSe2Zang et al. (2021); Bélanger et al. (2022); Zang et al. (2022); Pan et al. (2020); Pan and Das Sarma (2022); Biderang et al. (2022); Kiese et al. (2022). Furthermore, the pair-density-wave state has been shown to appear in weak out-of-plane Zeeman fieldWu et al. (2022). It is expected that some of those symmetry-broken states may also appear in the present $t$-$J$-$U$ based approach. In particular, in analogy with the cuprates one would expect that close to half-filling antiferromagnetic phase should become stable. As argued in Refs. Pan et al., 2020; Zang et al., 2021 the DM interaction acts as an anisotropy that favors in-plane spin ordering and should select the subset of $120^{\circ}$ antiferromagnetic states. In the considered model the DM interaction is affected by the hopping phase [cf. Eq. 4], which in turn can be tuned by the displacement field. Therefore, for particular range of the displacement fields also a tendency towards ferromagnetism can appear. It would be interesting to analyze the interplay between the magnetic/charge ordering and the mixed superconducting pairing considered here. In particular, the appearance of charge density waves may lead to modulation of the pairing amplitudes in space as has been shown by us in recent years for the analogical approach as applied to cupratesZegrodnik and Spałek (2018). Such study requires much more involved calculations and is beyond the scope of the present analysis. We should be see progress along this line soon.

Finally, the spin-triplet $p\pm ip$ paired state is interesting when it comes to its topological features. Namely, if a gate-defined nanowire structure based on tWSe2 could be realized, the $p\pm ip$ pairing could give rise to Majorana bound states at the edges of such nanowire. What is very attractive in such concept is that all the ingredients required for such effect can appear in a single material. This is in opposition to other considered systems in which subsequent factors like spin-orbit coupling or superconducting pairing are introduced with the use of proximity effect in a hybrid systemOreg et al. (2010); Lutchyn et al. (2010).

V Acknowledgement

Discussions with P. Wójcik are gratefully acknowledged. This research was founded by National Science Centre, Poland (NCN) according to decision 2021/42/E/ST3/00128. For the purpose of Open Access, the author has applied a CC-BY public copyright licence to any Author Accepted Manuscript (AAM) version arising from this submission.

Appendix A Expectation value of the $t$-$J$-$U$ Hamiltonian in the correlated state

Here we show how to express the expectation value of our $t$-$J$-$U$ Hamiltonian in the variational correlated state defined by Eq. (6) within the zeroth order diagrammatic expansion (DE-GWF method). As shown in Ref. Bünemann et al., 2012, it is convenient to impose the following condition on the correlation operator

[TABLE]

where $\hat{d}^{HF}_{i}=\hat{n}^{HF}_{i\uparrow}\hat{n}^{HF}_{i\downarrow}$, $\hat{n}^{HF}_{i\sigma}=\hat{n}_{i\sigma}-n_{s}$, and $x$ is yet another variational parameter. One can express the variational parameters $\lambda_{\Gamma}$ as a function of $x$ in the following manner

[TABLE]

therefore, in fact we are left with only one variational parameter. In such case the expectation values of subsequent terms of our Hamiltonian take the form

[TABLE]

[TABLE]

[TABLE]

where q=\lambda_{s}\big{(}\lambda_{d}n_{s}+\lambda_{\emptyset}(1-n_{s})\big{)}. Since all the four-operator expectation values which appear on the right-hand side of Eq. (17) are in the non-correlated state, one can carry out the standard Wick’s theorem decomposition in order to express them in terms of the mean-field parameters $S_{ij}^{\sigma\bar{\sigma}}$, $P_{ij\sigma}$. For the sake of clarity in Fig. 6 we provide the band filling dependence of the renormalization factor $q$ and the variational parameters $\lambda_{s}$, $\lambda_{d}$, and $\lambda_{\emptyset}$, appearing in expressions (16-18). It should be noted that the Hartree-Fock calculations are equaivalent to the situation for which $\lambda_{\Gamma}\equiv 1$, $q\equiv 1$.

Appendix B The effective Hamiltonian scheme

From the minimization condition of the ground state energy (9) one can derive the effective Hamiltonian, which for the case of superconducting phase has the form

[TABLE]

where the effective hopping, effective chemical potential, and the effective superconducting gap parameters are defined through the corresponding relations

[TABLE]

where $\mathcal{F}=E_{G}-\mu_{G}n_{s}$ with $\mu_{G}$ being the chemical potential determined in the correlated state. By using Eqs. (16-18) one can express the ground state energy, $E_{G}$ in the correlated state with the use of the mean-field parameters, variational parameters, as well as the band filling. This allows to calculate the values of the effective parameters from Eq. (20). Next, the effective Hamiltonian can be diagonalyzed by introducing the generalized Bogolubov de Gennes transformation and the self-consistent equations for $P_{ij\sigma}$, $S^{\sigma\bar{\sigma}}_{ij}$, $\mu_{G}$ can be derived in a standard manner. Additionally, the procedure of solving the self-consistent equations has to be supplemented with the minimization of the ground state energy with respect to the variational parameter $x$.

Appendix C Gap symmetry in the $\mathbf{k}$-space

Here we show explicitly the $\mathbf{k}$-dependency of the superconducting gap for selected pairing symmetries. It should be noted that in order to obtain the $\mathbf{k}$-dependant superconducting gaps one can transform the real-space pairing amplitudes given by equations (11) in the following manner

[TABLE]

where $j(\theta)$ is the nearest neighboring lattice site of site $i$, while $\theta$ corresponds to the angle between the positive half-$x$ axis and the $\mathbf{R}_{ij}=\mathbf{R}_{i}-\mathbf{R}_{j}$ vector. The summation is over the angles $\theta\in\{0,\;\pi/3,\;2\pi/3,\;\pi,\;4\pi/3,\;5\pi/3\}$, which correspond to six subsequent nearest neighbor lattice sites, $j$. In analogy to the real space gap amplitudes [cf. Eq. (LABEL:eq:gap_symmetries_spin)], also here we can define the singlet and triplet superconducting gaps in $\mathbf{k}$-space.

[TABLE]

In Fig. 7 we show the absolute value of the superconducting gaps for spin-singlet symmetries ($d\pm id$-$wave$ and $extended$ $s$-$wave$) and two spin-triplet symmetries ($p\pm ip$-$wave$ and $f$-$wave$), all as functions of momentum, $\mathbf{k}$. Additionally in Fig. 7 (f) we show the Fermi surface for one exemplary band filling value, $n=0.9$, for which the mixture of $d+id$ and $p-ip$ pairing has been shown to be stable in the considered model (cf. Fig. 2).

As one can see from the figures, the $d\pm id$ and $p\pm ip$ pairing symmetries provide a fully gapped situation at the Fermi surface in contradiction to the $extended$ $s$-$wave$ and $f$-$wave$ symmetries. This leads to the optimization of the condensation energy for the case of the $d\pm id$ and $p\pm ip$ mixing.

Appendix D The effect of DM interaction on the paired state

For the sake of completeness, in Fig. 8 we show the results obtained for the case when the spin-splitting is still present in the electronic structure, however, the exchange interaction term is taken in the standard form without the DM term. Namely,

[TABLE]

As one can see in such case the triplet component is still present, however it is significantly smaller with respect to the results presented in Fig. 2 where the DM term has been taken into account.