Introduction

The pairing of electrons in a superconductor is among the most intriguing effects in the study of collective phenomena. One strong driving force in the exploration of superconductors is the quest for ever higher critical temperatures1,2 close to ambient conditions with numerous high impact technological applications. So far the highest critical temperatures were either achieved in conventional, BCS-like superconductors but only under the application of immense pressure (with critical temperatures reaching room temperature) or at ambient pressure via an unconventional superconducting pairing mechanism3 (with critical temperatures reaching approximately 100 K). Developing a deepened understanding of the latter, raising the hope to relief the high pressure condition, is thus of utmost importance. At the same time, combining superconductivity with non-trivial topology is a promising route for quantum information sciences as such topological superconductors may harbor robust edge states at domain boundaries with topological properties advantageous to computing applications4.

However, realizing and controlling topological superconductors proves difficult to this date, with only a few candidate materials currently being suggested, e.g.,5,6,7,8,9. A new direction in the study of superconductivity opened up recently in twisted moiré quantum materials, i.e. two-dimensional van der Waals materials being stacked at a relative twist angle10,11,12,13,14,15,16,17,18,19,20,21,22,23,24,25,26,27,28. In these systems kinetic energy scales can be tuned by the twist angle allowing to promote the relative relevance of potential, spin-orbit coupling or other energy scales17. Indeed, topological properties as well as superconductivity were already demonstrated in these highly versatile systems and as a consequence they could provide an excellent opportunity to engineer topological superconductors. In particular, in moiré transition metal dichalcogenides, strong spin-orbit coupling or excitonic physics were experimentally explored. The observation of a Mott insulating state29,30 and other fascinating collective phenomena such as generalized Wigner crystals30, stripe phases31 and quantum anomalous Hall insulators32 confirmed the relevance of many-body interactions, and demonstrated the importance of their extended range. Moreover, near the insulating state, indications of superconductivity were found in homobilayer twisted WSe233.

Here, we explore this idea for moiré transition metal dichalcogenides (see Fig. 1a)17,29,33,34,35,36,37,38,39,40,41,42 by analyzing the Fermi surface instabilities of twisted hetero-bilayers of WX2/MoX\({}_{2}^{\prime}\) (X,X’ = S,Se) away from half filling of the moiré band. We unveil an exotic superconducting state near Van Hove filling described by form factors with eight zero crossings, arising from the extended range of interactions in these materials. We show that the superconducting ground state is formed by a chiral configuration, which is characterized by a full gap on the Fermi surface and non-trivial topology with winding number \(| {{{\mathcal{N}}}}| =4\). We argue that this type of topological superconductivity leads to distinct experimental signatures in quantum Hall transport measurements and elevates twisted hetero-bilayers of TMDs to prime candidates for experimental scrutiny of topological superconductivity.

Fig. 1: Correlated phase diagram of hetero-bilayer TMDs.
figure 1

a Sketch of a twisted WSe2/MoS2 bilayer with small twist angle θ. The resulting effective moiré potential34,39 is indicated by the contour shading at the top. The four red lines represent the chiral edge modes of the topological superconducting state with \({{{\mathcal{N}}}}=4\). b Dispersion of the energy band ϵk along high-symmetry lines and density of states (DOS). We explore filling levels indicated by the red band between half filling (μ ~ 0) and Van Hove filling (μ = −5.5 meV, ~1/4 filling). c Nearly nested Fermi surface close to the Van Hove energy. Numbered open circles show the momentum resolution of the Fermi surface employed in the FRG approach. d Extended interaction parameters as function of distance rn/a, a is the moiré lattice spacing34,39. The gray area marks the explored range of interactions, which are tunable by the environment. Inset: sketch of the interactions on the triangular moiré lattice. e Phase diagram as function of nearest-neighbor repulsion V1 and filling controlled by the chemical potential μ extracted from the Fermi-liquid instabilities in the FRG flow for U = 4t1, V2/V1 = 0.357 and V3/V1 = 0.260 (see Fig. 1d). The background color encodes our estimate for the critical temperature (see methods). We find an instability towards a valley density wave (VDW) near Van Hove filling μ ≈ − 5.5 meV. It is flanked by an i-wave superconducting instability (iSC) for small V1/U and by a g-wave superconducting instability for larger V1/U. In the ground state, the g-wave superconducting instability forms a chiral g + ig state (g + ig SC) characterized by winding number \({{{\mathcal{N}}}}=4\). The regime colored in black, does not show any instability within our numerical accuracy. f Filling-dependent critical temperature near Van Hove filling. Below the gray band, we cannot resolve any instability within our numerical accuracy. g Dependence of the critical temperature on the extended interactions tuned via V1/U at fixed V2/V1, V3/V1 slightly away from Van Hove filling.

