We consider the low-energy Hamiltonian for the metallic bands of monolayer NbSe2<\/p>\n
$$\\hat{H}={\\sum} _{{{\\bf{k}}},\\sigma }{\\xi }_{{{\\bf{k}}},\\sigma }{\\hat{c}}_{{{\\bf{k}}},\\sigma }^{{\\dagger} }{\\hat{c}}_{{{\\bf{k}}},\\sigma }+\\hat{V},$$<\/p>\n
\n (7)\n <\/p>\n
where the first term describes independent Bloch electrons in the bands crossing the Fermi level with energies \u03bek,\u03c3, and \\({\\hat{c}}_{{{\\bf{k}}},\\sigma }^{({\\dagger} )}\\) are the creation and annihilation operators of Bloch electrons with crystal momentum k. The second term is the statically screened Coulomb interaction which, in the position representation (1\u00a0\u2261\u00a0r,\u00a0\u03c3), is defined by a Dyson equation for its matrix elements \\(V(1,2)=\\left\\langle (1,2)\\right\\vert \\hat{V}\\left\\vert (1,2)\\right\\rangle\\) as<\/p>\n
$$V(1,2)={V}_{0}(1,2)+\\int\\,d(3,4){V}_{0}(1,3)P(3,4)V(4,2),$$<\/p>\n
\n (8)\n <\/p>\n
with V0(1,\u00a02) the matrix elements of the unscreened interaction and P the polarisability of the interacting system. We consider a trigonal 2D lattice with periodic boundary conditions composed of N prismatic unit cells of area \u03a9 and height L with sample volume \\({{\\mathcal{V}}}=N\\Omega L\\). The derivation of Eq.\u00a0(2)<\/a> from Eq.\u00a0(8)<\/a> can be found in the Supplementary Material.<\/p>\n We treat the screening at the RPA level, where the full polarisability in Eq. (8<\/a>) is replaced with the non-interacting polarisability tensor given by Eq. (4<\/a>) where<\/p>\n $${\\chi }_{\\sigma }({{\\bf{p}}}+{{\\bf{q}}},{{\\bf{p}}};\\omega )=\\frac{1}{N\\Omega }\\frac{f({\\xi }_{{{\\bf{p}}}+{{\\bf{q}}},\\sigma })-f({\\xi }_{{{\\bf{p}}},\\sigma })}{\\omega+{\\xi }_{{{\\bf{p}}}+{{\\bf{q}}},\\sigma }-{\\xi }_{{{\\bf{p}}},\\sigma }+i\\eta },$$<\/p>\n \n (9)\n <\/p>\n f is the Fermi function evaluated from the chemical potential \u03bc and the overlaps \\({{\\mathcal{F}}}\\) are defined in Eq. (5<\/a>). For momenta \\({{\\bf{k}}}\/{{{\\bf{k}}}}^{{\\prime} }\\) outside the first Brillouin zone, Eq. (5<\/a>) is modified by decomposing \\({{\\bf{k}}}={{\\mathcal{P}}}({{\\bf{k}}})+{{\\mathcal{Q}}}({{\\bf{k}}})\\) (see Supplementary Fig.\u00a03<\/a>), where \\({{\\mathcal{P}}}\\) projects onto the reciprocal lattice, \\({{\\mathcal{Q}}}\\) projects onto the associated crystal momentum and using \\({{{\\bf{u}}}}_{{{\\bf{k}}},\\sigma }({{\\bf{r}}})={e}^{i{{\\mathcal{P}}}({{\\bf{k}}})\\cdot r}{{{\\bf{u}}}}_{{{\\mathcal{Q}}}({{\\bf{k}}}),\\sigma }({{\\bf{r}}})\\). We express the interaction \\(\\hat{V}\\) entering Eq. (7<\/a>) in the Bloch basis. It holds for the statically screened Coulomb interaction between Bloch states (see the\u00a0Supplementary Material<\/a> for the full derivation),<\/p>\n $$\t\\left\\langle {{\\mathcal{Q}}}({{\\bf{k}}}+{{\\bf{q}}}),\\sigma ;{{\\mathcal{Q}}}({{{\\bf{k}}}}^{{\\prime} }-{{\\bf{q}}}),{\\sigma }^{{\\prime} }\\right\\vert \\hat{V}\\left\\vert {{\\bf{k}}},\\sigma ;{{{\\bf{k}}}}^{{\\prime} },{\\sigma }^{{\\prime} }\\right\\rangle=\\\\ \t \\frac{1}{N\\Omega }{\\sum} _{{{\\bf{G}}},{{{\\bf{G}}}}^{{\\prime} }}{V}_{{{\\bf{G}}},{{{\\bf{G}}}}^{{\\prime} }}^{2{{\\rm{D}}},{{\\rm{RPA}}}}({{\\bf{q}}};{0}^{+}){{{\\mathcal{F}}}}_{{{\\bf{k}}}+{{\\bf{q}}},{{\\bf{k}}}}^{\\sigma }(-{{\\bf{G}}}){{{\\mathcal{F}}}}_{{{{\\bf{k}}}}^{{\\prime} }-{{\\bf{q}}},{{{\\bf{k}}}}^{{\\prime} }}^{{\\sigma }^{{\\prime} }}({{{\\bf{G}}}}^{{\\prime} }),$$<\/p>\n \n (10)\n <\/p>\n where \\({{\\bf{k}}},{{{\\bf{k}}}}^{{\\prime} }\\) and q are in-plane momenta restricted to the first Brillouin zone. Noting that the dominant contribution to the matrix element is due to the first few nearest-neighbour sites (see Supplementary Section\u00a0II D<\/a> for a discussion) we introduce an effective interaction potential<\/p>\n $$\t\\left\\langle {{\\mathcal{Q}}}({{\\bf{k}}}+{{\\bf{q}}}),\\sigma ;{{\\mathcal{Q}}}({{{\\bf{k}}}}^{{\\prime} }-{{\\bf{q}}}),{\\sigma }^{{\\prime} }\\right\\vert \\hat{V}\\left\\vert {{\\bf{k}}},\\sigma ;{{{\\bf{k}}}}^{{\\prime} },{\\sigma }^{{\\prime} }\\right\\rangle=\\\\ \t \\approx \\frac{1}{N\\Omega }\\left({\\sum} _{{{\\bf{G}}},{{{\\bf{G}}}}^{{\\prime} }}{V}_{{{\\bf{G}}},{{{\\bf{G}}}}^{{\\prime} }}^{2{{\\rm{D}}},{{\\rm{RPA}}}}({{\\bf{q}}})\\right){{{\\mathcal{F}}}}_{{{\\bf{k}}}+{{\\bf{q}}},{{\\bf{k}}}}^{\\sigma }({{\\bf{0}}}){{{\\mathcal{F}}}}_{{{{\\bf{k}}}}^{{\\prime} }-{{\\bf{q}}},{{{\\bf{k}}}}^{{\\prime} }}^{{\\sigma }^{{\\prime} }}({{{\\bf{0}}}}^{{\\prime} })\\\\ \t=:\\frac{1}{N\\Omega }{V}^{{{\\rm{RPA}}}}({{\\bf{q}}}){{{\\mathcal{F}}}}_{{{\\bf{k}}}+{{\\bf{q}}},{{\\bf{k}}}}^{\\sigma }({{\\bf{0}}}){{{\\mathcal{F}}}}_{{{{\\bf{k}}}}^{{\\prime} }-{{\\bf{q}}},{{{\\bf{k}}}}^{{\\prime} }}^{{\\sigma }^{{\\prime} }}({{{\\bf{0}}}}^{{\\prime} }),$$<\/p>\n \n (11)\n <\/p>\n where we restrict \\({{\\bf{G}}},{{{\\bf{G}}}}^{{\\prime} }\\) to the range where \\({{\\mathcal{F}}}({{\\bf{G}}})\\approx {{\\mathcal{F}}}({{\\bf{0}}})\\) holds, i.e., up to third nearest neighbouring Brillouin zones. This effective interaction potential enables a discussion of the screened interaction on the lattice in terms of a simple q-dependent quantity shown in Fig.\u00a02<\/a>, but is not used for the calculation of the pairing, where we keep the dependence of the Bloch factors \\({{\\mathcal{F}}}\\) on G.<\/p>\n To access the possible superconducting instabilities described by Eq. (6<\/a>), we solve a linearised form of the gap equation. First, we project the gap equations on the Fermi surfaces by introducing a local density of states \u03c1k,\u03c3 along the Fermi surfaces. We restrict the allowed momenta \\({{{\\bf{k}}}}^{{\\prime} }\\) in Eq. (6<\/a>) to a range around the Fermi surfaces given by a cutoff \u00a0\u00b1\u00a0\u039b in energy. Next, we convert the sum over the transverse component \\({{{\\bf{k}}}}^{{\\prime} }-{{{\\bf{k}}}}_{{{\\rm{F}}}}\\) into an integral in energy over this range [\u00a0\u2212\u00a0\u039b,\u00a0\u039b]. Assuming that the gap and the interaction vary slowly in this thin shell around each Fermi surface, we can approximate \\({\\Delta }_{{{{\\bf{k}}}}^{{\\prime} },\\sigma }\\approx {\\Delta }_{{{{\\bf{k}}}}_{{{\\rm{F}}}},\\sigma }\\) and \\({V}_{{{\\bf{k}}},{{{\\bf{k}}}}^{{\\prime} },\\sigma }\\approx {V}_{{{\\bf{k}}},{{{\\bf{k}}}}_{{{\\rm{F}}}},\\sigma }\\). The resulting gap equation with \\({{\\bf{k}}},{{{\\bf{k}}}}^{{\\prime} }\\) on the Fermi surfaces takes the form<\/p>\n $${\\Delta }_{{{\\bf{k}}},\\sigma }=-{\\sum} _{{{{\\bf{k}}}}^{{\\prime} }\\in {{\\rm{FS}}}}{V}_{{{\\bf{k}}},{{{\\bf{k}}}}^{{\\prime} },\\sigma }\\int_{-\\Lambda }^{\\Lambda }d\\xi {\\rho }_{{{{\\bf{k}}}}^{{\\prime} },\\sigma }(\\xi )\\Pi ({E}_{{{{\\bf{k}}}}^{{\\prime} }\\sigma }){\\Delta }_{{{{\\bf{k}}}}^{{\\prime} },\\sigma },$$<\/p>\n \n (12)\n <\/p>\n with the pairing function \\(\\Pi (E)=\\tanh (\\beta E\/2)\/2E\\). Close to Tc we have \u03a0(E)\u00a0\u2248\u00a0\u03a0(\u03be) and Eq. (12<\/a>) becomes a set of linear equations. Instead of starting from the statically screened interaction and applying the mean-field approximation, an alternative approach to the problem is provided by investigating the pair-scattering vertex where the problem also simplifies to the real and static part of the effective pairing interaction18<\/a>,55<\/a>. The problem of finding the shape of the order parameter along the Fermi surfaces is equivalent to determining the kernel of a matrix \\({{\\mathbb{M}}}_{\\sigma }\\), whose indices are the momenta kF on the Fermi surface, as<\/p>\n $${{\\mathbb{M}}}_{\\sigma }{{{\\mathbf{\\Delta }}}}_{\\sigma }=0,$$<\/p>\n \n (13)\n <\/p>\n where, \\(\\forall {{\\bf{k}}},{{{\\bf{k}}}}^{{\\prime} }\\in {{\\rm{FS}}}\\),<\/p>\n $${({{\\mathbb{M}}}_{\\sigma })}_{{{\\bf{k}}},{{{\\bf{k}}}}^{{\\prime} }}={\\delta }_{{{\\bf{k}}},{{{\\bf{k}}}}^{{\\prime} }}+{V}_{{{\\bf{k}}},{{{\\bf{k}}}}^{{\\prime} },\\sigma }{\\rho }_{{{{\\bf{k}}}}^{{\\prime} },\\sigma }{\\alpha }^{0}(T,\\Lambda ),$$<\/p>\n \n (14)\n <\/p>\n $${\\alpha }^{0}(T,\\Lambda )=\\int_{0}^{\\Lambda }d\\xi \\frac{\\tanh (\\beta \\xi \/2)}{\\xi }.$$<\/p>\n \n (15)\n <\/p>\n Here, we assumed a flat density of states \u03c1k,\u2191(\u03be)\u00a0\u2248\u00a0\u03c1k,\u2191 within the energy range \u00a0\u00b1\u00a0\u039b around the chemical potential. We rewrite \\({{\\mathbb{M}}}_{\\sigma }\\) by factoring out the energy integral \u03b10 in the second summand as<\/p>\n $${{\\mathbb{M}}}_{\\sigma }={\\mathbb{1}}+{\\alpha }^{0}(T,\\Lambda ){{\\mathbb{U}}}_{\\sigma },$$<\/p>\n \n (16)\n <\/p>\n where \\({{\\mathbb{U}}}_{\\sigma }\\) is the pairing matrix containing the interaction and local density of states along the Fermi surfaces. The possible superconducting pairings solving Eq. (13<\/a>) are the eigenvectors of \\({{\\mathbb{U}}}_{\\uparrow }\\) with the corresponding eigenvalues \u03bbi fulfilling<\/p>\n $${\\lambda }_{i}=-1\/{\\alpha }^{0}({T}_{{{\\rm{c}}},i},\\Lambda ),$$<\/p>\n \n (17)\n <\/p>\n and Tc,i the critical temperature of the respective instability. The cutoff \u039b sets an energy scale against which to measure the gaps and Tc, similar to the Debye frequency in the usual BCS-theory1<\/a>. In general, not all eigenvectors of \\({{\\mathbb{U}}}_{\\sigma }\\) found in this way correspond to physically stable situations67<\/a>. Similarly, the requirement of vanishing \u0394k,\u03c3 close to Tc,i is only fulfilled if no other superconducting phase with higher Tc is present. As such, we expect only the solution with the highest Tc to be realised. Due to the monotonic nature of \u03b10(T,\u00a0\u039b) as a function of both T and \u039b, the solution with the highest critical temperature always corresponds to the lowest negative eigenvalue \u03bb of the pairing matrix. As \u039b is a model parameter, we classify and discuss the strength of the different instabilities by their eigenvalues \u03bbi.<\/p>\n To investigate how the gaps evolve at temperatures below Tc, we solve Eq. (12<\/a>) without linearizing. Since we want to use the approximations \\({\\Delta }_{{{{\\bf{k}}}}^{{\\prime} },\\sigma }\\approx {\\Delta }_{{{{\\bf{k}}}}_{{{\\rm{F}}}},\\sigma }\\) and \\({V}_{{{\\bf{k}}},{{{\\bf{k}}}}^{{\\prime} },\\sigma }\\approx {V}_{{{\\bf{k}}},{{{\\bf{k}}}}_{{{\\rm{F}}}},\\sigma }\\), we fix the cutoff to the value of \u039b\u00a0=\u00a00.1 eV which, together with the rescaling of the screened interaction to \\({\\hat{V}}_{\\gamma }:=\\gamma \\hat{V}\\) as discussed in the main text, yields Tc\u00a0=\u00a02 K. This enables us to fit to the experimental Tc while retaining the shape of the interaction as found from our calculation of the screened Coulomb interaction. Since the macroscopic polarisability P0,0 calculated from DFT (see Supplementary Section\u00a0I B<\/a>) is even lower at intermediate momentum transfers than our results from a tight-binding calculation, it is plausible that the screened interaction in absence of the substrate is even stronger than our \\(\\hat{V}\\). This is to be expected, since the inclusion of the Se p-orbitals, which contribute most noticeably near the M points, would further diminish the overlaps between states on the Fermi surface, resulting in lower screening.<\/p>\n We expand the gap for one spin species in terms of a restricted set of basis functions, corresponding to the irreducible representations of the symmetry group D3h, as<\/p>\n $${\\forall }_{{{\\bf{k}}}\\in \\pi }:\\quad {\\Delta }_{{{\\bf{k}}},\\uparrow }={\\sum} _{\\mu }{f}_{\\mu }^{\\pi }({{\\bf{k}}}){\\Delta }_{\\mu }^{\\pi },$$<\/p>\n \n (18)\n <\/p>\n where \u03bc labels basis functions and \u03c0 \u2208 {\u0393,\u00a0K} the Fermi surfaces22<\/a>,43<\/a>. We orthonormalise our basis functions with respect to the number N\u03c0 of momenta on each Fermi surface \u03c0 to satisfy<\/p>\n $${\\left\\langle {f}_{\\mu }^{\\pi },{f}_{{\\mu }^{{\\prime} }}^{\\pi }\\right\\rangle }_{\\pi }:=\\frac{1}{{N}_{\\pi }}{\\sum} _{{{\\bf{k}}}\\in \\pi }{f}_{{\\mu }^{{\\prime} }}^{\\pi }({{\\bf{k}}})\\;{f}_{\\mu }^{\\pi }({{\\bf{k}}})={\\delta }_{\\mu,{\\mu }^{{\\prime} }}.$$<\/p>\n \n (19)\n <\/p>\n Using this orthogonality, we reduce the dimensionality of the self-consistency problem drastically by projecting the gap equation into a form that couples only the expansion coefficients \\({\\Delta }_{\\mu }^{\\pi }\\)22<\/a>. It reads<\/p>\n $${\\Delta }_{\\mu }^{\\pi }=-\\frac{1}{{N}_{\\pi }}{\\sum} _{{\\pi }^{{\\prime} },{\\mu }^{{\\prime} }}{\\sum} _{\\begin{array}{c}{{\\bf{k}}}\\in \\pi \\\\ {{{\\bf{k}}}}^{{\\prime} }\\in {\\pi }^{{\\prime} }\\end{array}}{f}_{\\mu }^{\\pi }({{\\bf{k}}}){V}_{{{\\bf{k}}},{{{\\bf{k}}}}^{{\\prime} },\\uparrow }\\,{\\rho }_{{{{\\bf{k}}}}^{{\\prime} },\\uparrow }\\,{\\alpha }_{{{{\\bf{k}}}}^{{\\prime} }}^{{\\pi }^{{\\prime} }}\\,{f}_{{\\mu }^{{\\prime} }}^{{\\pi }^{{\\prime} }}({{{\\bf{k}}}}^{{\\prime} }){\\Delta }_{{\\mu }^{{\\prime} }}^{{\\pi }^{{\\prime} }},$$<\/p>\n \n (20)\n <\/p>\n where the pairing strength \\({\\alpha }_{{{\\bf{k}}}}^{\\pi }\\) is given by<\/p>\n $${\\alpha }_{{{\\bf{k}}}}^{\\pi }=\\int_{0}^{\\Lambda }d\\xi \\,\\frac{\\tanh (\\beta E(\\xi,{\\Delta }_{{{\\bf{k}}}}^{\\pi })\/2)}{E(\\xi,{\\Delta }_{{{\\bf{k}}}}^{\\pi })}.$$<\/p>\n \n (21)\n <\/p>\n At Tc, the pairing strength \\({\\alpha }_{{{\\bf{k}}}}^{\\pi }\\) becomes independent of k and \u03c0. In this case, the orthogonality Eq. (19<\/a>) yields a linearised gap equation for the expansion coefficients<\/p>\n $${\\Delta }_{\\mu }^{\\pi }=-{\\alpha }^{0}{\\sum} _{{\\pi }^{{\\prime} },{\\mu }^{{\\prime} }}{U}_{\\mu,{\\mu }^{{\\prime} }}^{\\pi,{\\pi }^{{\\prime} }}{\\Delta }_{\\mu }^{{\\pi }^{{\\prime} }},$$<\/p>\n \n (22)\n <\/p>\n where the pairing matrix \\({\\mathbb{U}}\\) in this representation is given by<\/p>\n $${U}_{\\mu,{\\mu }^{{\\prime} }}^{\\pi,{\\pi }^{{\\prime} }}=\\frac{1}{{N}_{\\pi }}{\\sum} _{{{\\bf{k}}}\\in \\pi,{{{\\bf{k}}}}^{{\\prime} }\\in {\\pi }^{{\\prime} }}{f}_{\\mu }^{\\pi }({{\\bf{k}}}){V}_{{{\\bf{k}}},{{{\\bf{k}}}}^{{\\prime} },\\uparrow }\\,{\\rho }_{{{{\\bf{k}}}}^{{\\prime} },\\uparrow }{f}_{{\\mu }^{{\\prime} }}^{{\\pi }^{{\\prime} }}({{{\\bf{k}}}}^{{\\prime} }).$$<\/p>\n \n (23)\n <\/p>\n To account for the mixing between different symmetries due to the k dependence of \\({\\alpha }_{{{\\bf{k}}}}^{\\pi }\\) below Tc, we insert a unity in Eq. (20<\/a>) to arrive at<\/p>\n $${\\Delta }_{\\mu }^{\\pi }=-{\\sum} _{{\\pi }^{{\\prime} },{\\mu }^{{\\prime} },{\\mu }^{{\\prime\\prime} }}{U}_{\\mu,{\\mu }^{{\\prime} }}^{\\pi,{\\pi }^{{\\prime} }}{\\Upsilon }_{{\\mu }^{{\\prime} },{\\mu }^{{\\prime\\prime} }}^{{\\pi }^{{\\prime} }}{\\Delta }_{{\\mu }^{{\\prime\\prime} }}^{{\\pi }^{{\\prime} }},$$<\/p>\n \n (24)\n <\/p>\n where the mixing of different basis functions is described by<\/p>\n $${\\Upsilon }_{\\mu,{\\mu }^{{\\prime} }}^{\\pi }=\\frac{1}{{N}_{\\pi }}{\\sum} _{{{\\bf{k}}}\\in \\pi }{f}_{\\mu }^{\\pi }({{\\bf{k}}}){\\alpha }_{{{\\bf{k}}}}^{\\pi }{f}_{{\\mu }^{{\\prime} }}^{\\pi }({{\\bf{k}}}).$$<\/p>\n \n (25)\n <\/p>\n For an incomplete set of basis functions this insertion of unity becomes approximate, but in our case still suffices to capture the weak mixing between basis functions induced by \\({\\alpha }_{{{\\bf{k}}}}^{\\pi }\\). As such, we employ Eq. (24<\/a>) to numerically determine the evolution of the solutions obtained from the linearised equations below their Tc,i. The basis functions we used and the temperature dependence curves for the respective expansion coefficients are shown in Supplementary Section\u00a0III A<\/a>.<\/p>\n We define a purely chiral solution as<\/p>\n $${\\Delta }_{{{\\rm{p}}}+{{\\rm{ip}}}}=\\frac{| | {\\Delta }_{{{\\rm{py}}}}| | }{| | {\\Delta }_{{{\\rm{px}}}}| | }{\\Delta }_{{{\\rm{px}}}}+i{\\Delta }_{{{\\rm{py}}}},$$<\/p>\n \n (26)\n <\/p>\n where \\(| | \\Delta | | :=\\sqrt{{\\sum }_{\\pi=G,K}{\\langle \\Delta,\\Delta \\rangle }_{\\pi }}\\), with the scalar product defined in Eq. (19<\/a>). The \u0394px component is renormalised so that both nematic solutions enter with the same weight. The ratio of admixture Rch of \u0394p+ip in the chiral solution \u0394ch in Fig.\u00a04<\/a>a is defined as<\/p>\n $${R}_{{{\\rm{ch}}}}:=\\frac{| {\\sum}_{\\pi=\\Gamma,K}{\\langle {\\Delta }_{{{\\rm{p}}}+{{\\rm{ip}}}},{\\Delta }_{{{\\rm{ch}}}}\\rangle }_{\\pi }| }{| | {\\Delta }_{{{\\rm{ch}}}}| | \\,| | {\\Delta }_{{{\\rm{p}}}+{{\\rm{ip}}}}| | }.$$<\/p>\n \n (27)\n <\/p>\n The square root of mean squared \u0394 shown in Fig.\u00a04<\/a>a for each gap symmetry is \\(\\sqrt{\\langle {\\Delta }^{2}\\rangle }\\equiv | | \\Delta | |\\).<\/p>\n The free energy of a superconductor with gap \u0394k\u03c3 can be written as<\/p>\n $$F=\\frac{1}{2}{\\sum} _{{{\\bf{k}}}\\sigma }\\left({\\xi }_{{{\\bf{k}}}\\sigma }-{E}_{{{\\bf{k}}}\\sigma }+\\Pi ({E}_{{{\\bf{k}}}\\sigma })| {\\Delta }_{{{\\bf{k}}}\\sigma }{| }^{2}\\right)\\\\ -{k}_{TB}T{\\sum} _{{{\\bf{k}}}\\sigma }{{\\rm{In}}} \\left(1+{e}^{-{E}_{{{\\bf{k}}}\\sigma }\/{k}_{TB}T}\\right).$$<\/p>\n \n (28)\n <\/p>\n Note that this free energy is written assuming that the gap obeys the self-consistent equations (see Supplementary Section\u00a0III<\/a>). It cannot be used to find the physical gaps by minimization with respect to \u0394k\u03c3, but does yield the correct result once a solution is found by solving the self-consistent gap equations67<\/a>. We therefore use it to rank the solutions by their free energy gain within the energy window of \u00a0\u00b1\u00a02\u2009meV where the contributions of different gap symmetries are distinguishable.<\/p>\n NbSe2 monolayers2 were grown on graphitised (bilayer graphene) on 6H-SiC(0001) by molecular beam epitaxy (MBE) at a base pressure of ~ 2 \u22c5 10\u221210\u2009mbar in our home-made ultrahigh-vacuum (UHV) MBE system. SiC wafers with resistivities \u03c1\u00a0~120\u2009\u03a9cm were used68<\/a>.<\/p>\n Reflective high-energy electron diffraction (RHEED) was used to monitor the growth of the NbSe layer2. During the growth, the bilayer graphene\/SiC substrate was kept at 580\u2009\u00b0C. High-purity Nb (99.99%) and Se (99.999%) were evaporated using an electron beam evaporator and a standard Knudsen cell, respectively. The Nb:Se flux ratio was kept at 1:30, while evaporating Se led to a pressure of ~ 5 \u22c5 10\u22129\u2009mbar (Se atmosphere). Samples were prepared using an evaporation time of \\(30\\,\\min\\) to obtain a coverage of \u00a0~\u00a01 ML. To minimise the presence of atomic defects, the evaporation of Se was subsequently kept for an additional 5 minutes. Atomic force microscopy (AFM) under ambient conditions was routinely used to optimise the morphology of the NbSe2 layers. The samples used for AFM characterization were not further used for the scanning tunnelling microscope (STM). Lastly, to transfer the samples from our UHV-MBE chamber to our STM, they were capped with a \u00a0~10\u2009nm amorphous film of Se, which was subsequently removed by annealing at the STM under UHV conditions at 280\u2009\u00b0C.<\/p>\n Scanning tunnelling microscopy\/spectroscopy experiments were carried out in a commercial low-temperature, high magnetic field (11 T) UHV-STM from Unisoku (model USM1300) operated at T\u00a0=\u00a00.34\u2009K unless otherwise stated. STS measurements (using Pt\/Ir tips) were performed using the lock-in technique with typical ac modulations of 20\u201340\u2009\u03bcV at 833\u2009Hz. All data were analyzed and rendered using WSxM software69<\/a>. To avoid tip artifacts in our STS measurements, the STM tips used for our experiments were previously calibrated on a Cu(111) surface. A tip was considered calibrated only when tunnelling spectroscopy performed on Cu(111) showed a sharp surface state onset at \u00a0\u22120.44\u2009eV followed by a clean and monotonic decay of the differential conductance signal (i.e., dI\/dV). In addition to this, we also inspected the differential conductance within \u00a0\u00b110\u2009mV to avoid strong variations around the Fermi energy.<\/p>\n For spectroscopy in the tunnelling regime with a normal metal tip, the differential conductance is given by<\/p>\n $$G(V)={\\sum}_{{{\\bf{k}}},\\sigma }{C}_{{{\\bf{k}}}\\sigma }\\int_{-\\infty }^{\\infty }\\frac{{{\\rm{d}}}E}{2\\pi }\\,{A}_{{{\\rm{s}}}}({{\\bf{k}}},\\sigma ;E)\\left(-\\frac{\\partial f(E+eV)}{\\partial E}\\right),$$<\/p>\n \n (29)\n <\/p>\n where As is the spectral function of the superconductor. The coupling Ck\u03c3 accounts for the spectral function of the tip and the tunnelling overlap between the tip\u2019s evanescent states and quasiparticle states from the superconductor with momentum k and spin \u03c3. We do not expect the coupling to vary much on any given Fermi surface \u03c0, and as such we approximate Ck,\u03c3\u00a0\u2248\u00a0C\u03c0. However, due to the different orbital composition and in-plane momenta of the quasiparticle states at the K, \\({K}^{{\\prime} }\\) valleys and at the \u0393 surface, cf. Fig.\u00a01<\/a>b, the coupling constant will be in general different. In particular, since the d2,0 orbital extends farthest out of plane and the STM favours states with small in-plane momentum63<\/a>,\u00a0as has been observed in previous STM measurements32<\/a>, we expect C\u0393 \u226b CK.<\/p>\n In the tunnelling regime, where the experimental data was taken, the lifetime of the quasiparticles is weakly affected by the coupling to the tip and the spectral function can be approximated with the local density of states of the quasiparticles, As(k,\u00a0\u03c3,\u00a0E)\u00a0=\u00a02\u03c0Dk,\u03c3(E). The tunnelling to each of the Fermi surfaces represents a distinct transport channel, whose strength is given by the sum over the local density of states<\/p>\n $${D}_{\\pi,\\sigma }(E)={\\sum}_{{{\\bf{k}}}\\in \\pi }{D}_{{{\\bf{k}}},\\sigma }(E).$$<\/p>\n \n (30)\n <\/p>\n The latter are modeled by assuming a BCS form factor modifying the local density of states \u03c1k,\u03c3 in the normal conducting state as<\/p>\n $${D}_{{{\\bf{k}}},\\sigma }(E)={\\rho }_{{{\\bf{k}}},\\sigma }{{\\rm{Re}}}\\sqrt{\\frac{{E}^{2}}{{E}^{2}-| {\\Delta }_{{{\\bf{k}}},\\sigma }{| }^{2}}}.$$<\/p>\n \n (31)\n <\/p>\n In the absence of magnetic fields, we can assume that, due to time reversal symmetry, D\u03c0,\u03c3\u00a0=\u00a0D\u03c0\/2 and C\u03c0,\u03c3\u00a0=\u00a0C\u03c0. As such, Eq. (29<\/a>) predicts for the total differential conductance a superposition of contributions from the distinct transport channels with weights provided by the couplings C\u03c0 as<\/p>\n $$G(V)\\approx {\\sum}_{\\pi }{C}_{\\pi }{G}_{\\pi }(V),$$<\/p>\n \n (32)\n <\/p>\n $${G}_{\\pi }(V)=\\int_{-\\infty }^{\\infty }{{\\rm{d}}}E\\,{D}_{\\pi }(E)\\left(-\\frac{\\partial f(E+eV)}{\\partial E}\\right).$$<\/p>\n \n (33)\n <\/p>\n For low temperature, the form of the gap on each separate Fermi surface results in characteristic signatures in the differential conductance of the respective transport channel. We use as a fit function this superposition of the differential conductance of the individual transport channels. For the gaps on the Fermi surfaces \u0394k,\u03c3, we use the form of the gaps as found from Eq. (24<\/a>) at the experimental temperature \\({T}_{\\exp }\\) indicated by stars in Fig.\u00a04<\/a>a. For single-band s-wave superconductors the ratio between the amplitude of the gap and the critical temperature is fixed as \\(\\Delta (T=0\\,{{\\rm{K}}})\/({k}_{{{\\rm{B}}}}{T}_{{{\\rm{c}}}})=\\pi {e}^{-{\\gamma }_{{{\\rm{e}}}}}\\approx 1.76\\), with Euler\u2019s constant \u03b3e1<\/a>. This ratio appears to be slightly violated in monolayer 1H-NbSe2 as the expected s-wave gap for the Tc\u00a0\u2248\u00a02\u2009K is \u0394(T\u00a0=\u00a00\u2009K)\u00a0\u2248\u00a00.3\u2009meV instead of the roughly 0.4\u22120.5\u2009meV gap we find in the spectral function, in agreement with earlier works33<\/a>. Such violations of the BCS-ratio are known to occur in multi-band or multi-patch superconductors and have been observed previously18<\/a>,70<\/a>. To account for this possibility, we allow for a fit of the amplitude of our solution by a common rescaling of all \u0394k,\u03c3 by a parameter A as \u0394k,\u03c3\u00a0\u2192\u00a0A\u0394k,\u03c3. The latter also takes care of the uncertainty of the actual experimentally realised Tc and the rescaling factor \u03b3 of the interaction.<\/p>\n To qualitatively account for additional sources of broadening in the experiment, we fit with an effective temperature Teff in the calculation of the derivative in Eq. (33<\/a>) which ranges between 1.3 and 1.4 times the recorded base temperature of the STM in Fig.\u00a04<\/a>d, e and is kept at 1.4 times the measurement temperature for Fig.\u00a04<\/a>f. To account for offsets in the calibration, we further allow for both a small constant offset G0 in the measured conductivity and V0 in the recorded voltage. The fit coefficients for the other traces shown in Fig.\u00a04<\/a>d, f only differ by small quantitative changes from the ones reported in Table\u00a01<\/a> for Fig.\u00a04<\/a>e. For the fitting we use the trust region reflective algorithm as implemented in SciPy\u2019s \u201ccurve_fit\u201d routine. We consider a range of \u00a0\u00b1\u00a00.5\u2009mV containing the main coherence peaks, but not the satellite features due to Leggett modes33<\/a>, which are not accounted for in our theory. The resulting fit parameters are listed in Table\u00a01<\/a>. The best fits for the solutions associated with the \\({E}^{{\\prime} }\\) irreducible representation are obtained by considering strongly selective coupling of the tip to the \u0393 pockets.<\/p>\n