Tilted vortex cores and superconducting gap anisotropy in 2H-NbSe2

Superconducting vortex cores have been extensively studied for magnetic fields applied perpendicular to the surface by mapping the density of states (DOS) through Scanning Tunneling Microscopy (STM). Vortex core shapes are often linked to the superconducting gap anisotropy---quasiparticle states inside vortex cores extend along directions where the superconducting gap is smallest. The superconductor 2H-NbSe$_2$ crystallizes in a hexagonal structure and vortices give DOS maps with a sixfold star shape for magnetic fields perpendicular to the surface and the hexagonal plane. This has been associated to a hexagonal gap anisotropy located on quasi two-dimensional Fermi surface tubes oriented along the $c$ axis. The gap anisotropy in another, three-dimensional, pocket is unknown. However, the latter dominates the STM tunneling conductance. Here we measure DOS in magnetic fields parallel to the surface and perpendicular to the $c$ axis. We find patterns of stripes due to in-plane vortex cores running nearly parallel to the surface. The patterns change with the in-plane direction of the magnetic field, suggesting that the sixfold gap anisotropy is present over the whole Fermi surface. Due to a slight misalignment between the vector of the magnetic field and the surface, our images also show outgoing vortices. Their shape is successfully compared to detailed calculations of vortex cores in tilted fields. Their features merge with the patterns due to in plane vortices, suggesting that they exit at an angle with the surface. Measuring the DOS of vortex cores in highly tilted magnetic fields with STM can thus be used to study the superconducting gap structure.

T he superconducting compound 2H-NbSe 2 has been considered for a long time as a prototypical example of a material with anisotropic superconducting propertiesfor instance, the upper critical field is three times larger when the field is applied in-plane than when it is perpendicular to the planes of the hexagonal crystalline structure 1 . The Fermi surface of 2H-NbSe 2 has two tubes around the Γ and K points due to bands derived from Nb orbitals which are nearly two-dimensional and a small three-dimensional pocket centered on the Γ point due to bands derived from Se orbitals [1][2][3][4][5] . A charge density wave (CDW) opens below 30 K in 2H-NbSe 2 , coexisting with superconductivity below T c = 7.2 K 4 . Scanning Tunneling Microscopy (STM) measurements in magnetic fields perpendicular to the hexagonal planes provide Density of States (DOS) maps giving sixfold vortex cores and the atomic scale tunneling conductance has a sixfold modulation [6][7][8][9][10] . Photoemission experiments suggest that the CDW induces the sixfold anisotropic superconducting gap 3,4 . This is reinforced by STM measurements on the superconducting gap in the isostructural compound 2H-NbS 2 , which has a similar T c and no CDW and is isotropic in-plane 11 .
Photoemission experiments also suggest that the superconducting gap anisotropy is located in the Nb tubes and do not address the three dimensional Se-derived band 3,4 . On the other hand, penetration depth studies show that the gap at the Se sheet is large 12 . Furthermore, the Se sheet plays an essential role in most tunneling experiments 5 . Tunneling occurs preferentially through the last layer that consists of the hexagonal Se atom lattice 5,13 . It is however yet unclear what is the contribution of the electronic properties derived from Se to the superconducting gap anisotropy. To study this issue, measurements of the vortex core shape on a surface perpendicular to the c axis would be useful. However, no vortex imaging can then be made, because the surface is highly irregular due to the sheet-like structure of samples of 2H-NbSe 2 . We study in this work the usual surface of 2H-NbSe 2 , but we apply tilted magnetic fields.
Hess et al. measured the vortex lattice of 2H-NbSe 2 with STM in tilted magnetic fields 14,15 . They focused on the structure of the vortex lattice and found a distortion of the hexagonal lattice compatible with the anisotropy of the upper critical field, as well as a rotation of the orientation of the vortex lattice that is consistent with the anisotropic London theory [16][17][18] . When the magnetic field was close to being parallel to the hexagonal planes, they observed elongated vortex cores, instead of the sixfold cores seen in perpendicular fields, and a peculiar pattern of stripes. More recently, this pattern of stripes was addressed in STM measurements with the field exactly in-plane focusing on the influence of the structure of the vortex lattice and of screening currents 19,20 . Until now, there is no systematic study of the vortex core shape with the in-plane direction of the magnetic field and the relation of the observed stripe pattern with vortex core states remains unclear.
Here we visualize the DOS patterns and perform a detailed study as a function of the direction and bias voltage and show, thanks to a slight misorientation of the sample when the magnetic field is parallel to the surface, that the pattern of stripes is due to bound states of subsurface vortices. Furthermore, we find that the pattern changes between a two-fold vs. three-fold structure depending on the direction of the in-plane magnetic field with respect to the crystalline direction in the hexagonal plane. The images strongly change with the bias voltage, with the stripes disappearing when reaching the gap edge. We perform microscopic calculations of the DOS maps and their bias voltage dependence and find that the shape of in-plane vortex cores depends on the direction of the magnetic field with respect to the in-plane crystal lattice direction.