In a range of small twist angles, isolated and narrow moiré bands emerge in TMD hetero-bilayers of WX2/MoX\({}_{2}^{\prime}\) (X,X’ = S,Se)33,34,35,36,37,38,39. These flat bands are formed by the highest, spin-polarized valence band of WX2 and can be described by an extended triangular-lattice Hubbard model H = H0 + HI, which features an effective SU(2) valley symmetry34

$${H}_{0}=\mathop{\sum}\limits_{v=\pm }\left[\mathop{\sum}\limits_{ij}{t}_{i-j}{c}_{i,v}^{{\dagger} }{c}_{j,v}-\mu \mathop{\sum}\limits_{i}{c}_{i,v}^{{\dagger} }{c}_{i,v}\right]$$
(1)
$${H}_{{{{\rm{I}}}}}=U\mathop{\sum}\limits_{i}{n}_{i,+}{n}_{i,-}+\mathop{\sum}\limits_{ij}{V}_{i-j}{n}_{i}{n}_{j}.$$
(2)

Here, ni = ∑vni,v and \({n}_{i,v}={c}_{i,v}^{{\dagger} }{c}_{i,v}\) is the number of electrons on site i with valley index ± , \({c}_{i,v}^{({\dagger} )}\) are the corresponding annihilation (creation) operators. Ref. 34 derived this model for WSe2/MoSe2 and for WSe2/MoS2 with R- and H-stacking configurations. The hopping amplitudes tn for the nth-nearest neighbors depend on the twist angle and we consider typical values for vanishingly small twist angle t1 ≈ 2.5 meV, t2 ≈ −0.5 meV, t3 ≈ −0.25 meV34. The resulting moiré band ϵk features a Van Hove peak in the density of states near 1/4 filling (−5.5 meV), where the Fermi surface is approximately nested (Fig. 1b, c). In experiment, the filling can be adjusted, and Van Hove filling can be reached, by tuning the gate voltage, which we model here by varying the chemical potential μ between 1/4 and 1/2 filling. The interaction parameters U, Vn also depend on the twist angle, and on the dielectric environment so that the strength and range of interactions can be controlled43. First-principles calculations show that the extended interactions Vn are sizable in effective models for hetero-bilayer TMDs34. For our analysis, we use an intermediate interaction strength for the onsite interaction U = 4t1 and explore the effect of further-ranged interactions by varying V1/U [0, 0.5] with V2/V1 ≈ 0.357 and V3/V1 ≈ 0.26034,39 (Fig. 1d, Vn≥4 = 0). In a second step we also investigate the impact of an additional nearest-neighbor exchange interaction HJ = Ji, jSiSj to model strong-coupling effects.

Results

Correlated phase diagram

We study the correlated phases of hetero-bilayer TMDs that emerge out of a metallic state within an itinerant scenario using the fermionic functional renormalization group (FRG)44,45. This method has been successfully applied to scenarios which aim at a realistic modeling of various materials, see, e.g.46,47,48 for some recent contributions or the review49 for a more exhaustive list of references. The FRG resolves the competition between different ordering tendencies in an unbiased way and is employed to calculate the dressed, irreducible two-particle correlation function V(k1, k2, k3, k4) for electrons with momenta ki, i = 1…4, on the Fermi surface (Fig. 1c). It is valley-independent due to the SU(2) valley symmetry (see Methods). Upon lowering the temperature, V(k1, k2, k3, k4) develops sharp, localized peaks for characteristic momentum combinations, indicating long-ranged correlations in real space. This allows us to extract the temperature where a strongly-correlated state forms, as well as the symmetry and type of the strongest correlations (see Methods).

In our model for hetero-bilayer TMD moiré systems, instabilities near 1/4 filling occur due to the high density of states and approximate nesting, which leads to symmetry-broken ground states. We start with varying μ and Vn and calculate the phase diagram based on the two-particle correlation functions (Fig. 1e). Closest to Van Hove filling μ ≈ −5.5 meV, we find that correlations corresponding to a valley density wave (VDW) are strongest, which manifest themselves by peaks at the nesting momenta Qα, α = 1, 2, 3 in V(k1, k2, k3, k4), i.e. at k3 − k1 = Qα or k3 − k2 = Qα. This state is the analogue of a spin density wave50,51 considering that, here, the SU(2) symmetry belongs to a pseudo-spin formed by the valleys. The VDW instability is insensitive towards the inclusion of Vn in the explored range.

Moving the filling slightly away from Van Hove filling, we obtain a superconducting instability, which is indicated by diagonal peak positions, i.e. V(k1, k2, k3, k4) ≈ V(k1, − k1, k3, − k3), that correspond to electron pairs with a total momentum of zero k1 + k2 = k3 + k4 = 0 (Fig. 2b). Increasing the filling further reduces the critical temperature until it vanishes (Fig. 1f). Its origin lies in the strong particle-hole fluctuations due to nesting, which induce an attraction in the pairing channel (see Methods). The inclusion of Vn has a profound impact: it strongly affects the symmetry of the superconducting correlations, because it penalizes electrons to be simultaneously on neighboring sites, so that electron pairing is shifted to farther-distanced neighbors. As a result, the largest attraction is promoted to occur in a higher-harmonic channel.

