Spectroscopy of bulk and few-layer superconducting NbSe2 with van der Waals tunnel junctions

Tunnel junctions, an established platform for high resolution spectroscopy of superconductors, require defect-free insulating barriers; however, oxides, the most common barrier, can only grow on a limited selection of materials. We show that van der Waals tunnel barriers, fabricated by exfoliation and transfer of layered semiconductors, sustain stable currents with strong suppression of sub-gap tunneling. This allows us to measure the spectra of bulk (20 nm) and ultrathin (3- and 4-layer) NbSe2 devices at 70 mK. These exhibit two distinct superconducting gaps, the larger of which decreases monotonically with thickness and critical temperature. The spectra are analyzed using a two-band model incorporating depairing. In the bulk, the smaller gap exhibits strong depairing in in-plane magnetic fields, consistent with high out-of-plane Fermi velocity. In the few-layer devices, the large gap exhibits negligible depairing, consistent with out-of-plane spin locking due to Ising spin–orbit coupling. In the 3-layer device, the large gap persists beyond the Pauli limit.

S uperconductors of the transition metal dichalcogenide (TMD) family have seen a revival of interest subsequent to developments in device fabrication by mechanical exfoliation [1][2][3][4][5][6][7] . The isolation of ultrathin NbSe 2 1,2 has yielded indications of a Berezinskii-Kosterlitz-Thouless transition, which occurs in 2D superconductors. NbSe 2 1 , gated MoS 2 5,6 , and gated WS 2 7 devices also remain superconducting in in-plane magnetic fields well beyond the Pauli limit H p ¼ Δ= ffiffi ffi g p μ B À Á 8,9 . (Here, Δ is the superconducting energy gap, g is the Landé g-factor and μ B the Bohr magneton.) This is likely due to Ising spin-orbit coupling (ISOC): The broken inversion symmetry of the monolayer TMD in the plane is expected to lead to the formation of Cooper pairs whose constituent spins are locked in the out-of-plane direction, in a singlet configuration. Interestingly, zero-resistance states have been observed in parallel magnetic fields exceeding the Pauli limit even in few-layer devices, where inversion symmetry is recovered 1 . This suggests that the inter-layer coupling is not strong enough to overcome the out-of-plane spin-locking due to ISOC, perhaps in part due to the presence of spin-layer locking 10 . These previous studies 1-7 used in-plane electronic transport at high magnetic field and temperatures close to the critical temperature, T c (H = 0), to determine the upper critical field H c2|| , which depends on the magnitude of the spin-orbit field H SO . Tunneling spectroscopy can provide complementary information. For example, tunneling and other probes (see ref. 11 and references therein) suggest that bulk NbSe 2 is a 2-band superconductor. The role of both bands in ultra-thin devices, and their response to magnetic field, can be addressed by tunnel spectroscopy but not by electronic transport (which addresses the condensate) or heat conductivity measurements 12 . The latter, specifically, would be shorted by the substrate. Tunneling can also probe the effect of magnetic field at the full temperature range and in fields ranging from zero to H c2|| .
To carry out tunneling measurements on devices of variable thickness, at variable magnetic field conditions, it is necessary to develop a device-based architecture suitable for integration with TMDs. Oxide-based tunnel barriers, such as those used since the days of Giaever 13 , have now reached technological maturity; however, there is a limited number of oxides which form highquality insulating, non-magnetic barriers, and they do not grow well on all surfaces. It is therefore of interest to explore alternatives based on van der Waals (vdW) materials 14 , ultrathin layers of which can be precisely positioned on many surfaces 15 . Indeed, such barriers have proven effective when integrated with graphene [16][17][18] and appear to be promising candidates for integration with TMD superconductors 19 .
In this work we study tunneling devices utilizing ultrathin TMDs as barriers. We find that these devices exhibit a hard gap with strong sub-gap conductance suppression. In the spectra of bulk NbSe 2 , we clearly resolve the presence of two superconducting order parameters, which exhibit a distinct evolution upon application of in-plane magnetic fields. Spectra of ultrathin NbSe 2 are stable against the application of magnetic fields, even beyond the Pauli limit 8,9 .

