Introduction

Collective phenomena such as superconductivity and charge- or spin-density modulations have been found to co-exist in many materials, the interplay among these correlations being discussed in terms of their competition or mutual enhancement1,2,3,4,5,6,7,8,9,10,11,12,13,14,15,16,17,18,19,20,21,22. The layered transition-metal dichalcogenide 2H-NbSe2 is a prominent example undergoing a transition into a charge-ordered phase at T ≤33 K, and into the superconducting phase at Tc ≤ 7.2 K. Both orderings are driven by electron-phonon coupling23,24,25,26,27,28,29,30,31 with a strong momentum dependence14,17,20,32. The superconducting gap turns out to have an anisotropic distribution over the Fermi surface (FS), whereby NbSe2 has been interpreted either as a multigap or a highly anisotropic superconductor32,33,34,35,36.

A clear understanding of the structure of the superconducting gap is indeed presently missing, and is one of the key unsolved questions about NbSe2. Knowing whether superconducting anisotropy favors the stability of the superconducting state (like, e.g., in MgB2) is interesting on itself. More importantly, accessing the momentum-dependent structure of the superconducting gap in the CDW phase would allow one to unveil the link between superconducting condensation and CDW ordering in this material, resolving the long-standing debate on their competition.

In this paper, we present clear evidence of anisotropic superconductivity in NbSe2 by measuring the tunneling density of states (DOS) with unprecedented resolution combined with state-of-the-art first-principles simulations of the superconducting state based on density functional theory for superconductors (SCDFT). We emphasize that, unlike conventional approaches to superconductivity37,38,39, SCDFT does not rely on any empirical parameter, which allows us to elucidate the strength and distribution of the superconducting gap.

Results

STM spectroscopy

An atomically resolved constant-current STM image is shown in Fig. 1a. Besides the atomic corrugation of the terminating Se layer, the image reveals the modulation of the local DOS imposed by the CDW. As the CDW is incommensurate with the atomic lattice, (periodicity \(\vert{{{{{{\bf{q}}}}}}}_{\tiny{{{{\mbox{CDW}}}}}}\vert\lesssim 2\pi /\left(3a\right)\), where a is the in-plane lattice constant), different local patterns appear across the surface. The area in the orange rectangle shows a triangular pattern as a result of the maximum of the CDW being centered on a hollow site (hollow-centered, HC pattern). This continuously transforms into the chalcogen-centered (CC) pattern highlighted in the blue rectangle. Differential-conductance (dI/dV) spectra in different regions of the CDW are shown in Fig. 1b (for more spectra, see Supplementary Notes 4–6). All spectra exhibit a fully gapped region flanked by a broad peak distribution in the range 2.1 mV < Δ/e < 3.1 mV, with peculiar intensity variations, which can only be resolved by employing a superconducting tip.

Fig. 1: Scanning tunneling microscopy measurements.
figure 1

a Constant-current image of a clean NbSe2 surface (recorded with a Nb tip at a set point of 10 mV, 100 pA; the scale bar is 2 nm) to illustrate the smooth evolution of the CDW with respect to the atomic lattice. The frames highlight regions where the CDW maxima fall on a hollow site (orange) or on a Se atom (blue). b Differential-conductance dI/dV spectra taken at the positions indicated by stars in (a), recorded in constant-height mode (feedback was opened at the position of the blue spectrum for both traces with a set point of 5 mV, 250 pA and a modulation of 15 μVrms).

Full size image

Simulations - structural and normal-state properties

The CDW instability is captured by density functional theory (DFT) calculations in the undistorted cell (Ω1), where it manifests itself in imaginary phonon frequencies occurring over a broad range of momenta centered at qCDW40. As discussed in23,28,29, the CDW can not be viewed as a purely electronic Peierls’ instability driven by FS nesting, but is instead a structural phase transition assisted by the enhancement of the electron-phonon coupling (EPC) at the characteristic wave vector. Additionally, the CDW transition temperature can be accurately reproduced in DFT calculations by including anharmonic effects in the lattice dynamics41. In this paper we carry out a DFT study of the superconducting state in the CDW phase, providing a first-principles interpretation of the STM data.