Fig. 2: Topological superconductivity in hetero-bilayer TMDs.
figure 2

a The phase winding of the superconducting gap in the g + ig state can be visualized by a skyrmion configuration constructed from the vector \({{{\bf{m}}}}=({{{\rm{Re}}}}{{{\Delta }}}_{{{{\bf{k}}}}},{{{\rm{Im}}}}{{{\Delta }}}_{{{{\bf{k}}}}},{\xi }_{{{{\bf{k}}}}})/{({\xi }_{{{{\bf{k}}}}}^{2}+| {{{\Delta }}}_{{{{\bf{k}}}}}{| }^{2})}^{1/2}\). m points up (down) at the highest (lowest) energies ξk = ϵk − μ and rotates \({{{\mathcal{N}}}}\) times in the plane in the vicinity of the Fermi level ξk = 0. The skyrmion configuration is used to calculate the winding number \({{{\mathcal{N}}}}\), which for a broad range of fillings around Van Hove filling is \(| {{{\mathcal{N}}}}| =4\), indicating an enhanced response in thermal and spin quantum Hall measurements. b Functional RG data of the irreducible two-particle correlation function V(k1, k2, k3, k4) near the instability temperature for incoming wave vectors k1, k2. Wave vectors are labeled by the patch points along the Fermi surface indicated in Fig. 1c. The outgoing wave vector k3 is fixed at patch no. 1 and k4 = k1 + k2 − k3 is given by momentum conservation. The sharp diagonal features occur at k1 = − k2, k3 = − k4, indicating the formation of long-ranged superconducting correlations. c Superconducting form factors g±(k) extracted from V(k1, k2, k3, k4) in b. They exhibit a large overlap with a linear combination of the second-nearest-neighbor lattice harmonics g1(k), g2(k) (solid gray lines) defined in the text, which belong to the two-dimensional irreducible representation E2 of the lattice symmetry group C6v. We classify them as g-wave form factors due to their eight nodes. d Absolute value and phase of the gap function on the Fermi surface. The chiral superposition Δk = Δ(g1(k) ± ig2(k)) fully gaps the Fermi surface, thereby minimizing the energy. Such a g + ig superconducting state breaks time-reversal symmetry and is topological with a four-fold phase winding along the Fermi surface \(| {{{\mathcal{N}}}}| =4\). e Stability of the g-wave superconduting state towards inclusion of J/U for V1/U = 0.2 at μ = 5.3075 meV. For growing values of the exchange interaction J, the nearest-neighbor harmonics d1 and d2 of E2 (defined in the text) start to contribute as indicated by the colored transition. They are pure d-wave form factors with only four nodes. For J/U 0.1, the contribution from d1 and d2 is negligible. f Example of the extracted form factors for J/U = 0.5 where d1, d2 and g1, g2 roughly contribute by equal amounts, showing the change in the number of nodes due to the admixture.

The symmetry of the pair correlations can be classified in terms of the irreducible representations of the lattice point group C6v by expanding the eigenfunctions of \(V({{{\bf{k}}}},-{{{\bf{k}}}},{{{\bf{k}}}}^{\prime} ,-{{{\bf{k}}}}^{\prime} )\) in lattice harmonics. Within an irreducible representation, lattice harmonics with the same symmetry but different angular-momentum form factors can mix and it depends on microscopic details which lattice harmonics are the strongest. The valley symmetry follows from the properties of the lattice harmonics under inversion: even (odd) lattice harmonics correspond to valley singlet (triplet) pairing. We obtain two valley-singlet pairing regimes with different spatial symmetries depending on Vn. For small Vn, we find a small regime with A2 symmetry (i-wave) in agreement with previous results for Vn = 052,53. However, for larger Vn (V1/U 0.15), including realistic values in twisted TMDs34,39, we unveil a large regime with E2 symmetry, which contains both, d- and g-wave form factors. The pair correlations in this regime are fitted well using the second-nearest-neighbor lattice harmonics \({g}_{1}({{{\bf{k}}}})=8/9[-\cos (3{k}_{x}/2)\cos (\sqrt{3}{k}_{y}/2)+\cos (\sqrt{3}{k}_{y})]\), \({g}_{2}({{{\bf{k}}}})=8/(3\sqrt{3})\sin (3{k}_{x}/2)\sin (\sqrt{3}{k}_{y}/2)\) (Fig. 2c). We can categorize them as g-wave based on their eight nodes and we note that the number of nodes on the Fermi surface depends on the chemical potential. It is four if the Fermi surface is closer to Γ (where we do not find a superconducting instability), and eight near Van Hove energy (where we do find one). The eight nodes are favored over the nearest-neighbor (d-wave) harmonics with four nodes due to the strong bare nearest-neighbor repulsion. That V1 drives this instability can also be seen at the critical temperature, which initially increases with V1 and then saturates (see Fig. 1g). The selection of the higher lattice harmonics due to V1 has measurable consequences for the topological properties of the superconducting state.