Results
Hard gap van-der-Waals tunnel junctions. We fabricate Normal-Insulator-Superconductor (NIS) tunnel junctions with either MoS 2 or WSe 2 -both van der Waals (vdW) materials-as the insulating barrier. The barrier material is placed on top of 2H-NbSe 2 (hereafter NbSe 2 ), a vdW superconductor with bulk T c ≈ 7.2 K, using the 'dry transfer' fabrication technique 14,15 . Crucially, such heterostructures can comprise NbSe 2 flakes of varying thicknesses, from bulk (≳6 layers) to few-layer, often within a single device. Figure 1a shows a typical junction consisting of a 20 nm-thick NbSe 2 flake partially covered by a 4-5 layer thick MoS 2 barrier (cf. Supplementary Note 1 and Supplementary  Fig. 1). The junction has an area A = 1.6 μm 2 , and we evaluate its transparency to be T $ 10 À8 (Supplementary Note 2 and Supplementary Fig. 2). Figure 1b shows the differential conductance G = dI/dV as a function of V obtained with the device, at T = 70 mK. This spectrum has two striking features: first, the very low sub-gap conductance (G 0 R N ≈ 1/500), which indicates all bands crossing the Fermi energy are fully gapped. Second, the intricate structure of the quasiparticle peak differs from a standard BCS density-of-states (DOS) by having a relatively low peak and a shoulder at lower energies. The latter feature can be clearly seen in the second derivative (Fig. 1c) where the slope separates into a double peak feature, similar to STS scans of bulk NbSe 2 11,20 . Based on this, the flake can be considered bulk in terms of the zero field superconducting properties, and is hence referred to as the bulk sample. Two-band model. Density functional theory calculations 21 , and ARPES data 22 show that five independent electronic bands cross the Fermi energy in NbSe 2 . Of these, four are Nb-derived bands with roughly cylindrical Fermi surfaces, centered at the Γ and K points. The fifth is derived from the Se p z orbitals, which give rise to a small ellipsoid pocket around the Γ point. Ref. 11 uses a twoband model to fit NbSe 2 tunneling data, which can be justified noting that the Se and Nb-derived bands differ in the density of states and value of the electron-phonon coupling parameter 23 . We follow ref. 11 in fitting our data using the same two-band model, which was developed in various forms by Suhl 24 , Schopohl 25 , and McMillan 26 (below "SSM").
The model entails a self-consistent solution to the coupled equations for the energy dependent order parameters Δ i (E) in the two bands i: Δ 0 i describes the intrinsic gap within each band i, that is generated by the electron-phonon coupling and by the scattering rates of quasiparticles between the bands Γ ij . The extension of the two-band model to include Abrikosov-Gor'kov depairing 27-30via the terms with Γ AG i -was done by Kaiser and Zuckermann 31 . Here, depairing is due to magnetic field; thus, Γ AG i are set to 0 when no magnetic field is applied. The DOS of each band is then given by where N i (E F ) is the DOS at the Fermi energy in the normal state in band i. The parameter α = 0.1, incorporates band-anisotropy. This anisotropy also affects the effective fit temperature (T = 0.44 K), which is higher than the sample temperature. Our fit indicates the presence of two independent order parameters. The larger, Δ 0 1 ¼ 1:26 ± 0:01 meV can be determined with high fidelity. The smaller order parameter Δ 0 2 , cannot be determined unambiguously, and could have any value between 0 and 0.3 meV. The suppressed intrinsic superconductivity in the second band is consistent with a band with small density of states and weak electron-phonon coupling. The Se-band, having these properties 23 , is thus a candidate for the 2nd band.
NbSe 2 in in-plane magnetic field: diffusive vs. ballistic regimes. We next investigate the evolution of the tunneling spectra vs. inplane magnetic field H || (model- Fig. 1d, data-Fig. 1e). The bulk sample appearing in Fig. 1 is thick enough to accommodate Meissner currents leading to orbital depairing. We track the peaks in d 2 I/dV 2 , which are associated with the two gaps. The high energy peak in d 2 I/dV 2 depends weakly on the field, whereas the low energy peak evolves nonlinearly towards lower energies. This trend persists up to H = 0.5 T which we interpret as H c1 . Bulk NbSe 2 is regarded as a clean-limit superconductor (coherence length ξ < l mean free path). This condition, however, might not hold equally for both bands, since the smaller gap is associated with a larger ξ. Indeed, the evolution of the low energy gap in H || does not agree with the relevant, Doppler shift model [32][33][34] . Instead, the magnetic field evolution of its d 2 I/dV 2 features towards lower energies is well-reproduced by the diffusive Kaiser-Zuckermann (KZ) model (Fig. 1d), assuming a depairing parameter Γ AG 2 quadratic in H || . For a thin sample of thickness d ≪ λ, hc 235 , with D i the diffusion coefficient, d ≈ 20 nm and the penetration length λ ≈ 230 nm 36 . We find Γ AG 2 = 640 μeV at 1 T, corresponding to D 2 = 40 cm 2 s −1 . The model therefore yields a very large value for the diffusion coefficient D 2 = v F l/3, indicating a large Fermi velocity. This lends further support to the identification of this feature with the Se-derived band, where v F is 4-5 times larger than in the Nb-derived bands 23 . In contrast, the high energy quasiparticle peak moves only slightly and linearly to higher energies at low H || (dashed line in Fig. 1e). From this, we estimate v F ≈ 5 × 10 4 m s −1 . This indicates that a comprehensive description of the field evolution of the full spectrum of bulk NbSe 2 would require a model with arbitrary disorder bridging clean and diffusive limits.
Reduced gap in ultra-thin NbSe 2 . Our methods allow us to carry out a straightforward comparison of the tunneling spectra from devices with different NbSe 2 thicknesses. Figure 2a shows differential conductance curves taken by tunneling into ultra-thin NbSe 2 flakes (3 layers, 4 layers), in comparison to the dI/dV of the bulk flake discussed above. The spectra of the thin devices are in good agreement with the SSM model (fits are shown in Supplementary Fig. 3, fit parameters are given in Supplementary  Table 1). We extract the values of Δ 0 1 from these fits, and  separately evaluate T c using the temperature dependence of the tunneling conductance (details in Supplementary Fig. 4). We find that Δ 0 1 increases with T c (Fig. 2b), however, deviating slightly from the BCS result Δ = 1.76k B T c . This deviation is plausible due to the multi-band nature of superconductivity in NbSe 2 . The T c values measured here are lower than those reported elsewhere 1,37 , likely due to higher disorder in our sample [38][39][40] . Note that, based on ref. 37 , it would seem that the T c dependence on number of layers seen by all works to date is not due to strain or other substrate effects.
Ultra-thin NbSe 2 in magnetic field. We now turn to the response of the ultrathin devices to in-plane magnetic fields. Figure 2 shows the tunneling spectra for the bulk device (Fig. 2a) and for the 4-layer and 3-layer devices (Fig. 2b and c, respectively). The KZ model, for all 3 devices, is shown in Fig. 2d-f, respectively. In the thin devices, unlike the bulk device, we find that the spectrum changes very little up to 3.5 T (which is the maximal field where a compensation coil keeps a zero perpendicular field). In both 3and 4-layer devices, there is a small reduction in the height of the quasiparticle peak; in the 3-layer device, the low energy spectrum exhibits a more intricate evolution (discussed below). Thinner samples have shorter scattering lengths and longer coherence lengths, and should thus be closer to the diffusive limit. Indeed, the KZ model fits our data well, and we use it to quantify the reduction of the peak height of the 3-and 4-layer devices (Fig. 2e,  f), finding that the depairing term Γ AG 1 % 0:5 μeV at 1 T (all model parameters are given in Supplementary Table 2).
Since orbital depairing is quadratic in sample thickness, we expect it to be diminished in the 3-and 4-layer devices, allowing us to probe the spin-dependent interaction. The interaction of the spin with magnetic field should lead to Zeeman splitting of the spectrum, which we do not see. This could be due to two mechanisms. First, spin-orbit scattering can effectively randomize the spin, giving rise to a depairing parameter given by where τ SO is the spin orbit scattering time 35 . Second, ISOC can align the spins in the out-of-plane direction with an effective field, H SO , and the depairing term Γ AG $ 2μ B H 2 jj =H SO 1,5 . The first scenario can be ruled out, since it yields τ SO < 50 fs, shorter than the scattering time 5 . The ISOC case is more likely, and the depairing energy of 0.5 μeV at 1 T (for the 3L device) yields H SO ≈ 200 T. Using H p = 4.9 T (extracted by setting Δ 1 = 0.4 meV) we can estimate H c2 ¼ ffiffiffiffiffiffiffiffiffiffiffiffiffi ffi H SO H p p ≈ 30 T, consistent with transport experiments 1 . We note that since the orbital term is not entirely suppressed, and we cannot estimate its contribution to the depairing, the estimate for H SO is a lower bound. Further details concerning the possible interpretation of the depairing term are given in Supplementary Note 3 and Supplementary Table 2.
The stability of the larger gap can even be demonstrated above the Pauli limit, H p = 4.9 T for the 3-layer device. In Fig. 3c we present the density of states of this sample above the Pauli limit by applying an in-plane field of 6.4 T. At this field, our measurement system did not allow us to compensate for angle misalignment leading to small component of perpendicular field (≈0.2 T) and possible vortex penetration. Nevertheless, the size of the gap remains unchanged. This lends further support to ISOC as the mechanism protecting NbSe 2 superconductivity at high parallel fields.
Anomalous secondary gap. The sub-gap spectrum of the 3-layer device appears to exhibit a secondary, well-formed small gap of Δ 2 = 50 μV, which is suppressed at H || = 1.2 T. As we show in Fig. 3c, f, the KZ model reproduces this data remarkably well. Here, too, the depairing term is quadratic in H || , with Γ AG 2 ≈ 13 μeV at 1 T. This value is too big to be interpreted in terms of orbital depairing. The observed depairing could alternatively be associated with spin-orbit-driven spin-flip scattering, with τ SO ≈ 1.5 ps. This depairing shows that unlike the band with larger gap, this band is not protected by ISOC. This difference can be explained by associating the smaller order parameter with the Se-derived band. The outer gap, appearing immune to depairing, would then be associated with the Nb-derived K-band. We note that this interpretation leaves open the question of the role of the Nbderived Γ-band. Addressing this question will require further spectral studies-in particular, of monolayer NbSe 2 where the Sederived band does not cross the Fermi energy.