Results
Vortex lattice in tilted magnetic fields. In Fig. 1a we show the vortex lattice at 0.6 T with varying polar angle θ. The magnetic field is tilted at a fixed azimuthal angle φ = 0°, along an in-plane crystalline axis. That is, along the nearest neighbor axis in the hexagonal plane of the atomic surface plane (given by Se atoms, see Supplementary Note 2 and refs. 8,10 ). For all tilts we observe a peak in the tunneling conductance due to Caroli-de Gennes-Matricon bound states at the center of the vortex cores, similar to the one observed at perpendicular magnetic fields [6][7][8]10 . However, there are remarkable differences in tilted fields.
The vortex core as observed in the zero bias tunneling conductance maps continuously increases its size when tilting the magnetic field (see also Supplementary Note 2) and acquires a two-fold shape at nearly in-plane magnetic fields. For these high tilt angles, we observe a pattern of stripes in the zero bias tunneling conductance maps (Fig. 1b).
In Fig. 1c we show the bias-voltage dependence of the tunneling conductance when crossing a stripe. We start from a gapped behavior in-between vortices and observe a small zerobias peak when we cross the stripe. The zero-bias peak splits into two peaks at nonzero bias when leaving the center of the stripe (Fig. 1c inset, see also Supplementary Fig. 3), as expected for Caroli-de Gennes-Matricon bound states [6][7][8]10 . Thus, the stripes result from vortex-core bound states from subsurface vortices oriented along the tilt direction. Accordingly, the stripes always follow the direction of the magnetic field.
In addition, the pattern of stripes changes with the in-plane azimuthal angle φ with respect to the in-plane hexagonal crystal lattice. For the field along a crystal axis, φ = 0°(left panel of Fig. 1b), the amount of stripes is about half the amount of stripes observed for φ = 30°(right panel of Fig. 1b, this is at 90°(or equivalently at 30°) from the nearest neighbor Se-Se direction of the surface atomic lattice, see Supplementary Note 3 and refs 8,10 ). The distance between vortices in the bulk is about 90 nm for a parallel magnetic field of 0.6 T. In Fig. 1b (left panel), when the magnetic field is along a crystalline axis (φ = 0°), we obtain 90 nm for the average distance between stripes. On the other hand, when the tilt of the magnetic field is in between crystalline axis (φ = 30°), we observe an average distance between stripes of 40 nm, about half the value found for φ = 0° (Fig. 1b, right panel). This suggests that there is one stripe per in-plane vortex for φ = 0°, but two stripes for φ = 30°per in-plane vortex. We can understand this if we consider that the in-plane vortex has a star shape, whose orientation changes with φ. As we will see below, a sixfold star vortex core whose shape is locked to the crystal lattice can provide such a behavior.
Vortex core in tilted magnetic fields. In Fig. 2 we compare the bias voltage dependence in perpendicular and tilted fields. In perpendicular magnetic fields, when increasing the bias voltage, we observe the same behavior as reported previously, showing the star shaped DOS associated to the in-plane anisotropy of the superconducting gap 8,10 . In tilted magnetic fields, the bias voltage dependence of the conductance maps is completely modified. The stripes observed at zero bias are no longer seen at higher bias (middle and right panels in Fig. 2). Instead, we observe broad dark-bright patterns along the tilt of the magnetic field.
Note that the images show several outgoing vortices, in addition to the stripes. Note also that this is very different from the observations made in refs 19,20 (at a single bias voltage and azimuthal angle and as a function of the value of the magnetic field) which show just stripes and no vortices. In our data, from the vortex density we obtain θ = 87.5°in  There is a small uncertainty in θ, which produces outgoing vortices and a pattern of stripes. Color scale of both figures is the same as in a. c The normalized tunneling conductance curves when crossing a stripe are shown in the main panel. Bottom curves are tunneling conductance data outside the stripe and top curves at the center of the stripe. In the lower right inset we show a zero bias conductance image with an outgoing vortex and the stripes. The red arrow provides the path giving the set of tunneling conductance curves of the main panel. Upper inset shows a zoom of the evolution of the conductance close to zero bias, with a dashed red line marking the splitting of the zero-bias peak. All data are taken at T = 150 mK COMMUNICATIONS PHYSICS | DOI: 10.1038/s42005-018-0028-1 Fig. 2 right column). These angles are within the possible misorientation with respect to a parallel magnetic field that can be obtained in our experiment (see Methods section). As we discuss in detail below, vortices might come out either perpendicular or at an angle to the surface. The shape of the vortices we observe in highly tilted fields, as well as their bias voltage dependence, is totally different than the shape of vortices in perpendicular fields. Their structures actually merge into the stripes when leaving the vortex centers, having the same two vs. three fold structure of rays that extend into the stripes far from vortex cores. They are also, as remarked above, considerably larger than vortices in perpendicular fields. Thus, we are observing vortices exiting the surface at an angle. Left panels show conductance maps with a field of 0.6 T parallel to the c axis. Middle and right panels show conductance maps at nearly parallel magnetic fields (θ ≈ 90°with an uncertainty of a few degrees). In the middle panels, the in-plane tilt is with the azimuthal angle along a crystalline direction (φ = 0) and in the right panels in between crystalline directions (φ = 30°). The in-plane direction of the tilt is shown by the red arrows in middle and right panels of a. Scale bars are of 30 nm in the left panels and 120 nm in middle and right panels. The color scale is given by the bar in a. The vortex density decreases for nearly parallel magnetic fields. The stripes at mostly visible at zero bias and follow the in-plane direction of the magnetic field in a. At high bias, the stripes are no longer present. Vortices present elongated shapes and are strongly in-plane asymmetric, with a clear yellow-blue pattern along the in-plane direction of the magnetic field (e, middle and right panels).
out that it is a formidable task to calculate the vortex core shape in tilted magnetic fields because one needs to obtain the DOS as a function of the position in three-dimensional space with a high spatial accuracy and then analyze the magnetic field and gap structure as a function of the tilt. This has never been achieved, to our knowledge. We start by building DOS maps of isolated vortices in an effective two-dimensional model (Supplementary Note 4) and calculate the DOS maps as a function of energy and the azimuthal angle of the magnetic field (Supplementary Note 5). We simulate a vortex tilted to the hexagonal plane by a two-fold distribution of the supercurrent along the tilt, using calculations of ref. 21 . We introduce in addition a shift q in the condensate momentum to account for the current distribution due to the inplane vs. out-of-plane crystalline anisotropy (Supplementary Note 5). In Fig. 3 we show the vortex-core shape from contours of constant order parameter around the vortex cores when tilting the magnetic field (a and b) with zero q and with nonzero q (c and d). With zero q, we see that the DOS pattern becomes elliptical. With a finite q we find that the DOS map of the vortex core develops a comet like shape, increasing considerably its size along the tilt of the magnetic field. The comet shape has a twofold structure perpendicular to the tilt when the tilt is along a crystalline axis (φ = 0°) and a triangular form when the tilt is at 30°to a crystalline axis (φ = 30°).
In Fig. 4 we show the energy dependence of the DOS with the field tilted along the nearest-neighbor direction of the vortices and perpendicular to this direction (middle and right columns, respectively). We see that the vortex core shape obtained from the DOS is highly energy dependent. The spatial extension of the DOS around a vortex core is considerably larger in tilted fields. This is consistent with the observed increase in the vortex core size in the tunneling conductance when comparing maps made at perpendicular field with maps at nearly parallel fields (Fig. 1a). For energies close to the gap edge, at 0.8 meV the spatial extension is considerably reduced and has more structure. Note that the model addresses a single vortex. Stripes result from vortex cores lying in-plane that reach the surface. Stripes of neighboring subsurface vortices extend to the surface and join each other. The shape of the order parameter extends along directions that change with the azimuthal angle φ (Fig. 3), indicating that along these directions states of neighboring vortices overlap and lead to the observed stripes. In Fig. 3c we observe that the vortex extends preferentially along its sides, leading to the two-fold structure seen in the left panel of Fig. 1b. In Fig. 3d (Fig. 1b, right panel), we observe instead that the vortex core extends along the tilt and its sides, leading to the additional stripe exiting from the center of vortex cores in the right panel of Fig. 1b. Note also that the model reproduces features of the bias voltage dependence of vortices exiting the surface at an angle. In particular, the increase of the DOS along the tilt for zero bias reverses for high bias. We can take in Fig. 4 the phase singularity (white point) as a reference to compare DOS maps for different bias. The tunneling conductance is high (bluish, high DOS) at low bias, but lower (yellowish, low DOS) at high bias along the direction of the magnetic field. Similarly, the vortices observed in Fig. 2 (middle and right columns) change their two-fold shape along the in-plane magnetic field.  In the left panels we show the DOS of a vortex with the field tilted along the a axis (φ = 0°, x-axis of the panels), calculated without taking into account surface currents. In the middle and right panels we show the result taking into account surface currents when the inplane tilt is φ = 0°(middle panels) and φ = 30°(right panels). The current is of qa = 0.008 in both cases. The direction of in-plane magnetic field and current is shown by the red arrows in a. The scale bar is shown in the left panel of (a) and is the same for all panels, of 30 nm size. The color map goes from yellow (minimum DOS) to the zero bias peak in the DOS in black. The white dots in each panel show the phase singularity point. The tilt strength is τ = 0.5 in all panels. Insets in middle and right panels of (a) display the core shape on a logarithmic gray scale together with the shape of the order parameter as an iso-contour of the DOS (red dashed lines) and an iso-contour of the DOS (red line).