Chiral superconductivity and quantized Hall responses

Due to the two-dimensional E2 symmetry, the superconducting gap has two components Δ1, Δ2 and additional symmetries besides U(1) can be broken depending on the configuration that forms the ground state54. The ground-state configuration is determined by minimizing the Landau energy functional

$${{{\mathcal{L}}}}\,=\,\alpha (| {{{\Delta }}}_{1}{| }^{2}\,\,+\,| {{{\Delta }}}_{2}{| }^{2})\,+\,\beta {(| {{{\Delta }}}_{1}{| }^{2}+| {{{\Delta }}}_{2}{| }^{2})}^{2}\,+\,\gamma | {{{\Delta }}}_{1}^{2}\,+\,{{{\Delta }}}_{2}^{2}{| }^{2}.$$
(3)

We use our FRG results as an input for the effective interaction \(V({{{\bf{k}}}},-{{{\bf{k}}}},{{{\bf{k}}}}^{\prime} ,-{{{\bf{k}}}}^{\prime} ){c}_{{{{\bf{k}}}}^{\prime} ,v}^{{\dagger} }{c}_{-{{{\bf{k}}}}^{\prime} ,v^{\prime} }^{{\dagger} }{c}_{-{{{\bf{k}}}},v^{\prime} }{c}_{k,v}\) close to the instability and perform a Hubbard-Stratonovich transformation to describe the coupling between electrons and pairing fields \({H}_{{{{\rm{P}}}}}={{{\Delta }}}_{1}^{* }({{{\bf{q}}}}){g}_{1}({{{\bf{k}}}}){c}_{{{{\bf{k}}}}+{{{\bf{q}}}},+}{c}_{-{{{\bf{k}}}},-}+{{{\Delta }}}_{2}^{* }({{{\bf{q}}}}){g}_{2}({{{\bf{k}}}}){c}_{{{{\bf{k}}}}+{{{\bf{q}}}},+}{c}_{-{{{\bf{k}}}},-}+\) c.c. Integrating out the electrons, we find in particular γ > 0. Thus, the chiral configuration Δ1 = iΔ2 minimizes the energy. Such a “g + ig” superconducting state breaks time-reversal symmetry and is topologically non-trivial. The Fermi surface is fully gapped as we can see from the quasiparticle energy \({E}_{{{{\bf{k}}}}}={({\xi }_{{{{\bf{k}}}}}^{2}+| {{{\Delta }}}_{{{{\bf{k}}}}}{| }^{2})}^{1/2}\), where ξk = ϵk − μ, Δk = Δ[g1(k) + ig2(k)], and g1, g2 are the FRG-extracted form factors (see Fig. 2c, d). In the superconducting gap Δk, an overall trivial phase of \({{\Delta }}\in {\mathbb{C}}\) can be removed, in contrast to the relative phase difference of π/2 between g1 and g2.

The topological properties can be classified by an integer invariant based on the Skyrmion number55,56,57

$${{{\mathcal{N}}}}=\frac{1}{4\pi }{\int}_{{{{\rm{BZ}}}}}{d}^{2}k\,{{{\bf{m}}}}\cdot \left(\frac{\partial {{{\bf{m}}}}}{\partial {k}_{x}}\times \frac{\partial {{{\bf{m}}}}}{\partial {k}_{y}}\right),$$
(4)

where the pseudo-spin vector is given by \({{{\bf{m}}}}=({{{\rm{Re}}}}{{{\Delta }}}_{{{{\bf{k}}}}},{{{\rm{Im}}}}{{{\Delta }}}_{{{{\bf{k}}}}},{\xi }_{{{{\bf{k}}}}})/{E}_{{{{\bf{k}}}}}\). The integral over the Brillouin zone is non-zero because m follows the phase winding of the superconducting gap around the Fermi surface (Fig. 2a). We calculate \({{{\mathcal{N}}}}\) for the entire range of fillings and find \(| {{{\mathcal{N}}}}| =4\) in the relevant regime where the superconducting instability occurs. Importantly, the high winding number \(| {{{\mathcal{N}}}}| =4\) implies stronger experimental signatures compared to other topological superconductors. Four chiral edge modes appear (as illustrated in Fig. 1a) and the quantized thermal and spin response is enhanced with the spin Hall conductance given by \({\sigma }_{xy}^{s}={{{\mathcal{N}}}}\hslash /(8\pi )\) and the thermal Hall conductance by \(\kappa ={{{\mathcal{N}}}}\pi {k}_{B}{3}^{2}/(6\hslash )\)58,59.