We approximated the incommensurate CDW vector by \({{{{{{\bf{q}}}}}}}_{\tiny{{{{\mbox{CDW}}}}}}=\left(1/3,0,0\right)2\pi /a\), performing simulations within a superstructure unit cell (Ω3×3) 3 × 3 × 1 times Ω1. The CDW reconstruction was obtained by relaxing the lattice along the phonon instabilities until all the modes were stable in a 2 × 2 × 4 q-grid of the Ω3×3 cell.

We found that the CDW formation induces small-amplitude atomic displacements reducing the minimum Nb–Se distance by ~1.0%. Remarkably, the dynamically stable reconstruction that we predict has a unit cell with two crystallographically inequivalent NbSe2 layers, L1 and L2. These are shown in Fig. 2c,d and correspond, respectively, to the HC and CC patterns observed in experiment. We note that the obtained structures closely resemble those suggested in refs. 17,18,42. The calculated CDW stabilization energy is of 3.7 meV per formula unit.

Fig. 2: Electronic properties and CDW structure.
figure 2

a Unfolded band structure and (b) electronic DOS of the CDW (green) and undistorted phase (black). CDW reconstruction of layers (c) L2, chalcogen-centered pattern, and (d) L1, hollow-centered pattern, in the Ω3×3 unit cell (blue frame). The modulation of the atomic positions is highlighted by using a 2.53 Å cutoff for the plotted Nb–Se bond-sticks.

Full size image

The CDW distortion is expected to induce relatively small changes of the electronic structure. Figure 2a,b compares the computed band structure and electronic DOS in the presence of the CDW and in the undistorted system. The two DOSs differ mainly by a broad dip developing in the CDW state slightly above the Fermi level, such that the reduction in the DOS at the Fermi level (NF) induced by the CDW is relatively small (about 10%). Due to the large size of the supercell used in the calculations (54 atoms), the BZ folding makes the band structure in the CDW state hard to interpret. We thus employed a technique of band unfolding43,44 from the supercell BZ to the larger BZ of the Ω1 unit cell. As one can see, the CDW and undistorted bands largely overlap. A closer inspection, however, reveals band gap openings at several points of the BZ path, especially along the MK line, as well as energy shifts at the K point throughout the entire valence-band width.

The CDW effect on the electronic structure can be seen more clearly by considering a horizontal cut through the BZ. Fig. 3a shows the CDW reconstruction of the FS at kz = 0 obtained from the unfolding procedure. The unperturbed FS sheets (blue curves) are included for direct comparison. A double-walled cylindrical sheet of Nb 4d bands and a pocket of Se 4p character are centered at the Γ point. A second pair of Nb 4d-derived FS cylinders is located at the K point. This last double sheet appears to be significantly broader and blurred, indicating a stronger effect of the CDW distortion. Missing segments of the FS crossing the MK line, signal the openings of the CDW gap, which occur in excellent agreement with the ARPES measurements of refs. 32,45.

Fig. 3: Fermi surface anisotropy of the electron-phonon coupling and superconducting gap.
figure 3

a kz = 0 cut of the NbSe2 Fermi surface obtained by unfolding the BZ. Color scale indicates the relative spectral weight. Blue curves represent the unperturbed FS sheets. b Electron–phonon coupling and (c) low-temperature superconducting gap computed in the unfolded BZ at random points on constant-kz slices (of \(0.1\frac{\pi }{c}\) width) and within ±0.1 mRy from the Fermi level. d Distribution of gap values computed over the full BZ (green) and the \(| {{{{{{\bf{k}}}}}}}_{z}| < 0.1\frac{\pi }{c}\) slice (red).

Full size image