Discussion
Our results show that TMD semiconductors transferred on top of NbSe 2 form stable tunnel barriers with a hard gap. We show that the SSM model, modified to include diffusive depairing, successfully reproduces the tunneling spectra in both the bulk and the ultrathin limits, in the presence of in-plane magnetic fields. This allows us to probe the effect of the spin and orbital degrees of freedom on the spectra, thereby differentiating between the responses of the different bands to the field. The large gaps, in the 3-and 4-layer devices, are remarkably stable to depairing by the in-plane field, exhibiting very small depairing energies (<1 μeV), which place a tight cap on the spin-depairing observed on this band, lending support to ISOC as the mechanism behind this stability. We suggest that our technique can be generalized to work with many other material systems, such as organic (super) conductors and other fragile systems which have hitherto not been investigated using tunneling spectroscopy.

Methods
Device fabrication and measurement. The vdW tunnel junctions were fabricated using the dry transfer technique 41 , carried out in a glove-box (nitrogen atmosphere). NbSe 2 flakes were cleaved using the scotch tape method, exfoliated on commercially available Gelfilm from Gelpack, and subsequently transferred to a SiO 2 substrate. MoS 2 and WSe 2 flakes were similarly exfoliated and thin flakes suitable for the formation of tunnel barriers were selected based on optical transparency. The barrier flake was then transferred and positioned on top of the NbSe 2 flake at room temperature. Ti/Au contacts and tunnel electrodes were fabricated using standard e-beam techniques. Prior to the evaporation of the ohmic contacts the sample was ion milled for 15 s. No such treatment was done with the evaporation of the tunnel electrodes. All transport measurements were done in a 3 He-4 He dilution refrigerator with a base temperature of 70 mK. The AC excitation voltage was modulated at 17 Hz; its amplitude was 15 μV at all temperatures for the bulk device and 10 μV for the few-layer devices. Measurement circuit details are provided in Supplementary Fig. 5.
Data availability. The data that support the findings of this study are available from the corresponding author upon request.