Discussion
Vortices in tilted fields close to the surface are very different in 2H-NbSe 2 than in other materials. In the isotropic superconducting material β−Bi 2 Pd, vortices at the surface have the same shape whatever the in-plane direction of the magnetic field 22 . On the other hand, highly anisotropic materials like cuprates show two-dimensional pancake vortices, with the same shape at the surface than for perpendicular fields 23,24 . In both cases, the vortex core in tilted magnetic fields just shows the anisotropic properties in the plane, either because vortices bend close to the surface in β−Bi 2 Pd or because vortices are fully confined within the layers (pancake vortices) in the cuprates. In 2H-NbSe 2 the observed DOS patterns are spatially highly structured at all bias voltages, indicating a more intricate situation.
In order to explain this situation, let us consider the energetics of vortices close to the surface in tilted magnetic fields and how this influences the three different systems mentioned here. We have to consider the balance between the elastic energy cost associated with vortex bending, which favors straight vortices exiting at an angle to the surface, and the cost associated with establishing strongly distorted current loops around vortices close to the surface which favors circular current loops and vortices exiting perpendicular to the surface 24,25 .
The shear modulus c 44 is given by c 44 ≈ (1 + k 2 λ 2 ) −1 with k being the wavevector for the distortion 25 . Close to the surface, the wavevector relevant for this problem is k = 1/a where a is the bending radius. a is given by the minimum between the intervortex distance a 0 and the penetration depth λ 26 . For the magnetic fields considered here, a 0 is of about 90 nm, always smaller than λ, thus a ≈ a 0 . On the other hand, the current distribution surrounding the vortex is determined by the anisotropy of the upper critical field and is elliptical, with its short axis along the direction where the upper critical field is higher (and the coherence lengths smaller) 25,27,28 .
In β−Bi 2 Pd, vortices exiting the surface have a circular shape and bend close to the surface 22 . The shear modulus of the vortex lattice is isotropic and the lattice is sufficiently soft that it is energetically favorable to exit perpendicular to the surface for all tilts. The upper critical field in β−Bi 2 Pd is low (0.6 T) and the coherence length of about 20 nm, so that vortices are larger than in 2H-NbSe 2 and the cuprates and certainly round because the upper critical field is nearly isotropic 29,30 . In β−Bi 2 Pd, λ ≈ 100 nm, which increases the elastic energy associated with vortex bending close to the surface as compared to materials with larger λ 22,30 . However, this is balanced by the gain in establishing current loops parallel to the surface with bent vortices 21 .
In the cuprates, tilted fields create arrangements of crossing pancake and Josephson vortex lattices whose structure strongly depends on the temperature, magnetic field and tilt [31][32][33][34] . Lattices of pancake vortices slightly shifted between layers can appear under some conditions, leading in essence to tilted lines of vortices that are expected to exit at an angle to the surface 21,[35][36][37] . However, at least for materials in which the c-axis coherence length is below the interlayer distance, there is no reason for vortex bending as no Abrikosov vortices are established between layers. Current loops of Abrikosov vortices are parallel to the surface, and the balance is rather played by the interaction between Abrikosov and Josephson vortex lattices.
Here in 2H-NbSe 2 , the upper critical field anisotropy is of a factor of three and the out of plane coherence length is well above the c-axis lattice constant 15,38,39 . This creates elliptical Abrikosov vortices (there are no Josephson vortices) with the short axis perpendicular to the surface. This helps establishing current loops close to the surface, as the vortex is anyhow shorter along the direction perpendicular to the surface. Tilted vortices can thus exit at an angle to the surface, neither parallel nor perpendicular to the surface. This leads to an intricate surface current pattern that modifies the vortex core shape at all bias, as shown by the vortices imaged in our experiment.
On the other hand, the appearance of stripes is also unique to 2H-NbSe 2 . As we show above, this is due to subsurface vortex cores. The zero bias peak due to Caroli-de Gennes-Matricon states is indeed a characteristic feature of 2H-NbSe 2 and it extends indeed over distances far from the vortex core 7,8,10 . In our experiments, the distance between stripes from subsurface vortices is comparable to the intervortex distance expected for the applied in-plane magnetic field, with each subsurface vortex providing a different number of stripes, depending on the inplane direction of the tilt, as shown in Fig. 1b and as discussed in connection with Fig. 3.
The combined appearance of stripes and vortices could be also due, in principle, to crossing perpendicular vortex and parallel vortex lattices 40 . Within such a situation, however, it seems quite difficult to obtain patterns of stripes with distances that change when varying the in-plane direction of the tilt.
We can conclude that the differences between 2H-NbSe 2 , β−Bi 2 Pd and the cuprates are due to a totally different DOS pattern for parallel magnetic fields. The pattern results from subsurface vortex cores and provides access to the shape of inplane vortices and thus to the gap anisotropy in the out-of-plane three-dimensional part of the Fermi surface.
The shape of the superconducting gap along the c axis should influence the dependence of the core shape with the in-plane direction of the tilt. As we show above, the main asymmetric features observed in our experiment and their dependence with the azimuthal angle can be explained by a sixfold anisotropic superconducting gap in 2H-NbSe 2 . This shows that the sixfold gap anisotropy is spread over the whole Fermi surface.
It is useful to compare with the situation in MgB 2 , where vortex cores have been studied for magnetic fields always applied perpendicular to the surface, but with surfaces out-as well as inplane 41,42 . An elongated vortex core shape was expected for fields in-plane, due to the two-dimensional part of the Fermi surface 43 . This was however not observed. Vortices are round out-and inplane, with their shape being dominated by the three-dimensional part of the Fermi surface, both the vortex core shape in the DOS [41][42][43] as well as in the magnetic field pattern 44 . It was concluded that tunneling is dominated in all cases by the threedimensional part of the Fermi surface. Tunneling in 2H-NbSe 2 is also dominated by the three dimensional part of the Fermi surface, here due to the stronger contribution to tunneling of the Se orbitals 5 . However, the situation is quite different. Whereas the three-dimensional pocket amounts for the largest part of the DOS at the Fermi level in MgB 2 , the Se pocket in 2H-NbSe 2 is small and amounts just for 20% of the DOS 5 . The gap anisotropy in 2H-NbSe 2 was until now just associated to the Nb two dimensional sheets. As we show here, the sixfold gap anisotropy remains out-of-plane, which suggests that it is present over the whole Fermi surface.
Interband interactions are strong in 2H-NbSe 2 5 . Although there are evidences for a two gap behavior from tunneling spectroscopy, penetration depth and thermal conductivity 5,10,12,45,46 , the obtained spread in gap values is of a factor of 1.5 or at most 2. Our measurements suggest that interactions are such that the sixfold anisotropy in the shape of the superconducting gap is inherited from the Nb twodimensional bands over to the whole Fermi surface.
STM in tilted magnetic fields is thus a sensitive probe of the superconducting electronic DOS, particularly in threedimensional superconductors with strongly anisotropic and sharply defined vortex core bound states. The main advantage with respect to studying different crystalline surfaces is that the direction of the magnetic field can be varied at will on a single surface. Efforts, still under development, have been also devoted to other techniques. Quasiparticle interference has been extended out of the surface plane by triangulating electronic wavefunctions from in-depth impurities or studying inter and intraband scattering [47][48][49][50][51][52] . Macroscopic specific heat and thermal conductivity measurements in rotating fields have been also performed with the aim to obtain the anisotropy of the superconducting gap 53 . These provide the spatially averaged DOS and are influenced by scattering effects, whereas STM provides directly the spatially resolved electronic DOS, i.e. in and around vortex cores.