As already mentioned, the complexity of the CDW state is increased by the fact that the two structural layers undergo different reconstructions. To account for this aspect, we considered projections of the FS quasiparticle states onto the layers. Selected cuts of the unfolded projected FS are shown in Fig. 4. Most of the states (displayed in white) are found to span the entire crystal, yet a significant fraction is strongly localized on one layer. L1 states are highlighted in red on the left panel of the figure, while L2 states are shown in blue on the right. We observe significant differences: (i) On the kz = 0 plane, the FS of L2 develops ring structures, unlike L1. (ii) At \({{{{{{\bf{k}}}}}}}_{z}\simeq \frac{2}{3}\frac{\pi }{c}\) the inner K-centered sheet is mainly composed by L2 states, whereas it has a stronger L1 character on the ΓKM plane. (iii) By projecting the Fermi-pocket resolved DOS onto the layers,

$${N}_{{{{{{\rm{F}}}}}}}^{iD}=\mathop{\sum}\limits_{n{{{{{\bf{k}}}}}}\in D}\delta \left({\varepsilon }_{n{{{{{\bf{k}}}}}}}\right){P}_{n{{{{{\bf{k}}}}}}}^{i},$$
(1)

where \({P}_{n{{{{{\bf{k}}}}}}}^{i}\) is a projection operator with i = L1, L2, and D refers to Γ- or K-centered sheets, we found that, remarkably, \({N}_{{{{{{\rm{F}}}}}}}^{{{{{{{\rm{L}}}}}}}_{1}{{{{{\boldsymbol{\Gamma }}}}}}}\) and \({N}_{{{{{{\rm{F}}}}}}}^{{{{{{{\rm{L}}}}}}}_{1}{{{{{\bf{K}}}}}}}\) are equally affected by the CDW (being reduced by 20%), whereas the DOS of L2 stays almost constant. The calculated values of \({N}_{{{{{{\rm{F}}}}}}}^{iD}\) are listed in Table 1.

Fig. 4: Fermi surface projected on the structural layers.
figure 4

Slices of the unfolded Fermi surface computed at different values of kz (0.1\(\frac{\pi }{c}\) thick), from kz = 0 (top) to kz = π/c (bottom), showing the projection of the electronic states onto the two structural layers L1 and L2. States localized on L1 (L2) are highlighted in red (blue).

Full size image
Table 1 Layer and Fermi surface resolved electron–phonon parameters. Layer projections of k-partial integrals of the Fermi DOS (\({N}_{{{{{{\rm{F}}}}}}}^{iD}\)), electron–phonon coupling (λiD) and their ratio (ViD), over the FS pockets at the Γ and K points (i = L1, L2; D indicates Γ- or K-centered sheets).
Full size table

The CDW modulation can strongly affect the lattice dynamics and, in turn, influence phonon-driven superconductivity in complex ways. In principle, the CDW can co-exist with and even enhance superconductivity. This latter scenario, in particular, has been put forth by recent experimental findings2,3,4,5,6,7,11,14,15,16,19,21,22,32. As a first step to clarify the CDW effect on the superconducting properties of NbSe2, we simulated the lattice dynamics in the CDW phase. We point out that this is a computationally demanding task because of the large size of the crystallographic unit cell and the accuracy required in simulating the small CDW distortion.

We found that the phonon DOSs in the CDW and undistorted structure have a very similar shape. The unstable phonon mode, appearing with imaginary frequency in the latter, has, in fact, low spectral weight. Nevertheless, the CDW transition significantly alters the EPC: since the EPC strength, \({\left|{g}_{n{{{{{\bf{k}}}}}}n^{\prime} {{{{{\bf{k}}}}}}^{\prime} }^{\nu }\right|}^{2}\), is inversely proportional to the phonon frequency37,46,47, a dynamical instability of the lattice shows up in a divergence of the coupling. The calculated α2F(ω) for the dynamically stable 3 × 3 × 1 CDW structure (shown in Supplementary Note 1) gives a total EPC of λ = 1.4, placing NbSe2 in the strong-coupling regime.

To assess the momentum anisotropy of the EPC, we calculated the k-resolved EPC parameters