The g + ig pairing state is robust with respect to the inclusion of small to intermediate nearest-neighbor exchange J. For example, for V1/U = 0.2 and μ ≈ −5.31 meV, g + ig pairing is dominant up to reasonably large values of J ≈ 0.1 (Fig. 2e). For larger values of J, d-wave contributions from the nearest-neighbor form factors \({d}_{1}({{{\bf{k}}}})=2(\cos {k}_{x}-\cos ({k}_{x}/2)\cos (\sqrt{3}{k}_{y}/2))\) and \({d}_{2}({{{\bf{k}}}})=3/\sqrt{3}\sin ({k}_{x}/2)\sin (\sqrt{3}{k}_{y}/2)\) start to mix with the previous g-wave ones g1, g2 (see Fig. 2e, f). This is expected when the attraction mediated by antiferromagnetic fluctuations from J overcomes the repulsion from Vn. For larger values of J and farther away from Van Hove filling, the d-wave form factors dominate (see Fig. 3a, b). Then, the superconducting ground state is a fully gapped d + id state with \(| {{{\mathcal{N}}}}| =2\), which can be seen from the phase of the superconducting state winding two times along the Fermi surface (see Fig. 3c, d).

Fig. 3: Effect of exchange coupling.
figure 3

a Filling-dependent critical temperature for a sizable exchange coupling J/t = 0.5. The exchange coupling generates a d-wave instability, which indicates a chiral d + id ground state. b Superconducting form factors d±(k) extracted from V(k1, k2, k3, k4) which exhibit a large overlap with a linear combination of nearest-neighbor lattice harmonics d1, d2 (solid gray lines) defined in the text. Their chiral superposition d ± id fully gaps the Fermi surface with the gap function Δ(k) shown in the lowest panels (c, d). The phase winds twice around the Fermi surface.

Discussion

Most numerical calculations in our work have been carried out for tight-binding and interaction parameters as explicitly suggested for the untwisted WSe2/MoS2 moiré system. However, our calculations with varied parameters suggest that the the correlated phase diagram is very similar in related systems that can be described by an extended Hubbard model on the triangular lattice with sizable longer-ranged density-density interactions near van Hove filling. Our results highlight twisted hetero-bilayers of TMDs as prime candidates for exotic topological superconducting states in two-dimensional materials. They allow — by moiré or substrate engineering — for a very high level of external control17 and our identification of topological g + ig superconductivity opens up pathways to interrogate this elusive phase of matter in a highly tunable setup. Exploiting the unparalleled level of control, these platforms provide the opportunity to scrutinize topological phase transitions using gating which, as we showed, drives the g + ig state into a density wave or a metallic state at either side of the topological superconductor. This can also shed light on related questions about chiral superconductivity in Van-Hove-doped graphene60,61,62,63 or about topological transitions, e.g., the nodal structure at the transition point from \(| {{{\mathcal{N}}}}| =2\) to \(| {{{\mathcal{N}}}}| =4\), which is intensely debated for NaxCoO264,65.

An intriguing avenue of future research concerns the relevance of disorder and finite size effects onto the g-wave superconducting state as moiré materials tend to form localized dislocations66. Similar to the d + id state, we expect the g + ig state to be robust against non-magnetic disorder67,68. However, twist-angle disorder can suppress the superconducting state and might be the reason why a clear observation of superconductivity is so far elusive in moiré TMDs (see however33). Another important pragmatic issue is the relevance of relaxation69,70. Here, we primarily investigated model parameters relevant to different-chalcogen heterobilayers, e.g., untwisted WSe2/MoS2, where the modeling of the moiré potential seems to be robust in comparison with experimental observations34,71. The description of the material as an effective triangular lattice Hubbard model should be valid even for other material combinations, but nevertheless it will be very interesting to consider the effects of lattice reconstruction, e.g., in untwisted WSe2/MoS2, and determine how the model parameters are tuned quantitatively in detail in future work.

From a theoretical angle, the recently developed real-space extension of the unbiased renormalization scheme used here might provide insights into these questions as well as clarify the relevance of physical boundaries72. Another interesting possibility is to investigate if non-local Coulomb interactions can also induce topological triplet superconductivity as discussed for ref. 47. Substrate engineering in the field of two-dimensional systems in general and tunable moiré materials in particular is another fascinating route to follow. Building on our work (highlighting the relevance of long-ranged interactions for exotic collective phases of matter) it will be interesting to generalize the analysis to substrates allowing to engineer long-ranged, logarithmic Keldysh-type of interactions73 and scrutinize their effects on collective phases of matter in moiré materials.

