Real-space anisotropy of the superconducting gap in the charge-density wave material 2H-NbSe2

We present a scanning tunneling microscopy (STM) and ab-initio study of the anisotropic superconductivity of 2H-NbSe2 in the charge-density-wave (CDW) phase. Differential-conductance spectra show a clear double-peak structure, which is well reproduced by density functional theory simulations enabling full k- and real-space resolution of the superconducting gap. The hollow-centered (HC) and chalcogen-centered (CC) CDW patterns observed in the experiment are mapped onto separate van der Waals layers with different electronic properties. We identify the CC layer as the high-gap region responsible for the main STM peak. Remarkably, this region belongs to the same Fermi surface sheet that is broken by the CDW gap opening. Simulations reveal a highly anisotropic distribution of the superconducting gap within single Fermi sheets, setting aside the proposed scenario of a two-gap superconductivity. Our results point to a spatially localized competition between superconductivity and CDW involving the HC regions of the crystal.


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 enhancement  . The layered transition-metal dichalcogenide 2H-NbSe 2 is a prominent example undergoing a transition into a chargeordered phase at T ≤33 K, and into the superconducting phase at T c ≤ 7.2 K. Both orderings are driven by electron-phonon coupling [23][24][25][26][27][28][29][30][31] with a strong momentum dependence 14,17,20,32 . The superconducting gap turns out to have an anisotropic distribution over the Fermi surface (FS), whereby NbSe 2 has been interpreted either as a multigap or a highly anisotropic superconductor [32][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 NbSe 2 . Knowing whether superconducting anisotropy favors the stability of the superconducting state (like, e.g., in MgB 2 ) 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 NbSe 2 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 superconductivity [37][38][39] , SCDFT does not rely on any empirical parameter, which allows us to elucidate the strength and distribution of the superconducting gap.

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 jq CDW jt2π= 3a ð Þ, 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.

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 q CDW 40 . As discussed in 23,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 dynamics 41 . 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 found that the CDW formation induces small-amplitude atomic displacements reducing the minimum Nb-Se distance bỹ 1.0%. Remarkably, the dynamically stable reconstruction that we predict has a unit cell with two crystallographically inequivalent NbSe 2 layers, L 1 and L 2 . 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.
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 (N F ) 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 unfolding 43,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 valenceband 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 k z = 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 .
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. L 1 states are highlighted in red on the left panel of the figure, while L 2 states are shown in blue on the right. We observe significant differences: (i) On the k z = 0 plane, the FS of L 2 develops ring structures, unlike L 1 . (ii) At k z ' 2 3 π c the inner Kcentered sheet is mainly composed by L 2 states, whereas it has a stronger L 1 character on the ΓKM plane. (iii) By projecting the Fermi-pocket resolved DOS onto the layers,  where P i nk is a projection operator with i = L 1 , L 2 , and D refers to Γor K-centered sheets, we found that, remarkably, N L1Γ F and N L1K F are equally affected by the CDW (being reduced by 20%), whereas the DOS of L 2 stays almost constant. The calculated values of N iD F are listed in Table 1.
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 findings [2][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   parameters. Layer projections of k-partial integrals of the Fermi DOS ( N iD F ), electron-phonon coupling (λ iD ) and their ratio (V iD ), over the FS pockets at the Γ and K points (i = L 1 , L 2 ; D indicates Γor K-centered sheets). properties of NbSe 2 , 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, g ν nkn 0 k 0 2 , is inversely proportional to the phonon frequency 37,46,47 , a dynamical instability of the lattice shows up in a divergence of the coupling. The calculated α 2 F(ω) for the dynamically stable 3 × 3 × 1 CDW structure (shown in Supplementary Note 1) gives a total EPC of λ = 1.4, placing NbSe 2 in the strong-coupling regime.
To assess the momentum anisotropy of the EPC, we calculated the k-resolved EPC parameters The values obtained for the unfolded k z = 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, and the corresponding DOS-independent parameters, The computed values of λ iD and V iD are collected in Table 1. Since V iD varies from a minimum of 0.30 (V L1Γ ) to a maximum of 0.39 (V L2K ), 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 L 2 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 SCDFT [48][49][50] with the state-of-the-art SPG functional 51 . The implementation of the method allowed us to make fully anisotropic calculations of the superconducting gap in both reciprocal and real space 52 . 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 T c of 10.2 K, is to be compared with an experimental value T exp c ¼ 7.2 K. Given the usual accuracy of SCDFT calculations (within~20% of T c ), 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 transition 41 , 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 μ * approximation 39,54 . This, in fact, gives T c ≃ 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 T c comes from including Coulomb interactions fully ab initio, thus retaining the energy dependence of the screened electron-electron interaction matrix elements V ε nk ; ε nk 0 ð Þ . In fact, our simulations show that the electronic states at the Fermi level are poorly screened, that is μ ¼ N F V 0; 0 ð Þ 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 MgB 2 , which has sheet-dependent anisotropy, the superconducting gap in momentum space of NbSe 2 , Δ nk , varies considerably in magnitude both between sheets and within a single sheet. We also note that whereas anisotropic effects in MgB 2 enhance T c 56-58 , that of NbSe 2 stays almost constant if anisotropy is washed out. Nevertheless, the degree of gap anisotropy is of essential importance for understanding thermodynamic and spectroscopic properties 27,32,34,36,45,[59][60][61] . Figure 3d shows the histograms of gap values over the full BZ (green) and the k z = 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 k z = 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 Δ = 1.85 meV, corresponding to 2Δ=k B T c = 4.3, which is in excellent agreement with experiments 59 . Our overestimation of T c naturally leads to an overestimation of the gap. However, multiplying Δ by the factor T exp c =T 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 measurements 33 . Figure 3c shows the unfolded gap structure on constant-k z slices of the BZ. We see a complex pattern of superconducting gaps that, like for other anisotropic superconductors [62][63][64] , mirrors the structure of the EPC (compare the k z = 0 slice with Fig. 3b). In the k z = 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 measurements 32 . 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 k z (see the slice at k z ¼ 2 3 π c ).

DISCUSSION
For a direct comparison with the experiment, we simulated STM images of bulk NbSe 2 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 differentialconductance spectra of layer L 1 , L 2 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 L 1 and L 2 (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, L 2 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 L 2 , whereas it almost disappears for L 1 . 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 A. Sanna et al. consider that the actual CDW is slightly incommensurate, and that the structural patterns on L 1 and L 2 exist, in fact, on the same surface, smoothly transforming into each other 42 . 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.
We mention that in a recent letter Heil and coworkers 65 used ab-initio Migdal-Eliashberg theory to simulate the superconducting properties and tunneling characteristics of NbS 2 . While the tunneling spectra show strong similarities with those we predict and measure for NbSe 2 , in the case of NbS 2 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 NbSe 2 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 MgB 2 , the gap anisotropy does not influence the stability of the superconducting state, that is perturbations affecting the gap anisotropy are not expected to reduce T c . (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 NbSe 2 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.

Crystal synthesis and characterization
Bulk crystals of NbSe 2 were grown by chemical vapor transport 32 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 DFT 66,67 using the local density approximation 68 . Atomic core states were treated in the normconserving pseudopotential approach as implemented in the Quantum Espresso code 69 . Phonon frequencies and EPC were computed by means of density functional perturbation theory 70 with an electronic temperature of 270 meV 40 . Further computational details are given in Supplementary Note 2. Fig. 5 Simulated real-space superconducting gap and STM spectra. a Simulated differential-conductance spectra of NbSe 2 for the bulk (green), layer L 1 (red), layer L 2 (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 exp c =T c to correct for the overestimation of T c . c Real-space view of the superconducting gap, showing stronger superconductivity on L 2 , especially at specific Nb sites (bright spots).
A. Sanna et al.