$${\lambda }_{n{{{{{\bf{k}}}}}}}=2\mathop{\sum}\limits_{n^{\prime} {{{{{\bf{q}}}}}},\nu }\frac{{\left|{g}_{n{{{{{\bf{k}}}}}}n^{\prime} {{{{{\bf{k}}}}}}+{{{{{\bf{q}}}}}}}^{\nu }\right|}^{2}}{{\omega }_{\nu {{{{{\bf{q}}}}}}}}\delta \left({\varepsilon }_{n^{\prime} {{{{{\bf{k}}}}}}+{{{{{\bf{q}}}}}}}\right).$$
(2)

The values obtained for the unfolded kz = 0 plane of the BZ are presented in Fig. 3b, where one can observe a complex pattern of hot spots (light) and weak-coupling (dark) regions. For a quantitative analysis, we examined the EPC contributions to Γ- and K-centered sheets arising from the two structural layers,

$${\lambda }^{iD}=\frac{1}{{N}_{{{{{{\rm{F}}}}}}}}\mathop{\sum}\limits_{n{{{{{\bf{k}}}}}}\in D}{\lambda }_{n{{{{{\bf{k}}}}}}}\delta \left({\varepsilon }_{n{{{{{\bf{k}}}}}}}\right){P}_{n{{{{{\bf{k}}}}}}}^{i},$$
(3)

and the corresponding DOS-independent parameters, \({V}^{iD}={\lambda }^{iD}/{N}_{{{{{{\rm{F}}}}}}}^{iD}\). The computed values of λiD and ViD are collected in Table 1. Since ViD varies from a minimum of 0.30 (\({V}^{{{{{{{\rm{L}}}}}}}_{1}{{{{{\boldsymbol{\Gamma }}}}}}}\)) to a maximum of 0.39 (\({V}^{{{{{{{\rm{L}}}}}}}_{2}{{{{{\bf{K}}}}}}}\)), we deduce that the layers contribute almost equally to the intrinsic coupling of the various FS pockets. λiD is then determined essentially by the value of the DOS at the Fermi level. We also note that the FS K-cylinder and the layer L2 can be identified as the regions with higher Fermi DOS and EPC in momentum and real space, respectively.

Simulation of the superconducting state

We simulated the superconducting state from first principles using SCDFT48,49,50 with the state-of-the-art SPG functional51. The implementation of the method allowed us to make fully anisotropic calculations of the superconducting gap in both reciprocal and real space52. The k-dependent SCDFT gap equation was solved on a mesh of random k points more densely distributed around the FS for higher accuracy (computational details are provided in Supplementary Note 2 and refs. 47,53). The estimated Tc of 10.2 K, is to be compared with an experimental value \({T}_{{{{{{\rm{c}}}}}}}^{\exp }=\) 7.2 K. Given the usual accuracy of SCDFT calculations (within ~20% of Tc), this error is likely to be ascribed to (i) the neglect of anharmonic corrections to the lattice dynamics, that are enhanced close to the CDW transition41, or (ii) limits on the precision in computing phonon properties and dielectric screening imposed by the large size of the CDW cell (see Supplementary Note 2 for an extended discussion). Nevertheless, our prediction is fairly accurate compared to that from conventional Eliashberg theory in the μ* approximation39,54. This, in fact, gives Tc ≃ 16 K under the usual assumption μ* = 0.11, while an exceedingly high value of the Coulomb pseudopotential (μ* = 0.28) is required to reproduce the experimental result. The improved SCDFT estimation of Tc comes from including Coulomb interactions fully ab initio, thus retaining the energy dependence of the screened electron-electron interaction matrix elements \(V\left({\varepsilon }_{n{{{{{\bf{k}}}}}}},{\varepsilon }_{n{{{{{\bf{k}}}}}}^{\prime} }\right)\). In fact, our simulations show that the electronic states at the Fermi level are poorly screened, that is \(\mu ={N}_{{{{{{\rm{F}}}}}}}V\left(0,0\right)\) is large. Moreover, because of a different orbital character, these states weakly interact with those far from the Fermi level, yielding a weak reduction of the Coulomb pseudopotential (high μ*)55.