Experimentally, Andreev reflection74 and Raman scattering75 provide way to distinguish the fully gapped chiral from an s-wave supercondutor. The prediction of g + ig topological superconductivity can also be verified using thermal or spin quantum Hall measurements which reveals the four-fold nature of the chiral and topologically protected edge modes. Domain walls between \({{{\mathcal{N}}}}=+4\) and −4 configurations must host eight propagating chiral modes55. Whether these edge modes can be utilized for future quantum information technologies requires additional investigation67. The option appears particularly intriguing with twisted hetero-bilayers of TMDs being so highly tunable and the energy scales on which the material properties can be altered being so low due to the flat bands.

Methods

Functional renormalization group approach

We have employed the functional renormalization group method to explore the phase diagram of our model44,45,49. Within the FRG, we choose the temperature as the flow parameter and use an approximation that neglects feed-back from the self-energy and three-particle vertices or higher. In this approximation, we obtain a renormalization group equation for the two-particle correlation function Γ(2p) that describes its evolution upon lowering the temperature. In an SU(2)-symmetric system, Γ(2p) can be expressed via a (pseudo-)spin-independent coupling function V as \({{{\Gamma }}}_{{s}_{1}{s}_{2}{s}_{3}{s}_{4}}^{(2p)}({{{{\bf{k}}}}}_{1},{{{{\bf{k}}}}}_{2},{{{{\bf{k}}}}}_{3},{{{{\bf{k}}}}}_{4})=V({{{{\bf{k}}}}}_{1},{{{{\bf{k}}}}}_{2},{{{{\bf{k}}}}}_{3},{{{{\bf{k}}}}}_{4}){\delta }_{{s}_{1}{s}_{3}}{\delta }_{{s}_{2}{s}_{4}}-V({{{{\bf{k}}}}}_{1},{{{{\bf{k}}}}}_{2},{{{{\bf{k}}}}}_{4},{{{{\bf{k}}}}}_{3}){\delta }_{{s}_{1}{s}_{4}}{\delta }_{{s}_{2}{s}_{3}}\), where si labels the (pseudo-)spin here given by the valley, and k1, k2 are incoming and k3, k4 outgoing momenta. Momentum conservation requires k1 + k2 = k3 + k4, so for brevity we will use V(k1, k2, k3) = V(k1, k2, k3, k1 + k2 − k3) in the following.

The RG equation for the temperature evolution of V(k1, k2, k3) can then be written as

$$\frac{d}{dT}V={\tau }_{{{{\rm{pp}}}}}+{\tau }_{{{{\rm{ph,d}}}}}+{\tau }_{{{{\rm{ph,cr}}}}}.$$
(5)

with contributions from the particle-particle (pp), the direct paricle-hole (ph,d), and the crossed particle-hole (ph,cr) channel on the right hand side. The corresponding diagrams are visualized in Fig. 4. They are given by

$${\tau }_{{{{\rm{pp}}}}}=-\frac{1}{2}{\int}_{{{{\rm{BZ}}}}}{d}^{2}kV({{{{\bf{k}}}}}_{1},{{{{\bf{k}}}}}_{2},{{{\bf{k}}}})L({{{\bf{k}}}},{{{{\bf{q}}}}}_{{{{\rm{pp}}}}})V({{{\bf{k}}}},{{{{\bf{q}}}}}_{{{{\rm{pp}}}}},{{{{\bf{k}}}}}_{3}),$$
(6)

where we used the short hand \({\int}_{{{{\rm{BZ}}}}}{d}^{2}k=-{A}_{{{{\rm{BZ}}}}}^{-1}\int {d}^{2}k\) and ABZ is the area of the Brillouin zone. The particle-hole contributions read

$$\begin{array}{l}{\tau }_{{{{\rm{ph,d}}}}}=\frac{1}{2}{\int}_{{{{\rm{BZ}}}}}{d}^{2}k\left[2V({{{{\bf{k}}}}}_{1},{{{\bf{k}}}},{{{{\bf{k}}}}}_{3})L({{{\bf{k}}}},{{{{\bf{q}}}}}_{{{{\rm{d}}}}})V({{{{\bf{q}}}}}_{{{{\rm{d}}}}},{{{{\bf{k}}}}}_{2},{{{\bf{k}}}})\right.\\\qquad\quad -\,V({{{\bf{k}}}},{{{{\bf{k}}}}}_{1},{{{{\bf{k}}}}}_{3})L({{{\bf{k}}}},{{{{\bf{q}}}}}_{{{{\rm{d}}}}})V({{{{\bf{q}}}}}_{{{{\rm{d}}}}},{{{{\bf{k}}}}}_{2},{{{\bf{k}}}})\\ \left.\qquad\quad-\,V({{{\bf{k}}}},{{{{\bf{k}}}}}_{1},{{{{\bf{k}}}}}_{3})L({{{\bf{k}}}},{{{{\bf{q}}}}}_{{{{\rm{d}}}}})V({{{{\bf{k}}}}}_{2},{{{{\bf{q}}}}}_{{{{\rm{d}}}}},{{{\bf{k}}}})\right],\end{array}$$
(7)