Methods
Model and calculations. We create a microscopic model with a sixfold anisotropic superconducting gap, by a triangular lattice and a tight-binding dispersion (see Supplementary Notes 4,5 for more information). To calculate the spatial dependence of the superconducting DOS we need a very large system size. In the lowenergy region the inter-level spacing is ΔE hv F Δk with Δk = 2π/(Na), where Na is the linear system size. In order to reach a resolution ΔE≲ 1 meV of the order of the superconducting gap, a total number of unit cells N 2 ≳200 000 is required. This sets the size of the Bogoliubov-de Gennes Hamiltonian to at least 400,000 × 400,000. Straight diagonalization is therefore not an option. Instead we use the method described in ref. 54 , the DOS is expanded on Chebyshev polynomials of the Hamiltonian and can thus be evaluated iteratively with low memory cost, even for very large systems. We use a finite lattice made of M concentric hexagons surrounding the site where the DOS is calculated. To reach the desired accuracy we perform the calculation with M = 1000, corresponding to a lattice of 3 003 001 sites. The Chebyshev expansion is truncated at order 4M and terminated using the Jackson kernel 55,56 . To obtain results in an applied magnetic field and as a function of position and energy we use a Lawrence-Doniach model following ref. 21 . We are not aware of an equivalent model for vortex latticeshence we consider only isolated vortices. In Supplementary Note 5 we solve the nonlinear model of ref. 21 to leading order in the tilt angle and get an analytical expression for the phase. The distortion is proportional to the tilt strength defined as τ ¼ ξ 2 ab =ð2ξ 2 c Þ tan 2 θ, with θ being the polar angle. The DOS maps and anisotropies of the order-parameter are determined self-consistently as a function of τ. We use tilt strengths of τ ≈ 0.5, which correspond to small θ ≈ 15°. We believe that this accounts well for our situation, with a three dimensional, strongly anisotropic superconductor. To further address this anisotropy, we introduce an additional phase in the order parameter, i.e., Φ(r) → Φ(r) − q ⋅ r, with q being a surface current. We take q up to about one third of the value corresponding to the critical current.
STM methods. We use a setup described in ref. 57 which consists of a dilution refrigerator with a STM thermally anchored to the mixing chamber reaching temperatures of 150 mK. The STM is located at the center of a three axis magnet. The three axis magnet consists of a long solenoid providing the z-axis magnetic field and a pair of split coils for x-and y-axis magnetic field components. We estimate the accuracy in the determination of the azimuthal and polar angles θ and φ at the magnetic fields used here (of order of a Tesla) to be about 4-5°. This uncertainty is composed of a slight misalignment of the surface of the sample with respect to the coil system and of remanent magnetic fields inside the coils, that can be in the range of the mT. We measure high quality 2H-NbSe 2 samples grown using the usual iodine vapor transport method. The in-plane tilt direction is corrected to obtain values for φ that are fixed to the atomically resolved hexagonal Se lattice obtained by STM imaging (see also Supplementary Fig. 5). We take tunneling conductance maps as a function of the bias voltage, as usual in vortex imaging using STM. Sample surface has been prepared by cleaving and the Au tip is prepared using a pad of Au as described in ref. 58 .