The computed superconducting state turns out to be extremely complex. Unlike MgB2, which has sheet-dependent anisotropy, the superconducting gap in momentum space of NbSe2, Δnk, varies considerably in magnitude both between sheets and within a single sheet. We also note that whereas anisotropic effects in MgB2 enhance Tc56,57,58, that of NbSe2 stays almost constant if anisotropy is washed out. Nevertheless, the degree of gap anisotropy is of essential importance for understanding thermodynamic and spectroscopic properties27,32,34,36,45,59,60,61. Figure 3d shows the histograms of gap values

$${{{{{\rm{d}}}}}}\left({{\Delta }}\right)=\frac{1}{{N}_{{{{{{\rm{F}}}}}}}}\mathop{\sum}\limits_{n,{{{{{\bf{k}}}}}}}\delta ({{{\Delta }}}_{n{{{{{\bf{k}}}}}}}-{{\Delta }})\delta ({\varepsilon }_{n{{{{{\bf{k}}}}}}})$$
(4)

over the full BZ (green) and the kz = 0 cut of the FS (red). The total distribution reveals a spread of gap values with the main peak at 2.3 meV and several smaller structures at lower values of Δ. While mid- and low-energy features of the curve are well reproduced by the kz = 0 slice of the CDW BZ, the source of the peak can be identified with a hot spot in Δnk at 2/3 of the ML line. The obtained average gap value is of \(\bar{{{\Delta }}}\) = 1.85 meV, corresponding to \(2\bar{{{\Delta }}}/{k}_{{{{{{\rm{B}}}}}}}{T}_{{{{{{\rm{c}}}}}}}\) = 4.3, which is in excellent agreement with experiments59. Our overestimation of Tc naturally leads to an overestimation of the gap. However, multiplying Δ by the factor \({T}_{{{{{{\rm{c}}}}}}}^{\exp }/{T}_{{{{{{\rm{c}}}}}}}\), the estimated gap is rescaled in between 0.7 and 1.7 meV, which is close to the range seen in Fig. 1b and provided by ARPES measurements33. Figure 3c shows the unfolded gap structure on constant-kz slices of the BZ. We see a complex pattern of superconducting gaps that, like for other anisotropic superconductors62,63,64, mirrors the structure of the EPC (compare the kz = 0 slice with Fig. 3b). In the kz = 0 plane, a high-gap region can be identified with the inner FS sheet of Nb character around the Γ point, while the Fermi arcs at K have, on average, a much smaller superconducting gap. This observation agrees with the inverse correlation between superconducting and CDW gap suggested by ARPES measurements32. Nevertheless, as already mentioned, the FS K-cylinder largely contributes to the superconducting condensation energy due to the overall higher DOS (Table 1) and a significant increase in the gap at large kz (see the slice at \({{{{{{\bf{k}}}}}}}_{z}=\frac{2}{3}\frac{\pi }{c}\)).

Discussion

For a direct comparison with the experiment, we simulated STM images of bulk NbSe2 by a superconducting Nb tip using the SCDFT gap. Specifically, the superconducting DOS of the sample was computed locally at the tip position from the real-space projection of Eq. (4) (see Supplementary Note 3 for an extended discussion). Figure 5a shows the computational differential-conductance spectra of layer L1, L2 and the full bulk, compared to the experimental spectrum averaged over the sample surface. The tunneling conductance calculated at the points corresponding to the experimental positions in Fig. 1a is also included. The related histograms of gap values are presented in Fig. 5b. Additionally, Fig. 5c provides a real-space image of the superconducting gap from a side view (bottom) and through horizontal cuts on the Nb planes of L1 and L2 (top). Simulations show a strong variation of the tunneling conductance with the scanning position, especially between the layers. Due to the higher DOS at the Fermi level, L2 has, in fact, a much larger average gap (see the high-energy peak in the distribution of gap values). Moreover, a double-peak structure is clearly visible in the differential conductance computed for L2, whereas it almost disappears for L1. Changes in the experimental spectra as a function of the tip position are much less pronounced, and two distinct peaks are always observed (see Fig. 1 and Supplementary Note 4). However, in comparing simulated and experimental spectra, one should consider that the actual CDW is slightly incommensurate, and that the structural patterns on L1 and L2 exist, in fact, on the same surface, smoothly transforming into each other42. This causes an additional degree of hybridization between the electronic states of the two layers, which is not accounted for in the simulations, as it would require using a prohibitively large cell. Additionally, simulations are done in the bulk, thus in absence of surface effects which may influence the gap anisotropy. Nevertheless, by comparing the black and green curves in Fig. 5, we see that the shape of the measured dI/dV tunneling spectrum is, overall, quite well reproduced, especially concerning the positions of the two peaks and the shoulder-like feature at ≈ 3 meV.