and

$${\tau }_{{{{\rm{ph,cr}}}}}=-\frac{1}{2}{\int}_{{{{\rm{BZ}}}}}{d}^{2}kV({{{\bf{k}}}},{{{{\bf{k}}}}}_{2},{{{{\bf{k}}}}}_{3})L({{{\bf{k}}}},{{{{\bf{q}}}}}_{{{{\rm{cr}}}}})V({{{{\bf{k}}}}}_{1},{{{{\bf{q}}}}}_{{{{\rm{cr}}}}},{{{\bf{k}}}}).$$
(8)

In these expressions, we introduced qpp = −k + k1 + k2, qd = k + k1 − k3, qcr = k + k2 − k3, and the loop kernel

$$L({{{\bf{k}}}},\pm {{{\bf{k}}}}\,+\,{{{{\bf{k}}}}}^{\prime})\,=\,\frac{d}{dT}\left[T\,\mathop{\sum}\limits_{{{{\rm{i}}}}\omega }{G}_{0}({{{\rm{i}}}}\omega ,{{{\bf{k}}}}){G}_{0}(\pm {{{\rm{i}}}}\omega ,\pm {{{\bf{k}}}}\,+\,{{{{\bf{k}}}}}^{\prime})\right],$$
(9)

with the free propagator \({G}_{0}({{{\rm{i}}}}\omega ,{{{\bf{k}}}})={[{{{\rm{i}}}}\omega -{\xi }_{{{{\bf{k}}}}}]}^{-1}\). In these expressions, we have neglected the (external) frequency dependence assuming that the strongest correlations occur for the lowest Matsubara frequencies.

Fig. 4: Loop diagrams.
figure 4

Contributing to the functional RG flow, cf. Eq. (5).

For the numerical implementation, we resolve the momentum dependence in a so-called patching scheme that divides the Fermi surface into N pieces based on equidistant angles and treats the radial dependence for a fixed angle as constant. This accurately describes the relevant momentum dependence, which is along the Fermi surface. In our numerical calculations, we have chosen between N = 48 and N = 96 patches, cf. Fig. 1b. Our results on the type of instability do not depend on this choice and the quantitative results for critical temperatures vary only mildly with N.

The initial condition for V(k1, k2, k3) at high temperatures is given by the Fourier transform of HI in Eq. (2). We set \({T}_{0}=\max ({\epsilon }_{{{{\bf{k}}}}})\) as starting temperature. We then calculate the temperature evolution of V(k1, k2, k3) according to Eq. (5) by solving the integro-differential equation. As described above, the development of strong correlations is signaled by a diverging V(k1, k2, k3) at a critical temperature Tc. Our numerical criterion to detect the divergence is a convex temperature dependence and max[V(k1, k2, k3)] exceeding 30t1. Tc would be the mean-field critical temperature in an RPA resummation, however, here the estimate is slightly improved due to the inclusion of the coupling between different channels. The Fermi liquid is stable within our numerical accuracy if no divergence occurs before Tl = 2  10−4t1 is reached. In the cases when correlated states develop, we can read off the type of correlations from the momentum structure of V(k1, k2, k3) at Tc. Up to an overall constant, this determines the effective interaction close to the instability and directly suggests the order-parameter corresponding to the instability. Following this procedure for an extended range of parameters, we obtain the presented phase diagrams. We have also checked that the qualitative features of the phase diagrams remain the same when we vary the overall interaction strength U. We show examples for U = 3t1 and U = 5t1 in Fig. 5.

Fig. 5: U dependence.
figure 5

Critical temperature for three values of the on-site interaction U and fixed Vn. The valley density wave (VDW) and g-wave pairing phase (g + ig) appear near Van Hove filling in all three cases.

Fig. 6: Flow of the effective interaction.
figure 6

a Representative scale-dependence of the maximal value of the effective interaction \({V}_{\max }=\max V({{{{\bf{k}}}}}_{1},{{{{\bf{k}}}}}_{2},{{{{\bf{k}}}}}_{3},{{{{\bf{k}}}}}_{4})\) in the parameter region of the g + ig SC instability. During the different stages of the RG flow, the different interaction channels (i.e. the pp-channel and the ph-channels) contribute. The initial increase of \({V}_{\max }\) for \(\log T/t\gtrsim -4\) is driven by the particle-hole fluctuations. Sufficiently far away from nesting, these fluctuations are not critical, i.e. \({V}_{\max }\) nearly saturates (cf. plateau for \(-4\gtrsim \log T/t\gtrsim -6\)). Eventually, particle-particle fluctuations lead to a divergence of the flow. b Interaction V with conventions as described in Fig. 2b for T = T0. ce Interaction V for the three choices of T as indicated in panel (a), exhibiting different dominant features. The strongest vertical and horizontal features in c,d correspond to valley fluctuations with momentum transfer of k1 − k3 = Qα and k2 − k3 = Qα (α = 1, 2, 3), diagonal features in e are typical for superconductivity with k1 + k2 = k3 + k4 = 0.

An advantage of the FRG approach over other methods such as mean-field or the random phase approximation is that it takes the coupling between different ordering channels into account on equal footing. In particular, it can resolve the generation of an attraction in the pairing channel due to particle-hole fluctuations. In this way the bosons mediating the pairing can be read off in an unbiased way, i.e. without choosing specific channels as relevant by hand. In our case, fluctuations due to nesting increase the interaction in the particle-hole channel via τph,d/cr upon lowering the temperature (ph diagrams in Fig. 4), which then feed back into the pairing channel via τpp (pp diagram in Fig. 4) and yield the effective pairing interaction shown in Fig. 2b. This can be visualized by snapshots of the interaction during the flow (Fig. 6), which first develops horizontal features typical for particle-hole fluctuations with momentum transfer Qα and which are associated with the VDW instability, Eventually, diagonal features of the effective interaction are induced that correspond to the superconducting pairing instability.

Landau functional

To extract the form factors of the superconducting instabilities, we diagonalize \(V({{{\bf{k}}}},-{{{\bf{k}}}},{{{\bf{k}}}}^{\prime} )\), keep the eigenfunction(s) with the largest eigenvalue and approximate it by lattice harmonics so that \(V({{{\bf{k}}}},-{{{\bf{k}}}},{{{\bf{k}}}}^{\prime} )\approx -\lambda {\sum }_{i = 1,2}\int d{{{\bf{k}}}}d{{{\bf{k}}}}^{\prime} {g}_{i}({{{\bf{k}}}}){g}_{i}({{{\bf{k}}}}^{\prime} ){c}_{{{{\bf{k}}}}^{\prime} +}^{{\dagger} }{c}_{-{{{\bf{k}}}}^{\prime} -}^{{\dagger} }{c}_{-{{{\bf{k}}}}-}{c}_{{{{\bf{k}}}}+}\) with λ > 0. We have used the extracted lattice harmonics to derive the Landau functional (3) from our microscopic model. The decisive prefactor of the term \(| {{{\Delta }}}_{1}^{2}+{{{\Delta }}}_{2}^{2}{| }^{2}\) is given by

$$\gamma =T\mathop{\sum}\limits_{{{{\rm{i}}}}\omega }{\int}_{{{{\rm{BZ}}}}}{d}^{2}k\frac{{g}_{1}{({{{\bf{k}}}})}^{2}{g}_{2}{({{{\bf{k}}}})}^{2}}{{({{{\rm{i}}}}\omega \,-\,{\xi }_{{{{\bf{k}}}}})}^{2}{({{{\rm{i}}}}\omega \,+\,{\xi }_{{{{\bf{k}}}}})}^{2}} \quad={\int}_{{{{\rm{BZ}}}}}{d}^{2}k\,{g}_{1}{({{{\bf{k}}}})}^{2}{g}_{2}{({{{\bf{k}}}})}^{2}\frac{1-2{n}_{{{{\rm{F}}}}}({\xi }_{{{{\bf{k}}}}})\,+\,2{\xi }_{{{{\bf{k}}}}}{n^{\prime}_{\rm{F}}} ({\xi }_{{{{\bf{k}}}}})}{4{\xi }_{{{{\bf{k}}}}}^{3}}$$
(10)

with the Fermi function nF. We have calculated γ numerically and found it to be positive in the considered range of μ and T. As an analytical estimate for γ, we can approximate the dispersion by ξ ≈ k2/(2m) − μ with density of states ρϵ, and the form factors by \({g}_{1}=\cos (n\varphi )\), \({g}_{2}=\sin (n\varphi )\) with \(\varphi =\arctan {k}_{y}/{k}_{x}\) and n = 4 for g + ig superconductivity (n = 2 for d + id and n = 1 for p + ip). With this simplification, we obtain

$$\gamma \approx \int \,d\varphi {\cos }^{2}(n\varphi ){\sin }^{2}(n\varphi )\,\int \,d\epsilon {\rho }_{\epsilon }\frac{\frac{1}{2}\,-\,{n}_{{{{\rm{F}}}}}(\epsilon )\,+\,\epsilon {n^{\prime}_{\rm{F}}} (\epsilon )}{2{\epsilon }^{3}} =\frac{m}{16\pi {T}^{2}}\int \frac{\sinh (x)-x}{4{x}^{3}(1\,+\,\cosh (x))}\approx 0.05\frac{m}{16\pi {T}^{2}},$$
(11)

which we can take as a rough estimate for γ if μ is away from the Van Hove energy. Right at the Van Hove energy, an additional logarithmic dependence on Tc emerges.