Fig. 5: Simulated real-space superconducting gap and STM spectra.
figure 5

a Simulated differential-conductance spectra of NbSe2 for the bulk (green), layer L1 (red), layer L2 (dark blue) and at the experimental positions of Fig. 1a (orange, light blue). Experimental spectrum (black) obtained by averaging over multiple data traces as in Supplementary Fig. 4b. Tunneling currents are computed from the real-space distributions of gap values shown in (b), where Δ is scaled by the factor \({T}_{{{{{{\rm{c}}}}}}}^{\exp }/{T}_{{{{{{\rm{c}}}}}}}\) to correct for the overestimation of Tc. c Real-space view of the superconducting gap, showing stronger superconductivity on L2, especially at specific Nb sites (bright spots).

Full size image

We mention that in a recent letter Heil and coworkers65 used ab-initio Migdal–Eliashberg theory to simulate the superconducting properties and tunneling characteristics of NbS2. While the tunneling spectra show strong similarities with those we predict and measure for NbSe2, in the case of NbS2 they are found to originate from a clear two-gap structure. If, on the one hand, these results demonstrate the importance of ab-initio methods for the interpretation of tunneling data, they also show that similar dichalcogenides can hide significantly different superconducting properties. It would be of interest to address this aspect by a detailed comparative ab-initio study.

To summarize, we presented high-resolution STM measurements of the differential conductance of NbSe2 and an ab initio fully anisotropic description of the CDW and superconducting state. The agreement between measured and computed superconducting properties appears to be good. From our analysis emerge the following key aspects: (i) Contrary to the textbook case of MgB2, the gap anisotropy does not influence the stability of the superconducting state, that is perturbations affecting the gap anisotropy are not expected to reduce Tc. (ii) The gap distribution on the FS shows rapid variations with k within each sheet. However, a high-gap structure can be clearly identified as responsible for the main STM coherence peak. Specifically, the peak originates from the same FS sheet that is broken by the occurrence of the CDW, i.e., the FS cylinder centered at the K point. (iii) HC regions of the crystal undergo a strong reduction of the Fermi DOS by the CDW, leading to a lower EPC compared to CC regions (whose Fermi DOS is almost unaffected).

Overall, our simulations indicate that superconducting and CDW orderings in NbSe2 compete in HC regions, whereas superconductivity is not suppressed in CC regions. We rule out, however, a two-gap scenario, providing evidence of strong superconducting gap anisotropy.

Methods

Crystal synthesis and characterization

Bulk crystals of NbSe2 were grown by chemical vapor transport32 and cleaved under ultra-high vacuum. Scanning tunneling microscopy and spectroscopy measurements were performed at 1.1 K with superconducting tips fabricated by indenting a clean NbTi wire into a Nb(111) crystal until the tip’s energy gap reached the bulk value (2Δtip ≈ 3.1 meV). Since the probed signal is a convolution of the DOS of substrate and tip, all the measured spectral features of the sample are shifted by the tip’s energy gap. Additionally, the differential-conductance signal is weighted by the energy-dependent tunneling probability.

Theoretical calculations

We calculated normal-state properties within DFT66,67 using the local density approximation68. Atomic core states were treated in the norm-conserving pseudopotential approach as implemented in the Quantum Espresso code69. Phonon frequencies and EPC were computed by means of density functional perturbation theory70 with an electronic temperature of 270 meV40. Further computational details are given in Supplementary Note 2.